計測自動制御学会東北支部 第 205 回研究集会 (9 November 2002)
資料番号 205-2

非ガウス型状態空間モデルによる時系列データの分解— 化学プロセスの異常検出への応用 —Decomposition of Time Series by non-Gaussian State Space Model— Application to Fault Detection in a Chemical Plant —◯ 武田直人, 山下善之
◯ Naoto Takeda and Yoshiyuki Yamashita
東北大学大学院

Tohoku University
キーワード : プロセス監視 (Process Monitoring), 時系列解析 (Time Series Analysis),
非ガウス型モデル (non-Gaussian Model), 異常検出 (Fault Detection)

連絡先 : 〒 980-8579 仙台市青葉区荒巻字青葉 07 東北大学大学院工学研究科 化学工学専攻
e-mail: {ntakeda |yyama}@pse.che.tohoku.ac.jp

1. Introduction

Abstract
A novel framework of fault detection of non-To improve process quality, productivity andThissafety, fault detection is an important problem of

method is based on the following three steps: (1)practical interest in many chemical process plants.

decomposition of the observed signal into trend

Many algorithms for fault detection have been de-and the residual component through the use ofveloped by many researchers. Statistical change

the non-Gaussian state space model, (2) model-

detection is a relatively old research field, but re-ing of the residual component through the use ofcently interest has increased again with the active

an AR model with time varying coefficients, andresearch for data mining. The problem of detecting(3) checking the modeling error of the AR modelwith a sequential probability ratio test. The pro-a change in the parameters of a static or dynamic
stochastic system is well summarized in book2)

posed framework can isolate the effects of noises

with key mathematical background. It includes

and process changes from the changes of opera-

Shewhart control charts, geometric moving aver-tional modes. To demonstrate the method, sim-age charts, finite moving average charts, CUSUM-

ulated data from a simple process system with atype algorithms, GLR detector and Bayes-type

level control are successfully analyzed.

charts. Chiang et. al. also summarizes many ap-

stationary chemical process is proposed.

proaches for fault detection and isolation, includ-

2. Method

ing statistical methods such as SPC, PCA, PLS,
FDA, CVA and soft computing methods3) .

2.1 Decomposition and Fault Detection
The proposed framework of the fault detection
has the following three steps.

Although these methods are useful, real plant
includes various problems to apply these methods• Trend modeling: Modeling the trend of theand none of the method is sufficient for all the pur-normal observation signals by non-Gaussianpose. Before applying the usual methods, decom-

state space model.

position of operational data series into trend signal• Residual modeling: Modeling the residual ofand the other components can solve the problem inthe trend model by time varying coefficientmany situations.

AR model.
• Fault detection: Decision making by sequen-

Decomposition of process signal into several

tial probability ratio test on the residual of thecomponents has been studied with a primary fo-

AR model.

cus on the problem of sensor validation. Many approaches have been investigated including cluster
analysis, pattern recognition, modeling individual

2.2 Trend Modelsensors, modeling a process including sensors, andstatistical analysis10) . A linear Gaussian model istrend and stationary covariances is described here.

one of the powerful approaches for the decompo-

Kitagawa proposed a method to decompose an ob-

sition with the help of Kalman filter and AIC crite-

served time series into local polynomial trend, sea-ria. However it has the limitation not to apply theprocess with the abrupt changes of process opera-

sonal, globally stationary autoregressive and observation error components6) . In this model, eachtional mode.

component is characterized by an unknown vari-

The problem of modeling a time series withance — white noise perturbed difference equation
In this paper, decomposition of the observed

constraint.
In general, the state space representation for thetime series into trend and the other component isinvestigated under the changes of process opera-

observations yn istional mode. To realize this facility, the trend of

xn = F xn−1 + Gv n ,the observed time series is modeled by the non-yn = Hxn + wn ,Gaussian state space model using the Gaussian(1)

mixture approximation. Subsequently, the resid-

where xn is a state vector, F, G and H are m ×

ual of the trend model is modeled by time varyingm, m × `, and 1 × m matrices, respectively. For

coefficient AR model. Finally, the residual of thepossibly non-Gaussian white noises v n and wn ,AR model is tested by sequential probability ratio

the recursive formulas for obtaining the densities

test to detect faults.

of the one-step-ahead prediction and filtering are–2–

given as follows7) :The best choice of the model can be made byminimizing the value of AIC defined byp(xn |Y n−1 )
=

Z ∞

−∞

(2)
+2(number of free parameters).(8)p(yn |xn )p(xn |Y n−1 )p(xn |Y n ) =,p(yn |Y n−1 )(3)where p(yn |Y n−1 ) is obtained by

Under the Gaussian assumption, the likelihood
is represented byp(yn |Y n−1 )L(θ, x0|0 )

(4)

Z

=

AIC = −2 log(maximized likelihood)p(xn |xn−1 )p(xn−1 |Y n−1 )dxn−1 ,p(yn |xn )p(xn |Y n−1 )dxn .

=

N
Y

(
1

(2πr n|n−1 )− 2 exp

n=1−r 2n|n2r n|n−1)(9).

For the dimension value m = 1, 2 and 3, Equa-

where r n|n is the one-step-ahead output prediction

tion (1) represents the polynomial trend, station-

error (r n|n = yn − Hxn|n−1 ), r n|n−1 is its error

ary AR, and seasonal effects component models,covariance matrix and θ = (τ, σ).

respectively.

The polynomial trend component

For the evaluation of the AIC criteria, initialis represented as a k−th order stochastically per-value of the variance of the state vector is cal-turbed difference equation

culated by applying the estimation filter in re-∇k tn = v n ,(5)where ∇ denotes the difference operator definedby ∇t = tn − tn−1 .

the Pokal-Ribiere-Polyak method is used in thecase study experiment.

non-Gaussian model:

The corresponding matrices and state vector
components are:If the noise densities are non-Gaussian, Kitagawa developed an algorithm for implementing the(k = 1)

non-Gaussian filter by approximating each density

tn = tn−1 + v n ,(6)xn = tn , F = G = H = 1.

function based on a continuous piecewise linear
function7) . The method can be applied to lower
order systems, but it is not suitable for higher or-(k = 2)der state space models due to the huge amounttn = 2tn−1 − tn−2 + v n ,of computation that is requires. Another practicalxn = (tn , tn−1 )t ,
ÃFverse. As for the optimization search algorithm,

=

2 −1
1 0

!

Ã

, G = Ht =

(7)

1
0

!

.

way of performing the computation is the use ofa Gaussian-sum filter by using Gaussian mixture
approximation to the related densities.
Let ϕi be properly defined Gaussian density.

Gaussian model:
If the system noise and observation noise have
Gaussian densities, that is v n ∼The following approximations can be used.

N (0, τ 2 ) andwn ∼ N (0, σ 2 ), the value of p(xn |Y n−1 ) andp(xn |y n ) can be effectively estimated through theutilization of the Kalman filter1) .

–3–p(v n ) =p(wn ) =

Kv
X
i=1
Kw
X
j=1

αi ϕi (v n ),(10)

βj ϕj (wn ),(11)

p(xn |Y n−1 ) =p(xn |Y n ) =

Ln
X
k=1
Mn
X

γkn ϕk (xn |Yn 1 ),(12)δ`n ϕ` (xn |Yn ).

A k-th order difference model for the AR coefficients is defined by∇k anj(13)= vnj .

(17)`=1

The state space representation of the time vary-

Substituting these approximations into Equations (2) and (3), the following one-step-ahead prediction algorithm is derived8) :p(xn |Y n−1 )=ing coefficient AR model becomes

(14)

Kv M
n−1
X
X

xn = F xn−1 + Gv n ,(18)

zn = H n xn + wn ,(19)

where H n is the km dimensional vector.

γi`,n ϕi` (xn |Y n−1 ).

For example, if k = 1, the state space represen-

i=1 `=1

Here, ϕi` (xn |Y n−1 ) is the one-step-ahead predic-

tation of the model is as follows:

tor of xn under the assumption of Gaussian noise.
Therefore, ϕi` (xn |Y n−1 ) can be estimated by ordinary Kalman filter. Here, one-step-ahead prediction with non-Gaussian white noises can be estimated.
Similarly, the following filtering algorithm is

xn = [an1 , an2 , · · · , anm ]T ,(20)

F n = I m,

(21)

Gn = I m ,(22)

H n = [zn−1 , zn−2 , · · · , zn−m ],(23)

where I m is the unit matrix of m × m.

obtained:p(xn |Y n )
=

(15)

Kw X
Ln
X

2.4 Fault Detection
By using the above mentioned decompositions,δjk,n ϕjk (xn |Y n ).

j=1 k=1

several kinds of fault detection schematics could

In the implementation of the algorithm, reduc-be constructed. As an example of the application

tion of the number of Gaussian components is con-of the decomposition, a method using the residual

sidered. It was observed that a relatively smallof the time varying coefficient AR model is pro-number of Gaussian densities can approximate alarge class of distributions9) .

posed in this section.
If the system is normal, the modeling error of thetime varying coefficient AR model can be consid-

2.3Time varying AR coefficient model

ered to have a normal Gaussian distribution. Then,

The stationary model was discussed in the previous section. In this section, time varying AR coefficient model is considered, where the coefficienta changes gradually with time4) . The model is
nj

zn =

m
X

anj zn−j + wn ,j=1

where wn is a white noise N (0, σa2 ).

a fault can be detected by Wold’s sequential probability ratio test (SPRT) on this error signal.
The sequential probability ratio at time n is defined by
λn = λn−1 + log(16)

P (ηn |H1 ),P (ηn |H0 )(24)

where ηn is the modeling error at time n. H0 is the
hypothesis that the system is normal, and H1 is the–4–hypothesis that the system is faulty. P (ηn |H0 ) andh [m]P (ηn |H1 ) are the probability density functions forthe normal data and abnormal data, respectively.
A fault can be detected by monitoring the meanvalue and the variance of λn with the following

λn















 











 













Fin [m3/s]

ma= λn−1 + ma (ηn −
),
(25)
2
1
η2 1
),
= λn−1 − n ( 2 − 1) − log(va2(26)2 va2

where ma and va2 are the mean and variance of ηn

Fout [m3/s]

in abnormal situation.

3. Case Study
3.1





equations:
λn

















Simple Tank System




To demonstrate the ability of the method, it is



u [-]

applied to a simple tank system with level control.




Figure 1 shows the flow diagram of the simulated

system. The tank is modeled with an ordinary dif-

ferential equation based on the mass balance.
dh
1
cu √
= Fin −
h
dt
A
A

(27)

Fig. 2

Simulation run for normal operation

Observation noises are added to each observation variable (h, Fin , Fout and u). The valve characteristics are defined by a first order lag system.
In this process, the liquid level h is controlled byfor the integration of the system equation, and thea PI controller, parameters of which are tuned bysampling time of the control is 5.0 second.

Cohen-Coon method. Runge-Kutta method is used
Figure 2 shows an example of the simulation.

F in

In this simulation, initial set point of h was 0.8
m. It was changed to 0.3 m at t = 1000 seconds,hand changed again to 0.6 m at t = 2500 seconds.

u
F out

Initial value of the input flow rate Fin was 0.01
m3 /S, and the value was changed to 0.005 beforet = 1500 seconds. Other variables, Fout and u, areFig. 1The example process

calculated from the system equations.

–5–
 
 

 

  








  

Estimated trend for h


 
 



  








 

  

            









              



 



 



 

               

  



 



 



 

                 
 


 

  Estimated trend for Fout

Fig. 3 Estimated trend for Fin in Figure 2(non-Gaussian model)

3.2


 
 

 

  


Trend Modeling

For the data set of Figure 2, non-Gaussian modeling is applied in order to make the trend model of


 
 





 
 


 
 



  
the normal operation. Based on the algorithm described in the previous section, Figure 3(a) is obtained for the estimated trend component of Fin .
This result shows that the Gaussian mixture approximation for system noise can handle both the

Estimated trend for u

slow shift and the abrupt shift in the process signal. In the same figure, graph (b) plots the residuals, which is the differences between the estimatedtrend and the observation. Figure 4 shows esti-


 
 





 
 


 
 



  


mations for the other measurements. As shown in
the graph, the residuals do not include any abruptchanges. These graphs show that the non-Gaussianmodel can decompose the observation signal verywell.
For comparison, the result of the decomposition of the observation signal by Gaussian modelFig. 4 Estimated trend for h, Fout , u in Figure 2(non-Gaussian model)

is shown in Figure 5. The figure shows that, Gaussian model cannot model the abrupt changes, suchas the change at 1500 seconds and relatively large
residual remains at the point. Also, compared to–6–
 
 

 

  


h [m]








 
 



  






  



























 











 













 

Fin [m3/s]                



Fig. 5model)

Fout [m3/s]

Estimated trend for Fig. 2 (Gaussian

















u [-]








  









 


 

 

Fig. 7




 

  

      
 

Faulty operation (valve stuck)   using the AIC criterion. Figure 6(a) shows that
AR coefficient changes slowly, which corresponds

Fig. 6 Estimated AR coefficient and modelingerror for the residualwith the fact that the residual of the trend model

the non-Gaussian model, the trend component iseling error of the AR model, seems to be a white

not as smooth.

noise (Figure 6(b)).

3.3

4. Fault Detection

Residual Modelingdoes not have strange values. As a result, the mod-The residual of the trend model is modeled byNow, let us consider about the fault in anothertime varying coefficient AR model. Figure 6(a)simulation run (Figure 7). Let us suppose that ashows the estimated coefficient of the model. The

valve stuck at 2000 seconds. In this simulation,modeling error, which is the difference between

initial value of the set point of the level was 0.6

the AR estimated value and the residual, is shown

m. It was changed to 0.4 m at 1000 seconds, andin Figure 6(b). In this calculation, model param-changed again to 0.7 m at 2500 seconds. Althougheters are set to m = 1, k = 1, which is decided

the valve stuck at 2000 seconds, the level h vir-–7–

 





 
 
tion noise is included after 2000 seconds. With




 
 



  
the existence of changes in the set-point of liquid level or changes of the process disturbance, themethod detects the real fault with considerable ac

 



 



 

 

curacy. This result shows that the proposed framework provides an excellent tool for fault detectionof non-stationary chemical process.



  

              

5. Conclusions
In this paper, we have proposed a novel methodFig.
8 Estimated trend for u in Fig. 7(non-Gaussian model)
30

for the decomposition of process time series data
into the trend and the residual component. The
method is based on the non-Gaussian state space

Large Variance
Small Variance

25
Tresholdmodel and can deal with the changes of operating

SPRT

20
15

conditions. As an application of the decomposi-

10tion, a method of fault detection is shown by using

5
0
0

500

1000

30

1500

2000

2500

3000

Large Average
Small Average

the modeling error of the time varying coefficient
AR model, which models the decomposed residual

25
Treshold

SPRT

20

from the trend.

15
10

The method is applied to a simple tank system

5
0
0

500

1000

1500

2000

2500

3000

with liquid level control. The example process hastime [s]several operational changes, such as the set-point
Fig. 9changes; additionally, a valve stick error is suc-

Result of Fault Detection by SPRT

cessfully detected by the sequential probability ratually seems to be controlled normally until 2500

tio test of the modeling error. The example shows

seconds.

that the proposed method add a powerful tool for

The proposed framewrok is applied to this faulty
data. Figures 8 and 9 show the estimated trend
and SPRT value for the fault detection respectively.
As shown in the figure 9, the SPRT value of thevariance suddenly becomes large at 2000 seconds,

the decomposition of process trend. Our method is
also useful for the detection of fault in the process
having changes in operation mode. Future direction will likely include the use of the other decomposed components for fault detection.

which corresponds to the occurrence of the valve
stick. A sudden change in the variation of residual is also shown in the Figure 8. In normal data,u should include both the process noise and observation noise, but in this data, only the observa-–8–References

10) R. Luo, M. Misra and D. M. Himmelblau,

1) B.D.O.Anderson and J.B.Moore, Optimal
Filtering, Englewood Cliffs, NJ, Prentice-“Sensor Fault Detection via Multiscale Analysis and Dynamic PCA,” Ind. Eng. Chem.
Res.,, 38, 1489–1495, 1999.

Hall, 1979.
2) M. Basseville and I. V. Nikiforov. Detection of Abrupt Changes —Theory and Application, Englewood Cliffs, NJ, Prentice-Hall
Inc., 1993.
3) L. H. Chiang, E. L. Russell and R. D. Braats,
Fault Detection and Diagnosis in Industrial
Systems, London, Springer, 2001.
4) X. Q. Jiang and G. Kitagawa, “A time varying
coefficient vector AR modeling of nonstationary covariance time series,” Signal Processing, 33, 315–331, 1993.
5) G. Kitagawa, “Changing spectrum estimation”, Journal of Sound and Vibration, 89,
433–445, 1983.
6) G. Kitagawa and W. Gersch, “A Smoothness
Priors — State Space Monitoring of Time
Series With Trend and Seasonality,” Journalof the American Statistical Association, 79,
378–389, 1984.
7) G. Kitagawa, “Non-Gaussian state spacemodeling of nonstationary time series,” J.
Am. Statist. Assoc, 76, 1032–1064, 1987.
8) G. Kitagawa, “Non-Gaussian seasonal adjustment,” Computers Math. Applic, 18, 503–
514, 1989.
9) G. Kitagawa, “The two-filter formula forsmoothing and an implementation of theGaussian-sum smoother,” Ann. Inst. Statist.
Math, 46, 605–623, 1994.

–9–