Design of multivariate alarm systems based on online calculation of variational directions

chemical engineering research and design 1 2 2 ( 2 0 1 7 ) 11–21

Contents lists available at ScienceDirect

Chemical Engineering Research and Design journal homepage: www.elsevier.com/locate/cherd

Design of multivariate alarm systems based on online calculation of variational directions夽 Kuang Chen a , Jiandong Wang b,∗ a b

College of Engineering, Peking University, Beijing, China College of Electrical Engineering and Automation, Shandong University of Science and Technology, Qingdao, China

a r t i c l e

i n f o

a b s t r a c t

Article history:

Alarm systems are critically important for safety and efficiency of industrial plants, but are

Received 26 August 2016

severely suffering from alarm overloading. This paper proposes a method to design a multi-

Received in revised form 10 March

variate alarm system based on variational directions of involved process variables, in order

2017

to alleviate the severity of alarm overloading. By using adaptive time scales, time gradients

Accepted 4 April 2017

of signals are extracted to calculate variational directions in an online manner. Adaptive

Available online 13 April 2017

time scales are determined from a function between time scales and volatilities of process variables. Alarms arise at the moment that a nominal relationship of variational directions

Keywords:

among process variables is invalidated. Numerical and industrial examples illustrate the

Industrial alarm systems

effectiveness of the proposed method.

Variational directions

© 2017 Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved.

Time gradients

1.

Introduction

Modern industrial plants are usually monitored by computerized alarm systems, which play a critically important role for safe and efficient operations of industrial plants (ISA, 2009; EEMUA, 2013). However, many existing industrial alarm systems are suffering from alarm overloading (Bransby and Jenkinson, 1998; Rothenberg, 2009; Hollifield et al., 2011), i.e., there are too many alarms appearing in a short time period to be promptly handled by plant operators. The phenomenon of alarm overloading is clearly revealed from Table 1 based on a study of 39 industrial plants from oil and gas, petrochemical, and power industries (Rothenberg, 2009). The performance metrics such average alarms per day for these industrial plants are much worse than the recommended performance benchmarks (the second column in Table 1) from the Engineering Equipment and Materials Users’ Association (EEMUA). Among the occurred alarms, there are many nuisance alarms that are not associated with any actual abnormalities. Due to “cry wolf” effects caused by nuisance

alarms, operators may have little trust on alarm systems and fail in noticing informative alarms that are associated with abnormal conditions. One main reason of nuisance alarms is that the most common way of arising alarms in process industries is to compare the measurements of a single process variable with a high or low alarm threshold, and related process variables are not taken into account in designing alarm systems (Wang et al., 2016). Thus, the design of multivariate alarm systems is an important research topic in alleviating the severity of alarm overloading. The essence of designing effective multivariate alarm systems is to exploit some features of involved multiple process variables to detect abnormal conditions in an online manner. Brooks et al. (2004) compared current data points to a pre-defined normal operating zone of multiple process variables via a geometric process control technique, and visualized dynamic alarm thresholds in parallel coordinates. Charbonnier et al. (2005) and Charbonnier and Portet (2012) extracted the trends of process variables for complex process monitoring and fault diagnosis. Izadi et al. (2009) and

夽 This work was supported by the National Natural Science Foundation of China under Grant No. 61433001 and the Research Fund for the Taishan Scholar Project of Shandong Province of China. ∗ Corresponding author. E-mail addresses: [email protected] (K. Chen), [email protected] (J. Wang). http://dx.doi.org/10.1016/j.cherd.2017.04.011 0263-8762/© 2017 Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved.

12

chemical engineering research and design 1 2 2 ( 2 0 1 7 ) 11–21

(a)

Table 1 – Cross-industry study (Rothenberg, 2009).

2

EEMUA

Oil-gas

PetroChem

Power

144 10

1200 220

1500 180

2000 350

X

0

Average alarms/day Peak alarms/10 min

−2 −4

1. Some online methods modify the extracted trends in the past after the most recent trends are updated, so that they are not implemented in a real online sense (see e.g., Keogh et al., 2001; Charbonnier et al., 2005; Charbonnier and Gentil, 2007; Maurya et al., 2010; Charbonnier and Portet, 2012). Taking the sliding window method in Keogh et al. (2001) as an example, this method fits a signal using a linear segment growing with the time until the fitting error becomes larger than a threshold, and then updates the linear segment by restarting the fitting procedure. As

0

50

100

150

200

(b) 0.01 t=1~100 t=1~200

0 dX

Kondaveeti et al. (2009) applied multivariate statistics such as principal component analysis to generate alarms more efficiently. Gupta et al. (2013) integrated multiple techniques such as wavelet analysis and principal component analysis to rationalize alarm thresholds of process variables to alleviate noise effects and detect faults earlier. Zhu et al. (2014) proposed dynamic alarm thresholds being adaptive for different pre-defined operating stages varying with time. Alrowaie et al. (2014) applied particle filters for multivariate nonlinear stochastic systems to design alarm thresholds. Zang and Li (2014) and Han et al. (2016) optimized alarm thresholds by minimizing false and missed alarm probabilities on the basis of joint probability densities of process variables in the normal and abnormal conditions. Variational directions of related process variables are exploited here as the features to detect abnormal conditions. This idea is based on an observation on how industrial plant operators find abnormal conditions. Let us take a variable frequency pump as an example. It can be easily observed that the outlet flow rate is always increasing or decreasing in a synchronized manner with the pump speed, when the pump works in normal conditions. This observation is guaranteed by physical laws of the pump operation. Hence, if the outlet flow and speed of a variable frequency pump are found moving in opposite directions, then there must be something wrong in the pump operation. Thus, by looking at the directions of variations in the outlet flow rate and the pump speed, industrial plant operators are able to judge whether the pump is in the abnormal condition or not. Variational directions can be obtained by qualitative trend analysis (QTA) methods that extract the increasing, decreasing or steady trends of signals. Several QTA methods have been proposed in literature. Signals are approximated by several piecewise linear segments (Keogh et al., 2001; Charbonnier et al., 2005; Charbonnier and Gentil, 2007; Charbonnier and Portet, 2012) or high-order polynomials (Savitzky and Golay, 1964; Maurya et al., 2010). Splines methods and wavelet methods decompose signals into different scale bases (Vedam et al., 1998; Flehmig et al., 1998; Flehmig and Marquardt, 2006). Clustering methods use a fuzzy C-means clustering technique to identify different trends (Melek et al., 2005). Primitives methods represent signals by a combination of primitives in series whose local trends are pre-defined (Janusz and Venkatasubramanian, 1991; Rengaswamy and Venkatasubramanian, 1995; Colomer et al., 2002). However, the above-mentioned methods may have two limitations:

−0.01 −0.02

0

50

100 t

150

200

Fig. 1 – (a) The time sequence plot of the signal x(t), (b) the time sequence plot of the difference signal xd (t) for t ∈ [1, 100] (upper) and that for t ∈ [1, 200] (lower). a numerical illustration, Fig. 1 shows the time sequence plots of a signal x(t) (upper subplot) and its backward difference xd (t) = xˆ (t) − xˆ (t − 1) (lower subplot), where xˆ (t) is an approximated signal of x(t) obtained by the sliding window method. The extracted trends of all samples before the time instant t = 100 are larger than zero; however, the trends are modified to take negative values at the time instant t = 200. The extracted trends before t = 100 and after t = 200 are erroneously inconsistent. 2. Few methods are able to adapt extracted trends to the coexisted sharp and smooth variations in signals. For sharp variations, narrow time windows or equivalently smaller time scales should be used in order to detect the changes of trends promptly; by contrast, for smooth variations, larger time windows or larger time scales should be exploited in order to reduce negative effects of noises on the trend extraction. As an illustration, Example 1 in Section 4 later illustrates this limitation for a widely used method, the Savitzky–Golay (S–G) filter (Savitzky and Golay, 1964). This paper proposes a new method to design a multivariate alarm system based on variational directions of involved process variables. Due to the above two limitations, the existing QTA methods are not adopted here to calculate the variational directions. Instead, a so-called adaptive time gradient (ATG) approach is formulated to calculate the time gradients of signals and the variational directions by using adaptive time scales. The key step of the ATG approach is to determine adaptive time scales from a function between time scales and volatilities of signals. The function parameters are estimated from training data sets. Compared to the existing QTA methods, the ATG approach is adaptive in time scales and is online in a real sense. Alarms are ready to arise, when the variational directions have a conflict with a nominal relationship of variational directions among process variables in normal conditions obtained from some process knowledge. To the best of our knowledge, this paper perhaps is the first one in literature to exploit the variational directions of related process variables in designing alarm systems. The rest of the paper is organized as follows. Sections 2 and 3 present the main idea and the details of the proposed method. Sections 4 and 5 respectively provide numerical and industrial examples to illustrate the effectiveness of the proposed method. Section 6 concludes the paper.

13

chemical engineering research and design 1 2 2 ( 2 0 1 7 ) 11–21

2.

The main idea

of x (t) is defined as the time gradient k(n) associated with a proper time scale C, i.e.,

This section introduces the main idea of the proposed method to design multivariate alarm systems based on variational directions.

Tx (n) = kC (n) .

2.2. 2.1.

Multivariate alarm systems

Time gradients

Consider a discrete-time univariate signal x (t) where t ∈ Z is the sampling time index, and Z is the set of non-negative integers. To calculate the time gradient in an online manner, the current time index is denoted as n. The time gradient at t = n is defined as the slope k (n) of a local linear regression model on the time interval t ∈ [n − L + 1, n], where L is the time window size, i.e., x (t) = k (n) t + b (n) ,

t ∈ [n − L + 1, n] .

n t=n−L+1

2 [x (t) − k (n) t − b (n)] K (n, t) ,

(2)

where K (n, t) is the kernel function, K (n, t) = n−t .

(3)

Here  ∈ (0, 1) is a constant forgetting factor so that n−t is no larger than 1 and decreases with t. Define a time scale C as the time duration that the weighting term n−t decreases from 1 to a small positive real number , i.e., C = .

(4)

A default value of  is 1%. The physical meaning of C is more transparent than the counterpart of , even though the two are closely related. A suitable range of C will be analyzed later in Section 3.1.1. Bakshi and Stephanopoulos (1994) pointed out that the resolution of a trend is closely related to the concept of time scales, i.e., if a trend takes a higher resolution, then a greater local change rate of the trend can be observed and vice versa. As a consequence, a proper range of time scales is required in analyzing qualitative trends of signals. This is achieved as follows. Setting the partial derivatives of G with respect to k (n) and b (n) to zeros and solving the corresponding equations yield the optimal solutions to minimize G as

 kˆ (n) =

tx (t) n−i t



n−t − t





x (t) n−i t



 

tn−t t

,

2 2 t n−t t n−t − tn−t t t 2 x (t) n−t t t n−t − tx (t) n−t t tn−t t t 2 2 t n−t t n−t − tn−t t t

 bˆ (n) =













Consider a multivariate alarm system involving multiple process variables, X:= [X1 , X2 , . . ., Xm ]. The data sequence of each variable Xi with i = 1, 2, . . ., m is represented by xi (t). The variational direction of Xi is defined as





sign TXi :=

⎧ +1 TXi ∈ (, +∞) ⎪ ⎨ ⎪ ⎩

0

TXi ∈ [−, ]

,

(7)

−1 TXi ∈ (−∞, −)

(1)

The loss function for the local regression model is G=

(6)





(5)

where TXi is the time gradient of Xi . Here  is the significance threshold to  belater determined in Section 3.1.2. In words, when sign TXi equals to +1, 0 and -1, Xi is increasing, unchanged (or having no evident trends) and decreasing in a statistical sense, respectively. If X is in normal conditions, then the variational directions of X1 , X2 , . . ., Xm are assumed to have a fixed relationship. Mathematically, the nominal relationship is represented by f (sign (TX1 ) , sign (TX2 ) , . . ., sign (TXm )) = 0.

(8)

By contrast, when X runs into abnormal conditions, the relationship in (8) is assumed to be no longer valid, i.e., f (sign (TX1 ) , sign (TX2 ) , . . ., sign (TXm )) = / 0. The nominal relationship in (8) is usually obtained based on some process knowledge of X. Under the above assumptions, our objective is to calculate time gradients and variational directions of X and to raise an alarm for X being in abnormal conditions in an online manner.

3.

The proposed method

The first part of this section proposes a so-called adaptive time gradient (ATG) approach to extract time gradients of a univariate signal, and the second part formulates a novel method to exploit the ATG approach to design a multivariate alarm system.

3.1.

Adaptive time gradient (ATG) approach

In this subsection, a function is first established to connect time scales with volatilities; second, parameters of the established function are determined from training data sets; finally, the steps of the proposed ATG approach are given.

,

3.1.1.

The function between time scales and volatilities

There are two basic assumptions as follows: where the index t ranges from n − L + 1 to n. The optimal estimate of x (n) is xˆ (n) = kˆ (n) t + bˆ (n). It is apparent that the time gradient k (n) is highly dependent on the forgetting factor . A higher value of , or equivalently a larger time scale C, implies that k (n) takes account of the information from more historical data samples. Hence, at the current time index n, the time gradient Tx

1. Signals are composed of different variations in various time scales. The smaller time scale corresponds to sharp variations, while the larger time scale leads to variation lasting for a longer period of time. 2. In order to adapt to time-varying variations in signals, adaptive time scales have to be used, and the current

14

chemical engineering research and design 1 2 2 ( 2 0 1 7 ) 11–21

time scale C (n) can be determined from the corresponding volatility V (n), i.e., C (n) = g (V (n)) .

(9)

Let us calculate the volatility of the training data set in order to fit the function (12). Define the time scale 1 Cp = Cmin + p − 1 for p = 1, 2, . . ., Cmax − Cmin + 1, and p =  Cp . We calculate the time gradient kCp (j) for each data sample s (j) with j = L, L + 1, . . ., l as that in (5)

Owing to the mapping relationship between  and C in (4), an equivalence of (9) is



 (n) = f (V (n)) .

max

k (t) −

t ∈ [n−C(n)+1,n]

min

k (t) .

t ∈ [n−C(n)+1,n]

(11)

A proper range [Cmin , Cmax ] of the time scale C for a process variable could be determined based on empirical knowledge as follows. The minimum time scale Cmin is oriented from a common sense that if an abnormality arises, it should last at least for a while before being regarded as an actual abnormality, instead of being treated as short-time random variations due to noises. Thus, Cmin is a boundary parameter in time to distinguish noises and trends in the finest resolution. The rationale of Cmin is similar to alarm delay timers (ISA, 2009; EEMUA, 2013), whose role is to raise an alarm when several consecutive samples of a process variable overpass an alarm threshold. A default value of Cmin here is taken as Cmin = 20, by taking the recommended value of 20 seconds for alarm delay timers (Wang and Chen, 2014). The maximum time scale Cmax is associated with trends in the roughest resolution. It is oriented from the influence of some events, at the time duration of a specific trend. Usually Cmax can be obtained from historical data sequences as the maximum duration of observed trends. When V (n) is small,  (n) approaches to 1/Cmax and C (n) is forced to be no larger than Cmax . In the opposite side, a higher value of V (n) corresponds to  (n) approaching 1/Cmin and C (n) is constrained to be no less than Cmin . It is a reasonable fact that a shorter (longer) time scale would be more suitable for a higher (lower) volatility case. Hence, (10) is formulated as  (n) = max (−ˇ0 · max (0, V (n) − V0 ) + 1/Cmax , 1/Cmin ) ,

(12)

where ˇ0 and V0 are two constants to be determined in the next subsection.

3.1.2. sets

ts (t) p

Determination of parameters from training data

This subsection determines the parameters ˇ0 and V0 in (12) and  in (7) from some training data sets.



l

Select a data segment s (t) := [s (1) , s (2) , . . ., s (l)] as the t=1 training data set with l ≥ L. The function (12) describes a quantitative relationship between time scales and volatilities. The user-selected range [Cmin , Cmax ] of time scales has been defined for the variations in signals of interest. Thus, fitting the function (12) is mainly for capturing the volatilities induced by noises. Hence, a requirement of s(t) is that the data segment has no notable variations on time scales less than Cmin in the view of industrial plant operators. A mild assumption should be mentioned that the variations on time scales less than Cmin are assumed to be absent, and the variations induced by noises in the training data set are expected to be persistent in the subsequent online data sequences.





j−t

t2 p t

(10)

The volatility is defined as the range of time gradients in the current time scale, i.e., V (n) =

j−t

t

kCp (j) =

j−t  t p



−

 t

j−t

 − t p

j−t

s (t) p





j−t

t

2

j−t

tp

(13)

,

t t p

where the index t ranges from j − L + 1 to j. Since the maximum meaningful time scale is Cmax , the window size L of the local regression model in (1) should be no larger than Cmax . Because the forgetting factor p is changing adaptively, L can be simply chosen as L = Cmax . Next, the volatility VCp (j) for s (j) on the time scale Cp is obtained as VCp (j) =

max

t ∈ [j−Cp +1,j]

kCp (t) −

min

(14)

kCp (t) .

t ∈ [j−Cp +1,j]

Repeat the calculation of kCp (j) and VCp (j) for all values of p and j so that both kCp (j) and VCp (j) are two dimensional matrices with j as the row index and p as the column index. Plot each row of VCp (j) versus Cp or p as the dot lines in Fig. 2, which is obtained from an industrial data set of a pump speed in Example 3 appeared later in Section 4. The envelope curve of all the dot lines marks the maximum volatility (Env) (Env) (Env) corresponding to a specific time scale Cp or p . Any Vp (Env)

data samples with volatility less than Vp time scale

(Env) Cp

on the specific

should be associated with a smoother trend (Env)

so that a longer time scale larger than Cp should be chosen. Then only the data samples with the highest volatility on each time scale are expected to be on the function curve   of (12). Hence, the envelope curve denoted as



(Env)

Cp

(Env)



, Vp

(Env)

p

(Env)

, Vp

or

is calculated from

⎧ 1 ⎪ ⎪ (Env) ⎪ Cmin + p − 1 ⎪  := p ⎪ ⎪ ⎪ ⎨ (Env)

Cp :=Cmin + p − 1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Vp(Env) = max V j ∈ [Cmax ,l]

(15)

,

(Env)

Cp

(j)

where p = 1, 2, . . ., Cmax − Cmin + 1. The least-squares method fit the function (12) to the envelope curve is used to   (Env)

p

(Env)

, Vp

with the optimal solution as:

⎧  (Env) V = VCmax −C +1 ⎪ ⎪ min ⎨ 0 Cmax −Cmin +1

2 .  (Env) (Env)  ⎪ f Vp , ˇ0 , V0 − p ⎪ ⎩ ˇ0 = argˇmin 0

(16)

p=1

A goodness of fit is defined as

Cmax −Cmin +1  =1−

(Env)

− ˆ p

(Env)

(Env) − ¯ p

p=1

p

p=1

p

Cmax −Cmin +1

(Env)

2

2 ,

(17)

15

chemical engineering research and design 1 2 2 ( 2 0 1 7 ) 11–21

(b)

(a) 1

500 Train Data Envelope Fitting Curve

Train Data Envelope Fitting Curve

450

0.95

400

350

300

λ

C

0.9

0.85

250

200

150

0.8

100

50

0.75

0

5

10

0

15

0

5

Volatility

10

15

Volatility

Fig. 2 – (a) The volatility VCp (j) versus p , and (b) the volatility VCp (j) versus Cp . (Env) ˆ p =f

where 1 Cmax −Cmin +1



(Env)

Vp

Cmax −Cmin +1

, ˇ0 , V0 ,

(Env) ¯ p =

and

(Env)

p . Generally  > 0.8 indicates i=1 that the quality of the fitting function is satisfactory. The parameter  in (12) is a threshold to tolerate small values of time gradients induced by noises. The variational direction is unchanged if the time gradient is inside the range [−, ]. The unchanged variational direction means that the volatility is low. The minimum volatility corresponds to the maximum time scale according to (12). Hence, the significance threshold  in (7) is determined as  = 2std (kCmax (j)) ,

C (n) = min

where kCmax (j) is from (13) and j ∈ [Cmax , l]. Eq. (18) is developed from a hypothesis test on mean values with the null hypothesis of zero time gradients by using a significance level 0.05.

Steps of the ATG approach

The proposed ATG approach consists of two parts, namely, the offline training of the unknown parameters {V0 , ˇ0 , } and the online calculating of time gradient Tx (n). Offline training. Determine a time scale range [Cmin , Cmax ] and select a training data set s (t) as described in Section 3.1.2. Obtain the unknown parameters {V0 , ˇ0 , } from (16) and (18). If  in (17) is less than 0.8, choose another training data set and repeat this part of offline training. Online calculating. The following steps are implemented when a new data sample x (n) is available. Step 1: Calculate the time gradient Tx (n) as that in (5) by using the time scale C (n − 1) or equivalently the forgetting factor  (n − 1), i.e., Tx (n) = kC(n−1) (n)



t

tx (t) (n − 1)



t

2

n−t

log () log ( (n − 1))

Cmax ,

 .

Step 3: Calculate the volatility for the updated time scale C (n), V (n) =

k (t) −

max

t ∈ [n−C(n)+1,n]

min

k (t) .

t ∈ [n−C(n)+1,n]

Step 4: Update the forgetting factor,  (n) = max (−ˇ0 · max (0, V (n) − V0 ) + 1/Cmax , 1/Cmin ) .

3.2.

3.1.3.



(18)

j

=

Step 2: Update the time scale as



t (n − 1)

t

(n − 1)

n−t



t

n−t

−



(n − 1)

t n−t

x (t) (n − 1)

−



t

n−t



t(n − 1)

where the index t ranges from n − Cmax + 1 to n.

t(n − 1)

2

t n−t

n−t

Design of multivariate alarm systems

For the multivariate systems X = [X1 , X2 , . . ., Xm ], the time gradients of Xi for i = 1, 2, . . ., m are calculated by the ATG approach in Section 3.1.3 in an online manner. Then, the variational direction of Xi is obtained from (7). Next, the multivariate alarm system is designed as follows. To compare the variational directions with the relationship in (8) in a systematic manner, a matrix R with size r × m is introduced as

⎡

sign (T1,X1 )

sign (T1,X2 )

···

sign (T1,Xm )

⎤

⎢ ⎥ ⎢ sign (T2,X1 ) sign (T2,X2 ) · · · sign (T2,Xm ) ⎥ ⎢ ⎥ R:= ⎢ ⎥. .. .. .. ⎢ ⎥ .. ⎣ ⎦ . . . . sign (Tr,X1 )

sign (Tr,X2 )

···

(19)

sign (Tr,Xm )

Each row of R represents one rule against the relationship in (8) under a specific circumstance. For example, if X := [X1 , X2 ] is in the normal condition when both variables are increasing, decreasing or unchanged simultaneously, then the relationship matrix R is

,

 R=

1

−1

−1

1

 .

(20)

16

chemical engineering research and design 1 2 2 ( 2 0 1 7 ) 11–21





sign (TX (n)) = sign (TX1 (n)) , sign (TX2 (n)) , . . ., sign (TXm (n)) .

(a) 1000 Signal

The variational directions of all process variables in the current time index n are represented by a vector

0 −1000

(21)

0

200

400

600

800

1000

600

800

1000

600

800

1000

(b)

4.

0.5 0 −0.5 −1 −1.5

k

The variational direction vector sign (TX (n)) is compared with the relationship matrix R to determine whether X(n) is in the normal or abnormal condition. In particular, if the variational direction vector is the same as one row of R, then X (n) is in the abnormal condition, and the alarm signal Xa (n) = 1; otherwise, X (n) is regarded as being in the normal condition, and no alarm arises, i.e., Xa (n) = 0.

0

200

400 (c)

C

100 80 60 40 20

Examples

0

200

400 t

Three examples are provided in this section. The first numerical example illustrates the effectiveness of the ATG approach. The second one designs a multivariate alarm system for a water tank system. The last example is an industrial case study for a variable-frequency pump.

Fig. 4 – The time sequence plots of signals: (a) x (t) (dot), (b) the time gradient kC (t), and (c) the time scale C (t) in Example 1.

0.4

Example 1. Consider a signal e(t),

for t ∈ [1, 300)

−1 ∗ (t − 300) + e(t),

for t ∈ [300, 1000]

0

where e(t) is the Gaussian white noise with zero mean and unit variance. Fig. 3 shows the result of the offline training part with Cmin = 20, Cmax = 100, and s (t) = x (1 : 200). The unknown parameters ˇ0 and V0 in (12) and  in (7) and the goodness of fit  in (17) are determined,

−0.2

−0.4 DX



x(t) =

Expected ATG S−G

0.2

−0.6

−0.8

V0 = 0.0481, ˇ0 = 0.3242,  = 0.0073,  = 0.9843.

−1

In particular,  = 0.9843 says that the quality of the fitting function is satisfactory. The time sequence plots of the signal x(t), the time gradient kC (t), and the time scale C(t) are shown in Fig. 4. The notable change occurring on t = 300 is well captured by the ATG approach, so that the time scale C (t) decreases to near Cmin . As x (t) keeps a steady descent trend after t = 300, C (t) grows to near Cmax and the time gradient kC (t) becomes smoother. The results in Fig. 4 are benchmarked by the noise-free counterpart in Fig. 5. As a comparison, a widely used method,

−1.2

−1.4 0

100

200

300

400

500

600

700

800

900

1000

t

Fig. 5 – The time sequence plots of the expected time gradient (circle), the estimated time gradient from the ATG approach (dot) and the counterpart from the S–G filter (triangle) in Example 1.

(a)

0.2

0.96 Train Data Envelope Fitting Curve

0.94

0

Estimation error

0.92

0.9

λ

0.88

−0.2

−0.4

−0.6 0.86 −0.8

0.84

0.82

−1 ATG

0.8

0.78 0

0.1

0.2

0.3

0.4

0.5

0.6

Volatility

Fig. 3 – The volatility VC (j) versus  in Example 1.

0.7

S−G

Fig. 6 – The boxplots of the estimation errors from the ATG approach and the counterpart from the S–G filter in Example 1.

17

chemical engineering research and design 1 2 2 ( 2 0 1 7 ) 11–21

(b) Signal

Signal

(a) 40 20 0

0

500 (c)

k

k

0

500 (e)

300 200 100 0 1 0.5 0

500 (d)

1000

0

500 (f)

1000

0

500 t

1000

100 C

C

In Out H

40 20 0

S2

40 20 0

0

1.5 1 0.5 0

1000

100 50

Fig. 7 – Schematic diagram of a water tank system in a laboratory.

200

1000

0 −0.5 −1 0

400

0

500 t

50

1000

Fig. 10 – The time sequence plots of signals: (a) Fo (t), (b) H (t), (c) the time gradient kFo (t), (d) the time gradient kH (t), (e) the time scale CFo (t) and (f) the time scale CH (t) in Example 2.

0

200

400

600

800

1000

the change of x (t) at t = 300 immediately, the ATG approach gives a much smoother and more accurate time gradient than the S–G filter. This observation is also supported by the boxplots of the estimation errors in Fig. 6. The superior performance of the ATG approach is owing to the adaptivity of time scales. In other words, the ATG approach is able to smooth (track) signals on a larger (smaller) time scale than the S–G filter when a trend with low (high) volatility is detected.

t

Fig. 8 – The time sequence plots of signals: (a) Fi , (b) Fo , (c) H and (d) S2 in Example 2. the Savitzky–Golay (S–G) filter (Savitzky and Golay, 1964), is applied. Since the S–G filter is associated with a fixed time window for computation of the time gradient, the lower bound Cmin is used as the time window size. It can be found that even though both the S–G filter and the ATG approach response to

Example 2. Consider a water tank system schematically shown in Fig. 7 in a laboratory at Peking University. The system can be described by a dynamic model

⎧ ⎨ S1 dH = Fi − Fo dt  ⎩ Fo = S2 ·

(a)

(22)

.

2gH + e (b)

0.96

0.96 Train Data Envelope Fitting Curve

0.94

Train Data Envelope Fitting Curve

0.94

0.92

0.9

0.9

0.88

0.88 λ

λ

0.92

0.86

0.86

0.84

0.84

0.82

0.82

0.8

0.8

0.78

0

0.1

0.2 Volatility

0.3

0.4

0.78

0

0.05

0.1 0.15 Volatility

0.2

0.25

Fig. 9 – (a) The volatility VC (j) versus  for Fo , and (b) the volatility VC (j) versus  for H in Example 2.

18

chemical engineering research and design 1 2 2 ( 2 0 1 7 ) 11–21

respectively. The two groups of unknown parameters ˇ0 and V0 in (12) and  in (7) and the goodness of fit  in (17) for Fo and H are respectively

(a)

Label

1

⎧ V0,Fo = 0.0444 ⎪ ⎪ ⎪ ⎪ ⎨ ˇ0,Fo = 0.7425

0.5 0 0

200

400

600

800

1000

(b)

⎪ Fo = 0.0134 ⎪ ⎪ ⎪ ⎩ Fo = 0.9919

Alarm

1 0.5 0 0

200

400

600

800

1000

Fig. 11 – The actual alarm signal (a) and the alarm signal Xa (t) from the designed multivariate alarm system (b) in Example 2. Here Fi and Fo respectively denote the inlet and outlet flows passing through two valves, H is the liquid level, S1 and S2 are the sectional areas of the tank and the outlet valve, respectively, g represents the gravity constant and e is the system noise. By adjusting the inlet flow Fi via a control valve, operators could change the tank level H to a desired position, while the outlet flow Fo varies with H. Thus, (22) implies that Fo and H increase or decrease simultaneously when the system is in normal conditions. Fig. 8 shows the time sequence plots of Fi , Fo and S2 in simulations. The inlet flow Fi is increased on t = 600, so that both Fo and H are increasing. A fault occurs such that S2 rapidly decays by manually blocking the outlet valve from t = 821 to t = 860. As a result, an abnormality occurs, because Fo decreases while H increases. Hence Fo and H are in the abnormal condition in the time period t ∈ [821, 860]. The proposed method is applied with Cmin = 20, Cmax = 100, s (t) = X (1 : 400, :), and the relationship matrix R in (20). Fig. 9 shows the results of the offline data training for Fo and H,

(a)

(b) 170

550

168

500

166

Amplitude

450 400 350 300 250

164

⎪ H = 0.0161 ⎪ ⎪ ⎪ ⎩

.

H = 0.9678

The fitness metrics Fo and H are close to 1, saying that the qualities of corresponding functions in (12) are satisfactory. The proposed method yields the time gradients kFo (t) and kH (t) and time scales CFo (t) and CH (t) as shown in Fig. 10. The red horizontal lines in the middle and bottom suplots in Fig. 10 denote the significance thresholds [−, ] in (7) and the range [Cmin , Cmax ] of time scales, respectively. As expected, the time scales CFo (t) and CH (t) take the maximum values when Fo (t) ad H(t) are steady, and decrease quickly when Fo (t) ad H(t) experience sharp variations. As presented in Fig. 11, the alarm system raises alarms in the time period t = [832, 875], which is close to the actual time interval t ∈ [821, 860] that the abnormality presents. Thus, with an expected detection delay between the detected and actual abnormalities less than Cmin = 20, the designed multivariate alarm system successfully detect the abnormality against the nominal relationship between Fo and H. Example 3. A variable-frequency feedwater pump, from a large-scale thermal power plant located at Shandong Province in China, is considered here. The current alarm system to monitor the operation of the feedwater pump raises an alarm when the inlet flow rate is larger than a high alarm threshold 560 t/h. Such a univariate alarm system sometimes results in false alarms, when the inlet flow rate of the feedwater pump increases beyond 560 t/h to supply more water to downstream devices. As shown in Fig. 12, even though the inlet flow rate overpasses the high alarm threshold 560 t/h, the relationship between the inlet flow rate and pre-pump electrical current are unchanged. Thus, the feedwater pump is in normal conditions, and the occurred alarms from the univariate alarm system are nuisance ones. Furthermore, the univariate alarm system may also lead to missed alarms, to be illustrated by the data sequences in Fig. 13 later. To avoid false and missed alarms from the current univariate alarm system, we formulate a three-variable alarm system X := [X1 , X2 , X3 ] for the pump, where X1 , X2 and X3 are the

162

(a) 600

160 X1

Prepump current (A)

600

,

⎧ V0,H = 0.0481 ⎪ ⎪ ⎪ ⎪ ⎨ ˇ0,H = 1.0961

400

158 200

156 (b)

200

6000

154 0

1 Time (sec)

152 450

2 4

x 10

X

2

150 500 550 Pump flow (t/h)

600

5000 4000 (c) 200

X3

Fig. 12 – (a) The time sequence plots of the inlet flow rate (the upper one) and pre-pump electrical current (the lower one), and (b) the scatter plot of the inlet flow rate and pre-pump electrical current with the high alarm threshold (red dash-dotted) in Example 3. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

150 100

0

500

1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 t

Fig. 13 – The time sequence plots of three process variables: (a) X1 , (b) X2 and (c) X3 in Example 3.

19

chemical engineering research and design 1 2 2 ( 2 0 1 7 ) 11–21

(a)

(b)

1

(c)

1

1

Train Data Envelope Fitting Curve

Train Data Envelope Fitting Curve

Train Data Envelope Fitting Curve

0.95

0.95

0.9

0.9

0.9

λ

λ

λ

0.95

0.85

0.85

0.85

0.8

0.8

0.8

0.75

0

2 4 Volatility

0.75

6

0

5 10 Volatility

15

0.75

0

0.5 Volatility

1

Fig. 14 – (a) The volatility VC (j) versus  for X1 , (b) the volatility VC (j) versus  for X2 , and (c) the volatility VC (j) versus  for X3 in Example 3. inlet flow rate, the pump speed and the pre-pump electrical current, respectively. Based on physical laws of pumps, the three process variables are expected to increase or decrease synchronously in normal conditions, so that the relationship matrix R is

 R=

1

1

−1 −1 −1

1

1

1

1

−1

−1

0

0

0

−1

0

1

−1

−1

−1

1

1

−1

−1

−1

−1

1

1

−1

0

1

1

−1

1

1

0

T .

Fig. 13 presents the time sequence plots of the three process variables, collected on January 24, 2014 with the sampling

(b)

(c)

5500

180

500

5000

160

400 300

0

2000

4000

6000

Signal

600 Signal

Signal

(a)

period 1 sec. By visualization, X is in the normal condition before t = 5035. After that, there are two intervals with obvious abnormalities: X2 increases while X1 and X3 keep decreasing during the interval t ∈ [5035, 5160], and X2 decreases while X1 and X3 keep increasing during the interval t ∈ [5195, 5239]. This is an offline detection of normal and abnormal conditions by visualization. Let us do so in an online manner via the proposed method. Note that the univariate alarm system currently used in the power plant is not able to detect such abnormalities, and thus results in missed alarms.

4500 4000

0

2000

(d)

4000

120

6000

0

−20

6000

0

2000

4000

−0.2

6000

400

C

200

4000

6000

4000

6000

C

400 C

400

t

2000 (i)

600

2000

0

(h) 600

0

6000

0

(g)

0

4000

0.2

600

200

6000

0.4

−10

4000

4000 (f)

k 0

2000

2000

0.6

k

2 k

10

0

0

(e)

4

−2

140

0

200

0

2000

4000 t

6000

0

0

2000 t

Fig. 15 – The time sequence plots of signals: (a) x1 (t), (b) x2 (t), (c) x3 (t), (d) the time gradient kX1 (t), (e) the time gradient kX2 (t), (f) the time gradient kX3 (t), (g) the time scale CX1 (t), (h) the time scale CX2 (t) and (i) the time scale CX3 (t) in Example 3.

20

chemical engineering research and design 1 2 2 ( 2 0 1 7 ) 11–21

1

1 0.8

0.8

0.4

0.4

0.2

Alarm

Alarm

0.6

0.6

0.2

0

−0.2

0 0

500

−0.4

1000 1500 2000 2500 3000 3500 4000 4500 5000 5500

Fig. 16 – The alarm signal Xa (t) for signals in Fig. 15 in Example 3.

−0.6 −0.8

The proposed method is applied by setting Cmin = 20, Cmax = 500 and s (t) = X (1 : 2000, :). Fig. 14 shows the offline training results. The unknown parameters and the goodness of fit for X1 , X2 and X3 are respectively

⎧ V0,X1 = 0.2096 ⎪ ⎪ ⎪ ⎪ ⎨ ˇ0,X = 0.0629 1

⎪ X1 = 0.1093 ⎪ ⎪ ⎪ ⎩ X1 = 0.9886

,

⎧ V0,X2 = 0.8608 ⎪ ⎪ ⎪ ⎪ ⎨ ˇ0,X = 0.0178 2

,

⎪ X2 = 0.4496 ⎪ ⎪ ⎪ ⎩ X2 = 0.9962

⎧ V0,X3 = 0.0319 ⎪ ⎪ ⎪ ⎪ ⎨ ˇ0,X = 0.3711 3

⎪ X3 = 0.0173 ⎪ ⎪ ⎪ ⎩

(23)

.

X3 = 0.9877

The fitness metrics X1 , X2 and X3 are close to 1, saying that the qualities of corresponding functions in (12) are satisfactory. Fig. 15 presents the time gradients and time scales for X1 , X2 and X3 . The red horizontal lines in the middle and bottom suplots in Fig. 15 denote the significance thresholds [−, ] in (7) and the range [Cmin , Cmax ] of time scales, respectively. Before t = 5000, the volatilities of the pump system are low so that large time scales are used to calculate the time gradients. The time scales decay rapidly when some sharp variations

(a)

500

4000

5000

165

4800

160

0

2000

4000

150

6000

0.02

0.5

0.01

k

1

0.1

k

0.2

0 2000

4000

−0.5

6000

0

2000

(g)

4000

6000

−0.01

200

4000 t

2000

6000

4000

6000

4000

6000

C

400

C

400

C

400

2000

6000

(i) 600

0

0

(h) 600

200

4000

0

600

0

2000 (f)

0.03

0

0

(e) 1.5

−0.1

6000

155

(d)

0

5000

(c) 170

4400

6000

4000

are detected around t = 5080. The designed multivariate alarm system yields the alarm signal Xa (t) in Fig. 16. The alarms arise at the time instants t = [5051, 5178] and t = [5215, 5239], which are close to the offline detection of abnormal conditions. The designed multivariate alarm system is able to avoid the false alarms illustrated in Fig. 12. By using the offline training results in (23), the proposed method yields the time sequence plots of the original signals, their time gradients and adaptive time scales as shown in Fig. 17. The increment of X1 over the original alarm threshold 560 t/h is accompanied by synchronized increments in X2 and X3 around t = 2750. Thus, no alarms occur for such a normal condition in Fig. 18. Note that the time scales in Fig. 17 are able to adapt to the sharp variations in signals around t = 2750, and to the smooth variations in other time periods.

4600 2000

3000

Fig. 18 – The alarm signal Xa (t) for signals in Fig. 17 in Example 3.

Signal

Signal

Signal

550

0

2000

5200

0.3

k

1000

(b)

600

450

−1 0

0

200

0

2000

4000 t

6000

0

0

2000 t

Fig. 17 – The time sequence plots of signals: (a) x1 (t), (b) x2 (t), (c) x3 (t), (d) the time gradient kX1 (t), (e) the time gradient kX2 (t), (f) the time gradient kX3 (t), (g) the time scale CX1 (t), (h) the time scale CX2 (t) and (i) the time scale CX3 (t) in Example 3.

chemical engineering research and design 1 2 2 ( 2 0 1 7 ) 11–21

5.

Conclusion

This paper proposed a method to design a multivariate alarm system based on variational directions of involved process variables. A so-called adaptive time gradient (ATG) approach was formulated to calculate the time gradients of signals using adaptive time scales, in order to calculate the variational directions afterwards. The adaptive time scales were determined from a function between time scales and volatilities of process variables, and the function parameters were estimated from training data sets. The designed multivariate alarm system raised alarms when the calculated variational directions were inconsistent with a prior-known relationship among variational directions of process variables in normal conditions. Numerical and industrial examples illustrated the effectiveness of the proposed method.

References Alrowaie, F., Gopaluni, R., Kwok, K., 2014. Alarm design for nonlinear stochastic systems. In: 11th World Congress on Intelligent Control and Automation (WCICA2014), Shenyang, China, pp. 473–479. Bakshi, B.R., Stephanopoulos, G., 1994. Representation of process trends. III. Multiscale extraction of trends from process data. Comput. Chem. Eng. 18, 267–302. Bransby, M., Jenkinson, J., 1998. The Management of Alarm Systems. Health and Safety Executive. Brooks, R., Thorpe, R., Wilson, J., 2004. A new method for defining and managing process alarms and for correcting process operation when an alarm occurs. J. Hazard. Mater. 115, 169–174. Charbonnier, S., Garcia-Beltan, C., Cadet, C., Gentil, S., 2005. Trends extraction and analysis for complex system monitoring and decision support. Eng. Appl. Artif. Intell. 18, 21–36. Charbonnier, S., Gentil, S., 2007. A trend-based alarm system to improve patient monitoring in intensive care units. Control Eng. Pract. 15, 1039–1050. Charbonnier, S., Portet, F., 2012. A self-tuning adaptive trend extraction method for process monitoring and diagnosis. J. Process Control 22, 1127–1138. Colomer, J., Melendez, J., Gamero, F.I., 2002. Qualitative representation of process trends for situation assessment based on cases. In: 15th IFAC World Congress, Barcelona, Spain, pp. 103–108. EEMUA, 2013. EEMUA-191: Alarm Systems – A Guide to Design, Management and Procurement. Engineering Equipment and Materials Users Association. Flehmig, F., Watzdorf, R., Marquardt, W., 1998. Identification of trends in process measurements using the wavelet transform. Comput. Chem. Eng. 22, S491–S496. Flehmig, F.u., Marquardt, W., 2006. Detection of multivariable trends in measured process quantities. J. Process Control 16, 947–957. Gupta, A., Giridhar, A., Venkatasubramanian, V., Reklaitis, G.V., 2013. Intelligent alarm management applied to continuous

21

pharmaceutical tablet manufacturing: an integrated approach. Ind. Eng. Chem. Res. 52, 12357–12368. Han, L., Gao, H., Xu, Y., Zhu, Q., 2016. Combining FAP, MAP and correlation analysis for multivariate alarm thresholds optimization in industrial process. J. Loss Prev. Process Ind. 40, 471–478. Hollifield, B., Habibi, E., Pinto, J., 2011. Alarm Management: A Comprehensive Guide. International Society of Automation. ISA, 2009. ANSI/ISA-18.2: Management of Alarm Systems for the Process Industries. International Society of Automation. Izadi, I., Shah, S., Shook, D., Chen, T., 2009. An introduction to alarm analysis and design. In: 7th IFAC Symposium on Fault Detection, Supervision and Safety of Technical Processes, Spain, pp. 645–650. Janusz, M.E., Venkatasubramanian, V., 1991. Automatic generation of qualitative descriptions of process trends for fault detection and diagnosis. Eng. Appl. Artif. Intell. 4, 329–339. Keogh, E., Chu, S., Hart, D., Pazzani, M., 2001. An online algorithm for segmenting time series. In: Proceedings IEEE International Conference on Data Mining (ICDM 2001), California, USA, pp. 289–296. Kondaveeti, S.R., Shah, S., Izadi, I., 2009. Application of multivariate statistics for efficient alarm generation. In: 7th IFAC Symposium on Fault Detection, Supervision and Safety of Technical Processes, Spain, pp. 657–662. Maurya, M.R., Paritosh, P.K., Rengaswamy, R., Venkatasubramanian, V., 2010. A framework for on-line trend extraction and fault diagnosis. Eng. Appl. Artif. Intell. 23, 950–960. Melek, W.W., Lu, Z., Kapps, A., Fraser, W.D., 2005. Comparison of trend detection algorithms in the analysis of physiological time-series data. IEEE Trans. Biomed. Eng. 52, 639–651. Rengaswamy, R., Venkatasubramanian, V., 1995. A syntactic pattern-recognition approach for process monitoring and fault diagnosis. Eng. Appl. Artif. Intell. 8, 35–51. Rothenberg, D.H., 2009. Alarm Management for Process Control: A Best-Practice Guide for Design, Implementation, and Use of Industrial Alarm Systems. Momentum Press, New Jersey. Savitzky, A., Golay, M.J., 1964. Smoothing and differentiation of data by simplified least squares procedures. Anal. Chem. 36, 1627–1639. Vedam, H., Venkatasubramanian, V., Bhalodia, M., 1998. A b-spline based method for data compression, process monitoring and diagnosis. Comput. Chem. Eng. 22, S827–S830. Wang, J., Chen, T., 2014. An online method to remove chattering and repeating alarms based on alarm durations and intervals. Comput. Chem. Eng. 67, 43–52. Wang, J., Yang, F., Chen, T., Shah, S.L., 2016. An overview of industrial alarm systems: main causes for alarm overloading, research status, and open problems. IEEE Trans. Autom. Sci. Eng. 3, 1045–1061. Zang, H., Li, H., 2014. Optimization of process alarm thresholds: A multidimensional kernel density estimation approach. Process Saf. Prog. 33, 292–298. Zhu, J., Shu, Y., Zhao, J., Yang, F., 2014. A dynamic alarm management strategy for chemical process transitions. J. Loss Prev. Process Ind. 30, 207–218.