Monitoring Nonstationary Dynamic Systems Using Cointegration and

Subscriber access provided by University of Pennsylvania Libraries

Article

Monitoring Nonstationary Dynamic Systems Using Cointegration And Common Trends Analysis Yuanling Lin, Uwe Kruger, and Qian Chen Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.7b00011 • Publication Date (Web): 10 Jul 2017 Downloaded from http://pubs.acs.org on July 16, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

MONITORING NONSTATIONARY DYNAMIC SYSTEMS USING COINTEGRATION AND COMMON TRENDS ANALYSIS

1,2

Yuanling Lin, 3Uwe Kruger, 1,2,••Qian Chen

1

2

College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, P.R. China

State Key Laboratory of Mechanics and Control of Mechanical Structures, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, P.R. China

3

Department of Biomedical Engineering, Jonsson Engineering Center, Rensselaer Polytechnic Institute, Troy NY12180-3590, U.S.A.

ABSTRACT This paper introduces a technique to monitor chemical processes that are driven by a set of serially correlated nonstationary and stationary factors. The approach relies on (i) identifying and separating the common stationary and nonstationary factors, (ii) modeling these factors using multivariate time series models and (iii) incorporating a compensation scheme to directly monitor these factors without being compromised by the effect of forecast recovery. Based on the residuals of the time series models, the technique yields two distinct test statistics to monitor both types of factors individually. Different from existing work, the paper highlights that the technique is sensitive to any fault condition and can extract and describe both, stationary and nonstationary trends. These benefits are illustrated by a simulation example and the application to an industrial semi-batch process describing nonstationary emptying and filling cycles.

KEYWORDS Common trend models, cointegration, latent stationary and nonstationary factors, forecast recovery, VARMA and VARIMA models.

•

Corresponding Author: Email: [email protected]

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 36

1 INTRODUCTION Industrial processes in the chemical and petrochemical industries often produce sizable data records that describe the current state of the process operation1,2. Measured variables include temperature and pressure readings, flow rate measurements, electrical readings for pumps and valves, concentration measurements through online analyzers etc. Through regulatory and advanced feedback control, the recorded variables, inherently, possess serial correlation3,4,5. In addition to that, the mean value of such variables may vary over time6,7,8,9. Typical causes of such nonstationary trends are changes in the operating condition, grade changes and variations in process feeds. Moreover, batch and semi-batch processes, for instance found in the pharmaceutical industry, typically produce nonstationary process variables. Based on their conceptual simplicity, multivariate statistical-based techniques have been developed and successfully applied to monitor such complex systems over the past few decades1. Core component technology includes principal component analysis, partial least squares and independent component analysis4,10,11,12. These methods assume that the process variables are random, have a joint multivariate distribution with a time-invariant mean and covariance matrix, and are not serially correlated. These assumptions, however, are rarely met in industrial practice. Although significant process has been made in handling serial correlation, adapting to slowly timevarying sample means and sample variances/covariances, and the presence of joint non-Gaussian distribution functions3,6,7,10,11,13,14, the issue of nonstationary process behavior has only been sporadically touched upon in the literature. Notable and recent exceptions include the work in Refs 8 and 15, which utilized stationary cointegration relationships16,17,18 for process monitoring.

Although successfully

employed to monitoring an industrial debutanizer process8, the proposed use of the stationary long-term equilibrium, manifested in the cointegration relationships, is insensitive to fault conditions that are orthogonal to the cointegrating vectors.

More precisely, existing work is either unable to handle

ACS Paragon - 2Plus - Environment

Page 3 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


nonstationary trends in the recorded variable set or relies on long-term equilibrium relationships that do not allow monitoring all latent stationary and nonstationary factors. The work in this article addresses the issue of monitoring stationary and nonstationary factors (i) by separating stationary and nonstationary components using the orthogonal Kasa decomposition19, (ii) by subsequently extracting latent nonstationary and stationary factors and (iii) by defining two test statistics to monitor these factors individually. The approach introduced in this article relies on common trend models20, which are developed on the basis of the cointegration theory. The core benefit of this approach is (i) that it is based on time-series models that are capable of capturing the serial correlation embedded within the recorded process variables and (ii) that it can extract latent sets of stationary and nonstationary factors from the recorded variable set. In sharp contrast, traditional multivariate statistical process control methods assume that the samples are drawn statistically independent from a joint multivariate distribution and have a time invariant joint distribution function1. Utilizing common trend models result in the definition of the following two multivariate models: •

a vector autoregressive integrated moving average, or VARIMA, model for describing the nonstationary trends that are embedded within the recorded variable set; and

•

a model that captures the equilibrium relationships by removing the nonstationary trends and is in the form of a vector autoregressive moving average, or VARMA, model.

Based on these models, the paper develops two test statistics to test the null hypothesis that the process is in-statistical-control against the alternative hypothesis that the process is out-of-statistical-control. These are: •

a Hotelling’s statistic based on the residuals of the VARIMA model, which includes a

residual compensation to prevent forecast recovery21,22, to monitor the latent nonstationary factors; and



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

•

Page 4 of 36

a Hotelling’s statistic based on the residuals of the VARMA model, which also includes a

residual compensation to extract the correct fault signature23, to monitor the latent stationary

factors. The next two subsections detail the assumptions imposed on the introduced monitoring scheme and summarize the contribution of this article.

1.1 Assumptions for monitoring approach In line with the standard assumptions for multivariate statistical process control, discussed above, we assume that the process operates at a specific operating point. Different to standard assumptions, however, the proposed framework does not assume that the process variables have a constant population mean and that each of the recorded data points are drawn independently. More precisely, the data structure of the recorded variables is composed of linear time-series models that may contain unit roots and are time invariant. This implies that the recorded variables can be serially correlated and can be stationary or nonstationary. If the process operates over a larger operating regime, e.g. possesses multiple operating points/grades, the identified models can be adapted to reflect changes in the variable interrelationships using a moving window approach for example. Although this is beyond the scope of this article, such adaptations are straightforward and well researched for multivariate statistical process control applications6,7,24. Examples that necessitate an adaptive implementation of the approach in this paper include chemical process systems that exhibit nonlinear steady-state relationships and transient changes between operating conditions/grade changes.

Moreover, if the process has operational

trajectories that are known or can be estimated, for example for batch processes, such trends can be removed prior to application of the introduced framework.


Page 5 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1.2 Contribution of this article Following from the preceding discussion, the difference between existing work on conventional multivariate statistical process control including the reported work on cointegration-based approaches (Refs 8 and 9) and the work presented herein is that existing work: •

does not consider detecting anomalous events in the orthogonal complementary space to that spanned by the cointegrating vectors, which Section 2 of this article exemplifies;

•

does not separate the recorded variable sets into a set of stationary and nonstationary trends;

•

does not take the inherent serial correlation of recorded process variables into consideration; and

•

does not address the problem of forecast recovery.

Conversely, the technique introduced in this article: •

simultaneously extracts stationary and nonstationary factors, given that the data space can be divided into two subspaces, i.e. the cointegrating space and its orthogonal complementary space – in contrast, the work in Ref 8 only focused on the detection of abnormal behaviors in the cointegrating space by ignoring the orthogonal complement space; and

•

removes serial correlations by applying the autoregressive filters and addresses the problem of forecast recovery by utilizing a compensation scheme – in contrast, the work in Ref 9, (i) neglects the fact that the presence of autocorrelation negatively impacts the Type I and II errors (this is because the assumption of independence between observations is violated, as described in Refs 4 and 5) and (ii) the differencing of nonstationary factors leads to the loss of information and the problem of forecast recovery.

2 LIMITATIONS OF EXISTING WORK ON COINTEGRATION This section presents a simulation example to illustrate that the work in Ref 8 may run into difficulties. With being an identity matrix of appropriate dimension, the simulation example has the following data

structure in the -domain:



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 36

() = () + ()

+ , + , () = ()

(1)

+ , () = ()

where:

0.2 −0.7 −0.6 0.2 0.5 0.1 0.3 0.1& 0.4 −0.7& % % = −0.7 0.4%, = 0.7 −0.3 0.8% 0.1 −0.6% −0.5 0.2 0.1% −0.5 $ 0.5 0.3 −0.4 0.6$ , = (

−1.2 0 0.2 0 ), =( ) 0 −0.3 , 0 −0.7

, = *

(2)

−0.4 0 0 0 −0.6 0+ 0 0 −0.5

The stationary random vectors ~-./, 0.250 and ~-./, 0.250 contain independently and identically distributed variables that have multivariate normal distributions of zero mean and a variance of 0.25. Moreover, the data structure in Eqs. (1) and (2) is defined in Subsection 3.1. From the above nonstationary process, a total of 3,000 samples were simulated. The first 1,500 samples served as reference samples to identify a cointegration model, whilst the remaining 1,500 samples were used to test the sensitivity of the cointegration model to multiple fault conditions. Utilizing the reference samples, the following cointegrating vectors were computed – discussed in Subsection 3.2:

1 = 23

3

0.146 −0.111 34 5 = −0.045 −0.218 0.008

0.133 0.017 0.142 0.086& % 0.077 −0.039% −0.066 0.074% 0.018 0.128$


(3)

Page 7 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


To demonstrate that utilizing the residuals constructed from the cointegration relationships among the five variables, manifested in Eq. (3), may yield insensitive monitoring statistics, the second half of the simulated data set was superimposed by a total of three distinct fault conditions, i.e. (7) = + Δ9 , : = 1,2,3:

Fault 1: A step-type fault condition that is superimposed on from data points 1501 to 3000, i.e. Δ; = (−8.445

3.106

1.614

−7.773 3.929).

Fault 2: A second step-type fault condition that is also superimposed on from data points 1551 to 3000, i.e. Δ; = (−0.688

Fault 3: A

sensor

Δ;4 = (0 0

bias

of

−5.181 9.871 0.290 6.305).

variable .? 7 (@)0 = >.1; ((@) + Δ9 )0 = >.?(@)0 + 1; Δ9 = / : = 1,2

(4)

[Figure 1 about here] For Eq. (4) to be equal to zero, the fault vector Δ9 must be orthogonal to the column space of 1, storing

the cointegrating vectors. More precisely, any fault condition that can be described as Δ = 1B C, where 1B is an orthogonal complement of 1, describing the fault direction, and C describes the fault magnitude,

is undetectable by the method proposed in Ref 8. Faults 1 and 2 were constructed to satisfy Eq. (4).



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 36

3 PRELIMINARIES This section presents some preliminaries required for developing the monitoring approach in Section 4. The next subsection describes the underlying data structure, whilst Subsections 3.2 and 3.3 summarize the basics of cointegration and common trend models, respectively. Finally, Section 3.4 discusses the issue of forecast recovery.

3.1 Data structure Figure 2 shows the data structure for modeling the recorded process variables, which comprises of a sum of linear combinations of a set of serially correlated stationary factors and linear combinations of a set of

serially correlated nonstationary factors. Defining as the time shift operator, the corresponding data

model has the following structure19,25,26:

() = () + () = 2 5 D

() E.

()

(5)

Here, the recorded variable set ∈ ℝH is a linear combination of a set of a random nonstationary factors, stored in the vector ∈ ℝIJ that is superimposed by a linear combination of a set of a random

stationary factors, stored in the vector ∈ ℝJ . Moreover, ∈ ℝH×IJ and ∈ ℝH×J are unknown

mixture matrices that define the contribution of the individual random vectors and upon the random vector . Without restriction of generality, we assume that L = M + M . [Figure 2 about here] The data models for the latent random vectors and in the -domain are of the following form: ACS Paragon - 8Plus - Environment

Page 9 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


()NO () = P () ()

(6)

() () = P () ()

and represent a nonstationary VARIMA model and a stationary VARMA model. In Eq. (6), () =

+ , + , + ⋯ + ,RIJ RIJ

and

P () = + P, + P, + ⋯ +

P,SIJ SIJ , () = + , + , + ⋯ + ,RJ RJ , P () = + P, + P, + ⋯ + P,SJ SJ , N = (1 − ) and T is the integration order. Moreover, the dimensions of the squared matrices ,9 , P,U and ,9 , P,U is M and M , respectively, and and are multivariable random

vectors that have a normal distributions with a mean of zero and the covariance matrices VWX ∈ ℝIJ ×IJ

and VWY ∈ ℝJ ×J , respectively. 3.2 Cointegration

Engle and Granger16 introduced the term cointegration to formulate the phenomenon that nonstationary

processes can have linear combinations that are stationary. More precisely, let ∈ ℝH be a random vector containing nonstationary variables, there is a vector 39 ∈ ℝH for which: Z9 = 3;9

(7)

is a stationary process and Z9 is the cointegration residual. It should be noted that for L nonstationary variables there are at most L − 1 different cointegrating vectors. If there are [, 1 ≤ [ ≤ L − 1, such

linearly independent vectors 39 (: = 1, ⋯ , [), is said to be cointegrated with a cointegrating rank of [, or there are [ cointegration relationships among the variable set .

cointegrating matrix 1 = 23

⋯ 3] 5 ∈ ℝH×] .


This allows defining the


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 36

Over the past few decades, econometricians and statisticians have developed a substantial body of work in developing the theory of cointegration and presenting numerous applications. This statistical

model framework allows estimating the cointegrating rank for and the cointegrating matrix. Specially,

two types of cointegration testing methods, the Engle-Granger two-step procedure16 and the Johansen procedure17, have been extensively used for constructing cointegration models.

3.3 Common trend models Given that cointegrated variables share common trends, Stock and Watson20 derived a common trends representation by establishing the connection between cointegration and common trends models. More

precisely, the random vector , which has the cointegrating rank [ , can be formed as the linear combinations of L − [ nonstationary trends, i.e., common trends, plus a set of stationary components: = 1;B ^ + _

(8)

where 1B ∈ ℝH×(H]) is the orthogonal complement of 1, ^ ∈ ℕH] is an random vector, storing a set of nonstationary variables that have no cointegration relationships and describe common trends, and _

represents the stationary component of .

Moreover, Escribano and Peña26 proved that , cointegrated with rank [ , is driven by L − [

common nonstationary and [ common stationary trends. This, in turn, implies that Eq. (8) is equivalent

to a common trends model, which can be represented in the form of a Kasa decomposition19: (@) = 1B 21;B 1B 5 1;B (@) + 121; 15 1; (@)

ACS Paragon Plus - 10 - Environment

(9)

Page 11 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where 1;B and 1; are the identified common nonstationary and stationary factors, respectively. For

the model in Eq. (9), 1B is referred to as the common trends loading matrix, and the vector 1;B

constitute the non-cointegrated part of . 3.4 Forecast recovery filter

Considering the VARIMA model in Eqs. (1) and (2), i.e. a(

1 0 −0.2 0 )+( ) b N () = () 0 1 0 0.7

and computing the difference between the @th sample and its prediction in the time domain yields: (@) = (@) + , (@ − 1) + , (@ − 2)

(10)

If there is a fault condition, assuming a step-type fault for simplicity, that arises at the @th index, we have: (7) (@) = (@) + ∆ = ( (@) + ∆ ) + , (@ − 1) + , (@ − 2)

(11)

Given that the diagonal matrix N contains the backward finite difference terms 1 − , i.e. N (7) (@ +

(7) d) = N( (@ + d) + ∆ ) = N (@ + d), it follows that (@ + d) → (@ + d) as d increases. This

phenomena is referred to as forecast recovery21,22 and implies that a step-type fault does not manifest

(7) itself in a step-type fault signature of the residual sequence (@ + d). This can be circumvented by

utilizing the following filter mechanism23:

(7) (7) (7) f (@) = ∑h9ij (:) f(@) = (7) (@) − f (@)


(12)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

More precisely, instead of using the samples (7) (@) , the filtered samples f(@) are used.

Page 12 of 36

This

mechanism was first used to remove the impact of utilizing VARMA filters upon estimating the fault magnitude in Ref 23.

4 MONITORING FRAMEWORK FOR NONSTATIONARY PROCESSES This section introduces a process monitoring scheme, which relies on the Kasa decomposition in Eq. (9) to construct common trends models.

The first step for constructing a monitoring model requires

estimating the number of nonstationary variables by using a unit root test27. This is followed by estimating the cointegrating matrix 1, using the cointegration test procedure17. Next, using Eq. (9)

allows extracting the stationary trends k = 1; and the nonstationary trends ^ = 1;B . Comparing Eqs.

(5) and (9), it follows that k ∝ , ^ ∝ , 121; 15 ∝ and 1B 21;B 1B 5 ∝ . Subsections 4.1 and 4.2

now describe how to utilize the serially correlated random vectors ^ and k to construct test statistics for testing the null hypothesis that the process is in-statistical-control against the alternative hypothesis that the process is out-of-statistical-control.

4.1 Monitoring stationary trends

Following from Eq. (6), the stationary vector contains serially correlated variables, as they are

essentially normally distributed filtered white noise sequence of zero mean. Consequently, k ∝ is

inherently serially correlated. Kruger et al. 4 and Xie et al. 5, however, proved that serial correlation among and between random variables can significantly alter the Type I and II errors of nonnegative quadratic test statistics, as it violates the assumption that the recorded samples have been drawn

independently, i.e. the @th sample is not independent of the (@ − 1)th sample. To address this issue, Refs 4 and 5 advocated identifying a time series model and design nonnegative test statistics on the basis of the model residuals.

Because of their simplicity, vector autoregressive or VAR models are

considered in this paper, as the estimation of the model parameters relate to the computationally efficient ACS Paragon Plus - 12 - Environment

Page 13 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Yule-Walker equation. With respect to Eq. (6), this implies that m = 0. The VAR model structure for

the random vector k is as follows:

J k(@) = ∑9i n9 k(@ − :) + o (@)

R

(13)

where p is the number of lagged terms, n9 , : = 1, ⋯ , p , are coefficient matrices and o denotes a random error vector that has a multivariate normal distributions with a mean of zero and the covariance

matrix Vo . It is important to note that the random vector o does not possess any serial correlation,

implying that its samples are independent and identically distributed. The next step is to construct a Hotelling’s statistic28 to monitoring the variation of the stationary components of : q o o = o; V

(14)

q o is the sample covariance of the random vector o , estimated from a total of r samples, i.e. Here, V ; q o = r ∑s V hi o (@)o (@). The critical value, or control limit, of the Hotelling’s statistic in Eq.

(14) for a significance of t is given by:

,u =

](s)(sv) s(s])

wu ([, r − [)

(15)

where wu ([, r − [) is the critical value of an w-distribution with [ and r − [ degrees of freedom and a

significance of t.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 36

4.2 Monitoring nonstationary trends

The extracted nonstationary trends ^ in the common trends model – Eq. (6) and Eq. (9) –, are assumed

to be integrated of order T = 1 and do not have a cointegration relationship. That is, its difference ∇^(@) ≡ ^(@) − ^(@ − 1) is stationary and can be represented in the form of a VAR model: IJ ∇^(@) = ∑Ui zU ∇^(@ − {) + 7 (@)

R

(16)

where p is the number of lagged terms, zU , { = 1, ⋯ , p , are coefficient matrices, and 7 denotes a random error vector that has a multivariate normal distributions with a mean of zero and the covariance matrix V7 . It should be noted that Eq. (16) is equal to Eq. (6) if m = 0. As before, the observations of

the random vector 7 are independent and identically distributed, and do not contain any serial

correlation. Finally, Eq. (16) can be reformulated in the form of a VARI model for the random vector ^: IJ ^(@) = ∑Ui |U ^(@ − {) + 7 (@)

R

v

(17)

where the coefficient matrices |U , { = 1, ⋯ , p , p + 1 are given by: ] + z if { = 1 |U = }zU − zU if { = 2, ⋯ , p , −zRIJ if { = p + 1

(18)

Since the random vector 7 does not possess any serial correlation and is stationary, the following Hotelling’s statistic can be constructed to monitor the nonstationary factors in : q7 7 = 7; V


(19)

Page 15 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


q7 is the sample covariance matrix of 7 , estimated from a total of r reference samples, i.e., Here, V

; q7 = r ∑s V hi 7 (@)7 (@). The critical value for rejecting the null hypothesis, or the control limit, of

the Hotelling’s statistic in Eq. (19) for a significance t is given by: ,u =

(H])(s)(sv) s(sHv])

wu (L − [, r − L + [)

(20)

where wu (L − [, r − L + [) is critical value of an w-distribution with L − [ and r − L + [ degrees of freedom and a significance of t.

5 SIMULATION EXAMPLE REVISITED This section revisits the simulation study in Section 2. Subsection 5.1 and 5.2 summarizes the identified monitoring models for stationary and nonstationary trends and show that the method introduced in Section 4 can detect the abnormal event in sharp contrast to existing work, which was unable to detect the first two fault conditions.

5.1 Identifying the dynamic monitoring models

The cointegrating matrix 1 is shown in Eq. (3) and, asymptotically, has an orthogonal column space to .

More precisely, this implies that = 1B is an orthogonal complement to 1. After computing the

stationary sequences of the random vector k, a VAR model was identified using the Yule-Walker

equation. Based on the BIC criterion29, the number of lagged terms for k was determined to be p = 1: 0.536 k(@) = *−0.027 −0.004

−0.039 −0.033 0.516 0.092+ k(@ − 1) + o (@), 0.062 0.479


(21)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 36

The estimated error covariance matrix of o is: 0.710 q o = *0.039 V 0.021

0.038 0.021 0.724 −0.077+. −0.077 0.766

(22)

Given how the data model for determining the 3000 samples was constructed in Eqs. (1) and (2), p = 1 is a correct estimate. To describe the nonstationary trends, a VARI model is identified for the nonstationary factors as described in Subsection 4.1. Again, applying the BIC criterion, the number of lagged terms for the random vector ^ order was determined to be 1. It follows from Eq. (1) and (2) that p = 1 is a correct estimate. The resultant VARI model is as follows: 0.338 0.659 0.048 ^(@) = ( ) ^(@ − 1) + ( −0.705 0.711 0.618

−0.049 ) ^(@ − 2) + 7 (@), 0.325

(23)

and the identified error covariance matrix for the random vector 7 is: q7 = ( 0.009 V −0.011

−0.011 ). 0.150

(24)

5.2 Application of identified monitoring scheme

Recall that Figure 1 confirmed that Faults 1 and 2 could not be detected using the Hotelling’s statistic. Given how the multiple fault conditions were injected, this result was expected, as Faults 1 and 2 were orthogonal to 1, i.e. Δ = 1B C, and hence k = 1; (7) = 1; ( + 1B C) = 1; .

In sharp contrast,

Figure 3 shows the statistic, which is based on the compensated VARI model residuals of the


Page 17 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


nonstationary factors, could detect the multiple fault conditions. In addition to that, each fault condition

had a different impact upon the statistic. Figures 4 (a) and 5 (a) show the theoretical impact of the

multiple faults to the residuals of the nonstationary and the stationary model, respectively.

To

demonstrate the necessity of applying the compensation filters, Figures 4 (b) and 5 (b) present the filtered residuals and outline that the fault signature and/or the magnitude of the fault impact are compromised. Conversely, utilizing the compensation terms allows correcting the negative impact of forecast recovery and effect on the fault magnitude for the nonstationary and stationary residuals sequences, which Figures 4 (c) and 5 (c) confirm, respectively. Finally, Figures 4 (d) and 5 (d) present the compensation term and confirm that they correctly describe the fault signature and magnitude for each fault condition. As expected, the compensated stationary residuals in Figure 5, extracted from the VAR model, only respond to the third fault condition, as Faults 1 and 2 are orthogonal to the cointegrating vectors. In contrast, the compensated stationary residuals in Figure 4, computed from the nonstationary VARI model, was sensitive each of the multiple fault conditions.

[Figures 3 to 5 about here]

6 APPLICATION STUDY TO AN INDUSTRIAL MELTER PROCESS Figure 6 presents a schematic diagram of the unit, which is part of a disposal unit for handling a powdery waste material that continuously enters the vessel. Raw glass is introduced in the form of frit every 120s. The binary composition is heated by four induction coils that are positioned around the vessel. The filling and heating produces an increasing liquid column until a desired height is reached at which the molten mixture is poured out through an exit funnel. After this, the next cycle of filling and heating begins. The recorded variable set includes readings from a total of eight temperature sensors, placed at various positions around the vessel, the power of each of the four induction coils and the voltage applied to the coil. A sampling time for each of the 13 variables is 5 minutes. A more detailed ACS Paragon Plus - 17 - Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 36

description of this process can be found in Refs 11 and 30. Subsections 6.1 and 6.2 describe how the two test statistics were constructed and applied to detect an abnormal process event, respectively.

[Figure 6 about here]


6.1 Identifying the dynamic monitoring models Figure 7 shows the recorded data sets, i.e. the 62.5 hours of reference data (750 samples) to identify the dynamic models and a second set of 25 hours (300 samples) that describe a fault condition. Analyzing the first two emptying and filling cycles showed that each of the 13 variables are nonstationary by applying the PP test27. Table 1 shows the results of the PP test statistic for which the critical value was

determined for a confidence of 99%. The next step was to determine the number of nonstationary

components using the Johansen test17, which revealed M = 5 followed by computing 1 ∈ ℝ4× and 1B ∈ ℝ4×= , which were determined to be: −0.286 −0.169 −0.026 0.573 0.412 0.094 1 = −0.198 −0.077 0.019 −0.454 −0.322 −0.002 −0.166

0.343 −0.067 0.127 −0.420 0.185 −0.135 −0.486 0.004 −0.272 0.140 −0.021 0.502 −0.225

0.453 0.219 −0.394 −0.207 0.037 0.120 0.042 −0.309 0.509 0.239 −0.309 −0.025 0.156

0.186 0.324 0.373 0.091 −0.122 −0.332 −0.495 0.175 0.148 −0.071 −0.127 −0.363 0.204

−0.338 0.396 0.183 −0.128 0.166 0.225 0.101 0.231 −0.161 0.092 −0.702 0.066 −0.061

0.325 −0.517 −0.127 0.005 −0.340 0.155 0.114 0.551 −0.139 0.112 −0.286 −0.009 0.205


−0.010 −0.195 0.113 −0.060 0.399 −0.720 0.335 0.219 0.276 0.124 −0.123 −0.020 −0.031

−0.057 0.359& −0.687% % 0.075 % 0.173% −0.141% −0.202% 0.489% −0.132% −0.023% 0.189% −0.060% −0.055$

(25)

Page 19 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


−0.180 −0.220 −0.344 −0.317 −0.245 −0.319 1B = −0.072 −0.342 −0.269 −0.309 −0.344 −0.307 −0.190

−0.430 −0.400 −0.048 −0.145 0.334 0.194 −0.503 0.022 0.249 0.016 −0.166 −0.192 0.323

−0.228 −0.116 −0.056 0.171 −0.349 0.024 −0.174 0.122 0.471 0.179 −0.074 0.166 −0.669

0.139 0.016 0.147 −0.517 0.122 0.274 0.106 0.301 0.287 −0.531 0.170 −0.228 −0.236

−0.225 0.034& −0.099% % −0.031 % −0.212% −0.141% −0.020% 0.039% 0.259% −0.491% −0.077% 0.630% 0.402$

[Table 1 about here]

The next step was to construct the VAR model for the 8 stationary and the VARI model for the 5 nonstationary components. Applying the BIC criterion29, the number of lagged term for the VAR and

VARI models were estimated to be p = 4 and p = 4, respectively. The appendix lists the resultant coefficient matrices for both models. The estimated error covariance matrices for o and 7 were:

and

0.269 −0.050 −0.024 q o = 0.002 V −0.012 −0.039 −0.011 0.012

−0.050 0.327 −0.054 −0.063 −0.006 0.000 0.016 −0.033

−0.024 −0.054 0.196 0.002 0.021 −0.017 −0.011 −0.017

0.002 −0.063 0.002 0.179 −0.014 0.037 −0.020 0.030

0.012 −0.003 −0.003 0.009 q7 = −0.013 V 0.006 −0.021 0.012 −0.020 0.011

−0.012 −0.001 0.021 −0.014 0.246 −0.030 −0.003 0.034

−0.013 0.006 0.032 0.037 0.024

−0.039 0.000 −0.017 0.037 −0.030 0.331 0.017 −0.027

−0.011 0.016 −0.011 −0.020 −0.003 0.017 0.201 −0.004

−0.021 −0.020 0.012 0.011& % 0.037 0.024%, 0.111 0.060% 0.060 0.114$

respectively.


0.012 −0.033& % −0.017% 0.030% , 0.034% −0.027% −0.004% 0.470$

(26)

(27)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 36

6.2 Application of identified monitoring scheme The recorded fault condition here was a developing crack in the melter vessel that was recorded over a

period of 25 hours (300 samples). Figures 8 (a) and 8 (b) show the and statistics describing the

variation of the stationary and nonstationary trends in , respectively. As before, the residuals were

generated by the identified VAR and VARI models. The residuals were also compensated using the approach discussed in Subsection 3.4 to avoid the undesired effect of forecast recovery21,22 and an alteration of estimated fault magnitude23. With respect to the control limits, which were determined for

a significance of t = 0.01, both figures highlight that the fault condition was detected around 40 samples into the data set.

For comparison, a PCA-based monitoring model for the original recorded data set was also applied. Figures 8 (c) and 8 (d) show that the PCA-based monitoring statistics only generate sporadic violations that do not reveal the presence of abnormal behaviors at the earlier stage of the developing crack.

Moreover, Figure 8 (e) shows the resulting Hotelling’s statistic based on a cointegration model identified using the technique discussed in Ref 8. It is important to note that the recorded data show strong nonstationary trends, as confirmed by the analysis in Table 1. The application of PCA, however, requires that the process variables are stationary. In addition to that, Figure 8 (e) shows that the technique in Ref 8 only shows a sporadic response to the fault condition. Conversely, the method introduced in Section 4 is tailored to handle nonstationary trends and is, hence, more suitable for the application to this semi-batch process, as exemplified by its increased sensitivity in detecting the fault condition.



Page 21 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


7 CONCLUDING SUMMARY This paper has developed a technique for monitoring chemical systems that are driven by dynamic nonstationary as well as stationary factors, which is a common phenomenon in industrial practice. Such behavior results from controller feedback, typically through standard regulatory PID and advanced supervisory controllers, changes in throughput or individual feeds and the presence of other measured or unmeasured disturbances. The method introduced in this article relies on extracting the common trends, which can be identified using the standard cointegration framework. This allows dividing the recorded variable set into a set of latent stationary and nonstationary factors. To model these serially correlated latent variable sets, the method utilizes VAR and VARI models to produce stationary and serially uncorrelated residuals that can be used to construct test statistics for process monitoring. To demonstrate that existing work, proposed to deal with nonstationary trends in recorded variable sets, is not capable to detect generic fault conditions, the article summarizes the application to a simulation example.

Existing work was not able to detect specifically designed fault conditions,

whereas the technique introduced in this paper could detect these abnormal events. To present a more realistic application, the paper summarizes the analysis of recorded data that stem from a semi-batch process. This process shows strong nonstationary trends, which result from emptying and filling cycles. The recorded set from this glass melter process also show the presence of a developing crack in the melter vessel, which the proposed method in this article was able to detect in its early stages. In contrast, a conventional PCA-based approach was only able to detect the latter stages of the developing crack but was not sensitive to the early incipient stages of this fault condition. For future work, we suggest to expand the work of this article to handle nonlinear steady-state relationships, for example, by utilizing nonlinear subspace identification31,32 instead of the VARMA and VARIMA models or by adapting the linear version using a moving window approach, which has received significant attention for conventional multivariate statistical process control applications6,7,24.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 36

ACHKNOWLEDGEMENT The authors gratefully acknowledge the support from a project funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD).

APPENDIX The appendix describes the dynamic models of the stationary and nonstationary factors for the industrial melter process in Section 6.

Dynamic VAR model of stationary factors

( + n + n + n4 4 + n )k() = o ()

where:

−0.833 0.004 0.012 0.067 n = −0.010 0.122 −0.105 −0.122

−0.082 −0.738 0.090 0.010 0.038 −0.104 −0.033 0.130

0.083 −0.013 −0.682 −0.047 0.071 −0.168 −0.019 −0.173

−0.131 −0.063 −0.014 −0.665 −0.040 −0.108 −0.015 0.014

−0.103 0.061 −0.017 −0.101 −0.568 0.062 0.006 −0.085

−0.030 0.070 −0.077 −0.039 0.097 −0.470 0.050 0.120

−0.009 0.084 −0.004 −0.010 −0.067 −0.016 −0.578 0.056

0.064 0.027& −0.070% % −0.043% 0.064% −0.040% −0.008% −0.435$

0.009 −0.068 0.059 −0.055 n = −0.021 −0.187 0.072 −0.007

−0.144 −0.056 0.018 0.044 −0.006 0.211 −0.008 −0.121

0.0097 0.0954 0.1635 0.0865 0.0616 0.0611 −0.1214 −0.0834

−0.052 0.194 −0.065 −0.161 −0.006 −0.036 −0.073 0.053

0.043 0.047 −0.018 0.005 −0.205 −0.088 −0.013 0.064

−0.046 0.117 −0.025 0.024 −0.025 −0.137 0.021 −0.052

0.003 −0.105 0.036 −0.010 0.036 −0.077 −0.197 −0.162

−0.039 −0.030& % 0.046% 0.003% 0.021% 0.070% 0.037% −0.090$


Page 23 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


−0.019 0.099 −0.031 −0.002 n4 = 0.058 0.035 −0.026 0.114

0.014 0.061 −0.055 −0.038 0.045 −0.047 0.017 −0.010

−0.128 0.038 0.011 −0.146 −0.034 0.120 0.103 −0.013 −0.034 −0.059 0.089 0.032 −0.194 0.020 0.203 0.093

0.029 −0.079 0.040 0.005 0.012 −0.061 0.007 0.041

0.041 0.021 0.048 −0.037 0.011 −0.058 −0.022 −0.054

0.094 −0.023 −0.052 −0.007 −0.038 0.002 −0.074 −0.014

−0.019 −0.014& −0.013% % −0.020% −0.077% −0.009% −0.046% −0.135$

0.014 −0.067 0.023 0.042 n = −0.052 0.073 −0.009 −0.054

0.147 −0.067 0.006 −0.047 −0.047 −0.074 0.064 0.099

0.044 0.044 0.042 −0.033 0.038 0.096 0.040 −0.092

0.043 −0.026 −0.065 0.082 −0.095 0.087 0.033 −0.013

0.036 −0.077 0.013 −0.003 0.039 −0.141 −0.004 −0.007

−0.031 0.077 −0.017 −0.023 0.087 0.046 −0.029 0.123

−0.055 0.057& % 0.018% 0.004% −0.008% 0.050% −0.047% −0.099$

0.093 −0.024 −0.018 −0.012 0.002 0.090 0.038 −0.071

Dynamic VARI model of nonstationary factors

( + | + | + |4 4 + | + |= = )^() = 7 ()

where

−1.456 0.213 | = 0.631 0.472 −0.302

−0.214 −1.183 0.215 0.837 0.877

0.267 0.038 −0.229 0.011 | = −0.538 −0.034 −1.351 0.072 0.310 −0.130 0.215 0.082 0.102 0.095 |4 = −0.223 −0.108 0.964 −0.546 0.523 −0.499

0.180 −0.024 0.045 −0.086 −0.036 −0.057& % −1.019 0.061 −0.172% −0.382 −0.678 −0.388% −0.314 0.054 −0.907$ −0.244 0.127 0.197 0.679 0.554

0.054 −0.058 −0.032 −0.287 −0.243

0.047 −0.100 −0.004 0.099& % −0.014 0.162% −0.084 0.412% −0.030 0.220$

−0.030 0.055 0.012 −0.021& % −0.074 −0.005% −0.271 0.021% −0.115 −0.171$



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

0.078 0.115 −0.107 0.065 | = −0.153 −0.023 −0.510 −0.445 −0.721 −0.400 |= =

−0.090 −0.072 0.085 0.037 0.273 −0.087 0.552 0.185 0.433 0.183

−0.005 0.036 −0.180 −0.022 −0.039

−0.016 −0.029 0.114 0.092 −0.040

0.003 −0.032 −0.004 −0.021& % 0.029 0.016% 0.089 0.080% 0.065 −0.006$ 0.016 0.007 0.016 0.002& % −0.055 0.007% −0.031 −0.088% 0.001 −0.044$


Page 24 of 36

Page 25 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


REFERENCES 1

Kruger, U.; Xie, L. Statistical Monitoring of Complex Multivariate Processes; John Wiley & Sons Ltd: West Sussex, UK, 2012.

2

Ge, Z.; Song, Z. Multivariate Statistical Process Control: Process Monitoring Methods and Applications; Springer: London, UK, 2012.

3

Ku, W.; Storer, R. H.; Georgakis, C. Disturbance Rejection and Isolation by Dynamic Principal Component Analysis. Chemometrics and Intelligent Laboratory Systems 1995, 30(1), 179-196.

4

Kruger, U.; Zhou, Y.; Irwin, G. W. Improved Principal Component Monitoring of Large-scale Processes. Journal of Process Control 2004, 14(8), 879-888.

5

Xie, L.; Kruger, U.; Lieftucht, D.; Littler, T.; Chen, Q.; Wang, S. Q. Statistical Monitoring of Dynamic Multivariate Processes – Part 1. Modeling Autocorrelation and Cross-correlation. Industrial & Engineering Chemistry Research 2006, 45(5), 1659-1676.

6

Wang, X.; Kruger, U.; Lennox, B. Recursive Partial Least Squares Algorithms for Monitoring Complex Industrial Processes. Control Engineering Practice 2003, 11(6), 613-632.

7

Wang, X.; Kruger, U.; Irwin, G. W. Process Monitoring Approach Using Fast Moving Window PCA. Industrial & Engineering Chemistry Research 2005, 44(15), 5691-5702.

8

Chen, Q.; Kruger, U.; Leung, A. Y. Cointegration Testing Method for Monitoring Nonstationary Processes. Industrial & Engineering Chemistry Research 2009, 48(7), 3533-3543.

9

Li, G.; Qin, S. J.; Yuan, T. Nonstationarity and cointegration tests for fault detection of dynamic processes. IFAC Proceedings Volumes. 2014; pp 10616-10621.

10 Kano, M.; Tanaka, S.; Hasebe, S.; Hashimoto, I.; Ohno, H. Monitoring Independent Components for Fault Detection. AIChE Journal 2003, 49(4), 969-976. 11 Liu, X.; Xie, L.; Kruger, U.; Litter, T.; Wang, S. Q. Statistical-based Monitoring of Multivariate Mon-Gaussian Systems. AIChE Journal 2008, 54(9), 2379-2391.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 36

12 Gunther, J. C.; Conner, J. S.; Seborg, D. E. Process Monitoring and Quality Variable Prediction Utilizing PLS in Industrial Fed-batch Cell Culture. Journal of Process Control 2009, 19(5), 914-921. 13 Treasure, R. J.; Kruger, U.; Cooper, J. E. Dynamic Multivariate Statistical Process Control Using Subspace Identification. Journal of Process Control 2004, 14(3), 279-292. 14 Li, W.; Qin, S. J. Consistent Dynamic PCA Based on Errors-in-variables Subspace Identification. Journal of Process Control 2001, 11(6), 661-678. 15 Shi, H.; Chen, Q.; Lin, Y. Fault Diagnosis of Non-stationary Engineering System Using Cointegration Coefficients Matrix. Journal of Vibration and Shock 2015, 34(1), 145-150. 16 Engle, R. F.; Granger, C. W. Co-integration and Error Correction: Representation, Estimation, and Testing. Econometrica: Journal of the Econometric Society 1987, 55(2), 251-276. 17 Johansen, S. Likelihood-based Inference in Cointegrated Vector Autoregressive Models; Oxford University Press: Oxford, 1995. 18 Pfaff, B. Analysis of Integrated and Cointegrated Time Series with R, 2nd ed.; Springer Science & Business Media: New York, 2008. 19 Kasa, K. Common Stochastic Trends in International Stock Markets. Journal of Monetary Economics 1992, 29(1), 95-124. 20 Stock, J. H.; Watson, M. W. Testing for Common Trends. Journal of the American Statistical Association 1988, 83(404), 1097-1107. 21 Superville, C. R.; Adams, B. M. An Evaluation of Forecast-based Quality Control Schemes. Communications in Statistics – Simulations and Computations 1994, 23(3), 645-661. 22 Apley, D. W.; Shi, J. The GLRT for Statistical Process Control of Autocorrelated Processes. IIE Transactions 1999, 31(12), 1123-1134. 23 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. Industrial & Engineering Chemistry Research 2006, 45(5), 1677-1688. ACS Paragon Plus - 26 - Environment

Page 27 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


24 Li, W.; Yue, H. H.; Valle-Cervantes, S.; Qin, S. J. Recursive PCA for adaptive process monitoring. Journal of Process Control 2000, 10(5), 471-486. 25 Gonzalo, J; Granger, C. Estimation of Common Long-memory Components in Cointegrated Systems. Journal of Business & Economic Statistics 1995, 13(1), 27-35. 26 Escribano, A.; Peña, D. Cointegration and Common Factors. Journal of Time Series Analysis 1994, 15(6), 577-586. 27 Phillips, P. C.; Perron, P. Testing for A Unit Root in Time Series Regression. Biometrika 1988, 75(2), 335-346. 28 Hotelling, H. The Generalization of Student’s Ratio. Annals of Mathematical Statistics 1931, 2, 360-378. 29 Hannan, E. J. The Estimation of the Order of An ARMA Process. The Annals of Statistics 1980, 8, 1071-1081. 30 Zeng, J.; Kruger, U.; Geluk, J.; Wang, X.; Xie, L. Detecting Abnormal Situations Using the Kullback-Leibler Divergence. Automatica 2014, 50(11), 2777-2786. 31 Lovera, M.; Gustafsson, T; Verhaegen, M. Recursive Subspace Identification of Linear and Nonlinear Wiener State-space Models. Automatica 2000, 36(11), 1639–1650. 32 Goethals, I.; Pelckmans, K.; Suykens, J. A. K.; De Moor, B. Subspace Identification of Hammerstein Systems Using Least Squares Support Vector Machines. IEEE Transactions on Automatic Control 2005, 50(10), 1509-1519.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 36

TABLE OF FIGURES Figure 1: (a) Simulated nonstationary sequences for simulation example. (b) Hotelling’s T statistic based on method proposed in Ref 8.

29

Figure 2: Graphical illustration of data structure used in this article.

29

Figure 3: T statistic based on compensated VARI model residual of nonstationary factors (simulation

example).

30

Figure 4: (a) Superimposed step-type fault signatures in cointegrating complement space. (b) Original VARI model residual of nonstationary factors (simulation example). (c) Compensated VARI model residuals of nonstationary trends (model residuals).

(d) Identified step-type fault

signatures – compensation terms in Eq. (12) (simulation example).

30

Figure 5: (a) Superimposed step-type fault signatures in cointegrating space. (b) Original VAR model residual of stationary factors (simulation example). (c) Compensated VAR model residuals of stationary trends (model residuals). (d) Identified step-type fault signatures – compensation terms in Eq. (12) (simulation example). Figure 6: Schematic diagram of glass melter process.

31 31

Figure 7: (a) Reference data to produce a monitoring model and (b) data showing the developing fault condition in glass melter process.

32

Figure 8: (a) T statistic based on compensated VAR model residuals of stationary factors (glass melter

process). (b) T statistic based on compensated VARI model residuals of nonstationary

factors (glass melter process). (c) Hotelling’s T statistic based on principal components (glass

melter process). (d) SPE statistic based on PCA model residuals (glass melter process). (e) Hotelling’s T statistic based on cointegration model (glass melter process).


33

Page 29 of 36

FIGURES

Nonstationary sequences

10 10

3

2

2

10

4

T

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


10 10

0

500

1000

1500

2000

2500

3000

10

1

0

-1

0

500

1000

Sample

Sample

(a)

(b)

1500

2000

Figure 1: (a) Simulated nonstationary sequences for simulation example. (b) Hotelling’s statistic based on method proposed in Ref 8.

factors nonstationary t1

u1

y1

t2

u2

y2

...

... tns

At

observations cointegrated

stationary

+

... us

Bu

=

yN

y

Figure 2: Graphical illustration of data structure used in this article.



10

2

2

10

4

T

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 36

10

10

0

-2

0

500

1000

1500

2000

2500

3000

Sample Figure 3: C statistic based on compensated VARI model residual of nonstationary factors (simulation example).

Step-type fault signatures

Filtered residuals 6.7

0 0.7

0

0

500

1000

1500

2000

2500

3000

0

0

0

0

0

500

1000

1500

2000

Sample

Sample

(a) Filtered and compensated residuals

(b) Compensation terms

2500

6.7

0

6.7

0 0.7

0

0

3000

500

1000

1500

2000

2500

3000

0.7

0

0

500

Sample

1000

1500

2000

2500

3000

Sample

(c)

(d)

Figure 4: (a) Superimposed step-type fault signatures in cointegrating complement space. (b) Original VARI model residual of nonstationary factors (simulation example). (c) Compensated VARI model residuals of nonstationary trends (model residuals). (d) Identified step-type fault signatures – compensation terms in Eq. (12) (simulation example).


Page 31 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Step-type fault signatures

Filtered residuals 2.6 3.2

0 0

5.9

0 0

-2.9

57.4

30.8

0

0 0

500

1000

1500

2000

2500

3000

0

500

1000

1500

2000

Sample

Sample

(a) Filtered and compensated residuals

(b) Compensation terms 2.6 3.2

0 0

2500

2.6 3.2

0 0

57.4

0 0

3000

57.4

0 500

1000

1500

2000

2500

3000

0

500

Sample

1000

1500

2000

Sample

(c)

(d)

Figure 5: (a) Superimposed step-type fault signatures in cointegrating space. (b) Original VAR model residual of stationary factors (simulation example). (c) Compensated VAR model residuals of stationary trends (model residuals). (d) Identified step-type fault signatures – compensation terms in Eq. (12) (simulation example).

Figure 6: Schematic diagram of glass melter process.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Recoreded time sequences

Recoreded time sequences

M1

M1

M2

M2

M3

M3

M4

M4

M5

M5

M6

M6

M7

M7

M8

M8

M9

M9

M10

M10

M11

M11

M12

M12

M13

M13 0

150

300

450

600

Page 32 of 36

750

0

Sample Number

100

200

300

Sample Number

(a)

(b)

Figure 7: (a) Reference data to produce a monitoring model and (b) data showing the developing fault condition in glass melter process.


Page 33 of 36

10

6

10

4

4

2

T

T

2

10

2

10

10

10

2

10

0

0

50

100

150

200

250

10

300

0

-2

0

50

100

150

Sample

200

250

300

200

250

300

Sample

(a)

(b)

60

10

3

50 10

2

SPE

2

40

T

30

10

1

20 10

0

10 0

0

50

100

150

200

250

300

10

-1

0

50

Sample

100

150

Sample

(c)

(d) 10 10

4

3

2

10

5

T

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


10 10 10

2

1

0

0

100

200

300

Sample (e)

Figure 8: (a) C statistic based on compensated VAR model residuals of stationary factors (glass melter process). (b) C statistic based on compensated VARI model residuals of nonstationary factors (glass melter process). (c) Hotelling’s statistic based on principal components (glass melter process). (d) SPE statistic based on PCA model residuals (glass melter process). (e) Hotelling’s statistic based on cointegration model (glass melter process).



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 36

LIST OF TABLES

Table 1: PP test results for recorded variables set (glass melter process)


35

Page 35 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TABLES Table 1: PP test results for recorded variables set (glass melter process)

critical value

test result

9

statistic -1.9926

-4.0090

nonstationary

7

-9.2068

-4.0093

stationary

9

-2.1424

-4.0090

nonstationary

Monitoring Nonstationary Dynamic Systems Using Cointegration and

Recommend Documents