Sensitivity analysis and sensitivity-based design for linear alarm filters

Control Engineering Practice 70 (2018) 29–39

Contents lists available at ScienceDirect

Control Engineering Practice journal homepage: www.elsevier.com/locate/conengprac

Sensitivity analysis and sensitivity-based design for linear alarm filters✩ Ying Xiong *, Yindi Jing, Tongwen Chen Department of Electrical and Computer Engineering, University of Alberta, Edmonton, AB, Canada T6G 1H9

a r t i c l e

i n f o

Keywords: Alarm monitoring Linear filters Sensitivity analysis Sensitivity-based design

a b s t r a c t This paper conducts sensitivity analysis and sensitivity-based design for linear filter alarm monitoring systems. Based on a derivative-based local sensitivity measure, models are proposed to assess the sensitivity of the system detection errors to changes in the trip point and to uncertainties in the collected data. Then, analytical expressions are derived to quantitatively evaluate the sensitivity of a general linear alarm filter with unknown data distributions. Subsequently, a new sensitivity-based linear filter design method is formulated to minimize a weighted sum of the detection errors subject to upper bounds on the system sensitivities. Extensive simulations with both Gaussian and industrial data are conducted to verify the analytical results and to show trade-offs between the detection errors and sensitivities of linear filter alarm system. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction With increasingly high requirements on safety, reliability and efficiency in industrial plants, the design of effective alarm monitoring systems has become a heated topic during the last decade. Generally speaking, to improve the accuracy of alarm monitoring systems is to reduce the two types of errors: the false alarm rate (FAR) and the miss alarm rate (MAR) (Adnan, Izadi, & Chen, 2011). There are two main research directions related with this error reduction task. The first one is the design of alarm triggering mechanism, including delay-timers (Zang, Yang, & Huang, 2015), dead-bands (Naghoosi, Izadi, & Chen, 2011) and filters (Cheng, Izadi, & Chen, 2011). The second direction is to analyze the connections among alarms to eliminate sequential alarms and ascertain rationalization suggestions, using correlation analysis (Yang, Shah, Xiao, & Chen, 2012), similarity analysis (Ahmed, Izadi, Chen, Joe, & Burton, 2013) and sequence pattern analysis (Lai & Chen, 2015). Due to the heavy computation burden and implementation cost in the second direction above, research topics on improved alarm triggering mechanism, especially the advanced filter design with a level-crossing alarm generation mechanism, are still among the most promising ones to reduce nuisance and false alarms. Moreover, linear alarm filters are of particular interest since they are not only easy to implement but also effective in performance enhancement. In the literature, firstly, moving average (MA) filters and moving variance filters were designed for known Gaussian distributions in Izadi, Shah, Shook, Kondaveeti, and Chen (2009). Later, a general optimal design ✩

framework for linear and quadratic alarm filters (Cheng et al., 2011) was proposed to tackle data with other distributions. Based on the design framework in Cheng et al. (2011), MA filters were analytically proved to be optimal among linear filters under log-concave and symmetric data distributions (Cheng, Izadi, & Chen, 2013). The majority of existing results for linear alarm filters still rely on the assumption that the trip point and the probability density functions (PDFs) of data under normal and abnormal conditions are known precisely. However, in reality, uncertainties and changes are inevitable in both the trip point and the PDF estimation. Specifically, computational errors exist in the trip point setting. For example, in the linear filter design (Cheng et al., 2013), the solution is found with iterative numerical methods (Cheng et al., 2011), which offer an approximate solution instead of the exactly optimal one due to limited precisions and limited processing power. As for the uncertainties in the PDFs, strictly speaking, the exact PDFs of process signals are unknown. The PDFs are actually estimated based on the collected data and preassumed probability models. But measurement noise in data collection is ubiquitous. Uncertainties in both the trip point and the PDFs will propagate into filter design and alarm performance evaluation. Thus, to understand how uncertainties affect the performance of an alarm system, it is important to quantitatively measure the sensitivity of the system behavior with respect to the trip point change and the PDF variations.

A preliminary version of this paper was presented in the American Control Conference (ACC), Seattle, WA, USA, May 24–26, 2017.

* Corresponding author.

E-mail addresses: [email protected] (Y. Xiong), [email protected] (Y. Jing), [email protected] (T. Chen). https://doi.org/10.1016/j.conengprac.2017.09.013 Received 3 May 2017; Received in revised form 30 July 2017; Accepted 18 September 2017 0967-0661/© 2017 Elsevier Ltd. All rights reserved.

Y. Xiong et al.

Control Engineering Practice 70 (2018) 29–39

In the literature, sensitivity measures can be generally classified into two types: local sensitivity measures and global sensitivity measures (Borgonovo & Plischke, 2016). Global sensitivity measures are usually obtained using statistical properties (such as variance ratio Borgonovo & Plischke, 2016) based on assigned known density distributions (Saltelli, Ratto, Andres, Campolongo, Cariboni et al., 2008). On the other hand, pre-given PDFs are unnecessary for local sensitivity measures (Borgonovo & Plischke, 2016). Besides, the two types of measures have different applications. Local sensitivity measures are widely used in evaluating output variations due to slight changes in inputs, while global ones deal with input variations over their entire ranges of interest. With the help of proper sensitivity measures, a sensitivity-based design problem can be formulated to meet diverse engineering needs. For example, in mechanical engineering, derivativebased local sensitivity measures are used to optimize the feasibility of a welded beam in Gunawan and Azarm (2004). In control engineering, the sensitivity of a control system to disturbances is usually reflected by a global sensitivity measure namely the sensitivity transfer function. For example, the weighted ∞ -norm of the sensitivity function was minimized to achieve optimal robust feedback controller in Zames and Francis (1983), and the 2-norm of the sensitivity function quantified the robustness of the assigned pole to parameter perturbations in a dynamical system in Abdelaziz (2015). In this work, a local sensitivity measure is chosen instead of global ones. Because results obtained with global sensitivity measures highly rely on pre-assumed and limited types of distribution models; without known distributions, the calculation of global sensitivity measures is mathematically troublesome, if not impossible. However, distributions of the collected data are usually unknown and have relatively large diversity for different processes. Moreover, among existing local sensitivity measures, the derivative-based elasticity measure (Borgonovo & Plischke, 2016; Karnavas, Sanchez, & Bahill, 1993) is adopted and customized for linear alarm filters. Specifically, the inputs of the elasticity are chosen as the trip point variation and the PDF offsets caused by noise in the collected data. And the Kullback–Leibler divergence (KLD), a popular distance measure between probability distributions in information theory and statistics (Kullback & Leibler, 1951), is applied to measure the PDF offsets between the ideal noise-free data and the collected noisy data. The MAR and FAR are chosen as the outputs. Without any knowledge on the collected data distribution, sensitivities of the detection errors to the trip point and KLD are calculated based on the Gaussian kernel estimation. In addition, based on derived sensitivity expressions, a new linear filter design method is proposed, which takes both detection errors and sensitivity measures into account. The design procedure is modeled as a constrained minimization problem, where the objective is to minimize the weighted summation of detection errors within sensitivity constraints. A grid search method is used to locate the solution. Simulations on the Gaussian and industrial data are conducted to verify the effectiveness of the proposed sensitivity analysis and sensitivitybased design. Also, simulation results provide insightful observations in detection errors and the sensitivity of linear alarm filters. The rest of the paper is organized as follows. Backgrounds on alarm monitoring with linear filters and sensitivity analysis models are illustrated in Section 2. Then in Section 3, the proposed sensitivity models and the detailed calculation procedures are presented, together with the sensitivity-based linear filter design. In Section 4, simulation results on data sets from Gaussian distribution and from an industrial plant are presented and studied. Section 5 concludes this paper.

discrete-time signals 𝒚 = [𝑦(1), 𝑦(2), …], written as 𝒚 = 𝐹 (𝒙). In the alarm monitoring area, linear filters show significant advantages over non-linear filters (Cheng et al., 2011) in easy application and simple computation. Therefore, linear filters are considered in this work, with the following form: 𝑦(𝑘) =

𝑁−1 ∑

(1)

𝜃𝑖 𝑥(𝑘 − 𝑖),

𝑖=0

where 𝑁 is the length of the linear filter and 𝜃𝑖 ’s are the filter coefficients. Without loss of generality, the summation of all the filter ∑ coefficients is set to be 1, i.e., 𝑁−1 𝑖=0 𝜃𝑖 = 1. With the typical level crossing alarm generation mechanism (Izadi et al., 2009), the filtered signal 𝑦(𝑘) is compared with the trip point 𝑦𝑡𝑝 to trigger an alarm, i.e., an alarm will be raised if 𝑦(𝑘) > 𝑦𝑡𝑝 . In this work, collected data under the normal and abnormal modes are denoted as 𝒙𝑛 = [𝑎1 , 𝑎2 , … , 𝑎𝑙𝑛 ] with length 𝑙𝑛 and 𝒙𝑎𝑏 = [𝑏1 , 𝑏2 , … , 𝑏𝑙𝑎𝑏 ] with length 𝑙𝑎𝑏 respectively. It is assumed that 𝑎𝑖 ’s are independent and identically distributed (i.i.d.), 𝑏𝑗 ’s are i.i.d., and 𝑎𝑖 ’s are independent to 𝑏𝑗 ’s. Let 𝑓𝑋,𝑛 (⋅) and 𝑓𝑋,𝑎𝑏 (⋅) be the PDFs of 𝑎𝑖 and 𝑏𝑗 respectively. Also, let 𝒚𝑛 = 𝐹 (𝒙𝑛 ) and 𝒚𝑎𝑏 = 𝐹 (𝒙𝑎𝑏 ) denote the filtered signals, and 𝑓𝑌 ,𝑛 (⋅) and 𝑓𝑌 ,𝑎𝑏 (⋅) represent the common PDF of each element in 𝒚𝑛 and 𝒚𝑎𝑏 respectively. The FAR (MAR) of the alarm system is defined as the probability that an alarm is raised (missed) when the system is under the normal (abnormal) mode. The goal in conventional linear filter design is to find the optimal combination of filter coefficients and the trip point with the lowest weighted sum of FAR and MAR. With the trip point 𝑦𝑡𝑝 , the cost function can be calculated as: Jconv = 𝑐1 FAR + 𝑐2 MAR +∞

= 𝑐1

∫𝑦𝑡𝑝

𝑦𝑡𝑝

𝑓𝑌 ,𝑛 (𝑦)d𝑦 + 𝑐2

∫−∞

𝑓𝑌 ,𝑎𝑏 (𝑦)d𝑦,

(2)

where 𝑐1 and 𝑐2 are positive weights on the FAR and MAR. From (2), it is obvious that the FAR, the MAR and the cost function can be significantly affected by the accuracy of the data PDFs and the trip point. Uncertainties in the estimated PDFs and the trip point will lead to performance change. Thus, there exists a strong need to quantify the sensitivity to these changes. 2.2. Sensitivity definition and properties In this work, the derivative-based elasticity, a local sensitivity measure, is adopted because it eliminates unit differences of input and output variations (Borgonovo & Apostolakis, 2001). More rationale for this measure can be found in Borgonovo and Plischke (2016). The measure is defined as: / d𝑂 𝐼 𝛥𝑂 𝛥𝐼 = , (3) 𝑆𝑂𝐼 = lim 𝛥𝐼→0 𝑂 𝐼 d𝐼 𝑂 where 𝑂 stands for the output of the sensitivity model and 𝐼 represents the input with uncertainties. The ratio 𝐼∕𝑂 is used to scale the changes in the input and output over their current values to obtain a normalized dimensionless value. 𝑆𝑂𝐼 is naturally a sensitivity measure of 𝑂 to 𝐼 representing how much the output 𝑂 changes with the input 𝐼. To complete the sensitivity modeling, it is necessary to specify inputs and outputs that are not only relevant but also capable of representing practical industrial needs. According to Borgonovo and Plischke (2016), model outputs should be variables of most interests to the decision maker. So the performance indices, the FAR and MAR, are chosen as outputs for sensitivity analysis. As for the inputs, once the filter design is completed, candidates are the trip point and the estimated PDFs. The sensitivity of the detection error over the trip point is straightforward to model by using the definition in (3):

2. Linear alarm filters and sensitivity models 2.1. Linear alarm filters Let 𝒙 = [𝑥(1), 𝑥(2), …] be a vector of discrete-time signals. A filter is essentially an operator that processes 𝒙 to produce another vector 30

Y. Xiong et al. 𝑦

𝑡𝑝 𝑆FAR =

Control Engineering Practice 70 (2018) 29–39

dFAR 𝑦𝑡𝑝 , d𝑦𝑡𝑝 FAR

𝑦

𝑡𝑝 𝑆MAR =

dMAR 𝑦𝑡𝑝 . d𝑦𝑡𝑝 MAR

(4)

𝑓𝑋,𝑛 (𝑥) =

The sensitivity over the estimated PDFs, however, needs further modeling since the PDF itself is a function not a scalar. For alarm systems, estimated PDFs are obtained from the collected data sets with random noise. Thus, the uncertainty in estimated PDFs can be represented by the distance between the PDF estimated from the noisy data and the real PDF. In this work, the KLD is applied to measure this distance. Specifically, let 𝑓𝑌 𝜖,𝑛 (𝑦) be the estimated PDF using data without noise, and let 𝑓𝑌 ,𝑛 (𝑦) represent the estimated PDF of raw data with noise. The KLD from 𝑓𝑌 𝜖,𝑛 (𝑦) to 𝑓𝑌 ,𝑛 (𝑦) is defined as: +∞

𝐷𝑘𝑙 (𝑓𝑌 ,𝑛 (𝑦) ∥ 𝑓𝑌 𝜖,𝑛 (𝑦)) =

∫−∞

𝑓𝑌 ,𝑛 (𝑦) ln

𝑓𝑌 ,𝑛 (𝑦) 𝑓𝑌 𝜖,𝑛 (𝑦)

d𝑦.

= lim

𝛥𝐷𝑘𝑙 →0

𝛥FAR 𝐷𝑘𝑙 , 𝛥𝐷𝑘𝑙 FAR

𝐷𝑘𝑙 𝑆MAR

= lim

𝛥𝐷𝑘𝑙 →0

𝛥MAR 𝐷𝑘𝑙 . 𝛥𝐷𝑘𝑙 MAR

(9)

The characteristic function of the estimated PDF under the normal mode, denoted with 𝜓𝑋,𝑛 (⋅), can be calculated as: +∞

𝜓𝑋,𝑛 (𝑡) =

∫−∞

𝑓𝑋,𝑛 (𝑥)e𝑗𝑥𝑡 d𝑥.

(10)

With the help of (9) and properties of the characteristic function under Gaussian distribution (Walck, 2007), it can be derived that 𝜓𝑋,𝑛 (𝑡) =

(5)

𝑙𝑛 1 ∑ − 21 𝑡2 +𝑗𝑡𝑎𝑘 e . 𝑙𝑛 𝑘=1

(11)

As 𝑎𝑖 ’s are i.i.d, with the linear filter in (1), the characteristic function of the filtered signal 𝜓𝑌 ,𝑛 (⋅) is derived as:

According to the sensitivity definition in (3), the sensitivities of the FAR and the MAR over the KLD are: 𝐷𝑘𝑙 𝑆FAR

𝑙𝑛 2 1 1 ∑ − (𝑥−𝑎2 𝑖 ) . e √ 𝑙𝑛 2𝜋 𝑖=1

𝜓𝑌 ,𝑛 (𝑡) =

(6)

𝑁−1 ∏

𝜓𝑋,𝑛 (𝜃𝑖 𝑡)

𝑖=0

( )𝑁 𝑁−1 𝑙𝑛 ∏ ∑ − 1 𝑡2 𝜃2 +𝑗𝑡𝑎 𝜃 1 𝑘 𝑖. = e 2 𝑖 𝑙𝑛 𝑖=0 𝑘=1

Next, properties of the sensitivity measures and insights they can provided are discussed. The elasticity of the FAR over the trip point, 𝑦𝑡𝑝 i.e., 𝑆FAR defined in (3), is used as an example.

(12)

To simplify the notation, let 𝒳𝑛𝑁 = {𝒙 ∣ 𝒙 = [𝑥1 , … , 𝑥𝑖 , … , 𝑥𝑁 ], 𝑥𝑖 ∈ 𝒙𝑛 }, then (12) can be rewritten as: ( )𝑁 1 2 2 ∑ 𝑇 1 𝜓𝑌 ,𝑛 (𝑡) = e− 2 𝑡 ‖𝜽‖ e𝑗𝑡𝜽𝒙 , (13) 𝑙𝑛 𝑁

1. Sign of the sensitivity measure: If the scaling term in the sensitivity measure (i.e., 𝑦𝑡𝑝 ∕FAR in (4)) is positive, the positive sign 𝑦𝑡𝑝 of the sensitivity measure 𝑆FAR means that the FAR increases with the trip point, and vise versa. Thus, the sign indicates the direction of the FAR changes with 𝑦𝑡𝑝 . If the scaling term is negative, the opposite holds. 2. Magnitude of the sensitivity measure: The magnitude is the percentage change ratio of the FAR over 𝑦𝑡𝑝 . A larger magnitude means that a more significant change occurs in the FAR with the same change in 𝑦𝑡𝑝 . For good robustness, a small sensitivity magnitude is desired. 3. Information from different sensitivity measures: Since different sensitivity measures reveal model behavior from different aspects, each measure deserves careful attention in terms of the sign, magnitude and trend. So, it is natural to analyze them together. A direct way is to define an overall sensitivity measure as: ∑ 𝑆overall = 𝑐𝑖,𝑗 𝑆𝑖𝑗 , (7)

𝒙∈𝒳𝑛

where ‖ ⋅ ‖ denotes the Euclidean norm of a vector. The PDF of the filtered data from the normal mode can thus be calculated as follows: +∞

1 e−𝑗𝑦𝑡 𝜓𝑌 ,𝑛 (𝑡)d𝑡 2𝜋 ∫−∞ ( )𝑁 ∑ − (𝑦−𝜽𝒙𝑇 )2 1 1 = 𝐞 2‖𝜽‖2 . √ 𝑙𝑛 𝑁 2𝜋‖𝜽‖

𝑓𝑌 ,𝑛 (𝑦) =

(14)

𝒙𝒏 ∈𝒳𝑛

According to (2), the FAR of the linear alarm filter can be calculated to be: ( )𝑁 +∞ ∑ − (𝑦−𝜽𝒙𝑇 )2 1 1 (15) FAR = e 2‖𝜽‖2 d𝑦. √ 𝑙𝑛 ∫𝑦𝑡𝑝 2𝜋‖𝜽‖ 𝒙∈𝒳𝑛𝑁 From the definition in (4), the sensitivity of the FAR over the trip point is eventually obtained as:

𝑖,𝑗

where 𝑖 ∈ {FAR, MAR} , 𝑗 ∈ {𝑦𝑡𝑝 , KLD} and 𝑐𝑖,𝑗 ’s are constants.

𝑇 2 ( )𝑁 ∑ − (𝑦𝑡𝑝 −𝜽𝒙 ) 𝑦𝑡𝑝 𝑦𝑡𝑝 1 𝑆FAR =− e 2‖𝜽‖2 . √ 𝑙𝑛 FAR 2𝜋‖𝜽‖ 𝑁

(16)

𝒙∈𝒳𝑛

3. Sensitivity analysis and sensitivity-based design for linear alarm filters

Via similar analysis, the sensitivity of the MAR over the trip point can be obtained as: 𝑇 2 ( )𝑁 ∑ − ( 𝑦𝑡𝑝 −𝜽𝒙 ) 𝑦𝑡𝑝 𝑦𝑡𝑝 1 2‖𝜽‖2 𝑆MAR = e , (17) √ 𝑙𝑎𝑏 MAR 2𝜋‖𝜽‖ 𝒙∈𝒳 𝑁

In this section, calculations of the sensitivity measures defined in (3) are provided for a generic linear alarm filter. The collected data under the normal and abnormal modes are 𝒙𝑛 with length 𝑙𝑛 and 𝒙𝑎𝑏 with length 𝑙𝑎𝑏 as introduced in Section 2. The coefficients of the linear filter 𝜽 = [𝜃0 , 𝜃1 , … , 𝜃𝑁−1 ] are assumed to be known though arbitrary.

𝑎𝑏

𝑁 = {𝒙 ∣ 𝒙 = [𝑥 , … , 𝑥 , … , 𝑥 ], 𝑥 ∈ 𝒙 }. where 𝒳𝑎𝑏 1 𝑖 𝑁 𝑖 𝑎𝑏

3.1. Sensitivity over trip point

3.2. Sensitivity over noise caused PDF offsets

Since the PDFs of the process signals are unknown, the Gaussian kernel function is applied to the PDF estimation using the collected data from the normal mode:

In this subsection, sensitivities of the performance over the PDF offsets are derived. In order to mimic the uncertainties in the estimated PDFs, the Gaussian noise contained in the collected data are denoted by 𝝐 𝒏 = [𝜖1,𝑛 , … , 𝜖𝑙𝑛 ,𝑛 ] and 𝝐 𝒂𝒃 = [𝜖1,𝑎𝑏 , … , 𝜖𝑙𝑎𝑏 ,𝑎𝑏 ] respectively under the normal and abnormal conditions. For each data set, the noise from different time instants are assumed to be i.i.d. with the same distribution 𝜖𝑛 ∼  (𝑚, 𝜎 2 ) and 𝜖𝑎𝑏 ∼  (𝑚, 𝜎 2 ). Thus, data before filtering is modeled as

𝑓𝑋,𝑛 (𝑥) =

𝑙𝑛 1 ∑ 𝐾(𝑥 − 𝑎𝑖 ), 𝑙𝑛 𝑖=1

(8)

where 𝐾(𝑥) is the Kernel function. In this work, the Kernel function is chosen as the standard Gaussian function (Parzen, 1962). Thus, 31

Y. Xiong et al.

Control Engineering Practice 70 (2018) 29–39

and 𝒙𝒂𝒃,𝒇 = 𝒙𝒂𝒃 + 𝝐 𝒂𝒃

𝒙𝒏,𝒇 = 𝒙𝒏 + 𝝐 𝒏

𝜕FAR = 𝜕𝑚

under the normal and abnormal modes respectively. It is reasonable to 𝜕FAR = 𝜕𝜎

assume that noise is independent of the collected data. Thus, under the

( (

normal mode, the output of the linear filter in (1) can be written as: 𝑦𝑛 (𝑘) =

𝑁−1 ∑

[ ] 𝜃𝑖 𝑥𝑛 (𝑘 − 𝑖) + 𝜖𝑛 (𝑘 − 𝑖) .

=

( ) − 12 (1+𝜎 2 )𝑡2 +𝑗𝑡 𝜽𝒙𝑇 +𝑚

e

,

(𝑦𝑡𝑝 − 𝜽𝒙𝑇 − 𝑚)e

∑ (𝑦𝑡𝑝 − 𝜽𝒙𝑇 − 𝑚)

(19)

(1 + 𝜎 2 )‖𝜃‖2

𝒙∈𝒳𝑛 𝑇

tion. Accordingly, the PDF of the filtered data with noises is obtained

(20)

𝒙∈𝒳𝑛

𝐷

𝑘𝑙 𝑆𝑀𝐴𝑅 =

Thus, the FAR can be expressed as follows: )𝑁

+∞

𝑓𝑌 𝜖,𝑛 (𝑦)d𝑦 ∫𝑦𝑡𝑝 ) ( ( )𝑁 ∑ 𝑦𝑡𝑝 − 𝜽𝒙𝑇 − 𝑚 1 1 1 , = − erf √ 2 2 𝑙𝑛 2(1 + 𝜎 2 )‖𝜽‖ 𝑁

(21)

As is shown in (20), both the mean and standard deviation of the

𝜕MAR = 𝜕𝜎

Gaussian noise offset the PDF of the filtered data. One method to

(

quantify the offset is to further introduce the mean and variance as two inputs to capture the sensitivity of the FAR over the uncertainty

×

(25)

] d𝑦.

(26)

𝐷𝑘𝑙 , MAR

(27)

∑

(𝑦𝑡𝑝 − 𝜽𝒙𝑇 − 𝑚)e

−

(𝑦𝑡𝑝 −𝜽𝒙𝑇 −𝑚)2 2(1+𝜎 2 )‖𝜽‖2

,

(28)

,

(29)

𝑁 𝒙∈𝒳𝑎𝑏

1 𝑙𝑎𝑏

)𝑁

∑

𝜎 1 √ (1 + 𝜎 2 ) 2𝜋(1 + 𝜎 2 )‖𝜽‖

(𝑦𝑡𝑝 − 𝜽𝒙𝑇 − 𝑚)e

−

(𝑦𝑡𝑝 −𝜽𝒙𝑇 −𝑚)2 2(1+𝜎 2 )‖𝜽‖2

( )𝑁 +∞ 𝑓 𝜕𝐷𝑘𝑙 𝑌𝑎𝑏 (𝑦) 1 1 = √ ∫−∞ 𝑓𝑌 (𝑦) 𝑙𝑎𝑏 𝜕𝑚 2𝜋(1 + 𝜎 2 )‖𝜽‖ 𝜖,𝑎𝑏

variance of the noise are considered separately (Xiong, Jing, & Chen, 2017). In this paper, a different method is taken with the help of the KLD to model the PDF estimation offset caused by both the mean and

×

standard deviation of the noise jointly.

∑ (𝑦𝑡𝑝 − 𝜽𝒙𝑇 − 𝑚) 𝑁 𝒙∈𝒳𝑎𝑏

From (6), the sensitivity of the FAR over the KLD can be expanded

(1 + 𝜎 2 )‖𝜃‖2

e

−

(𝑦𝑡𝑝 −𝜽𝒙𝑇 −𝑚)2 2(1+𝜎 2 )‖𝜽‖2

d𝑦,

(30)

+∞ ( )𝑁 𝑓𝑌𝑎𝑏 (𝑦) 𝜕𝐷𝑘𝑙 1 𝜎 = √ ∫−∞ 𝜕𝜎 𝑙𝑛 𝑓 2 2 2𝜋(1 + 𝜎 )(1 + 𝜎 )‖𝜽‖ 𝑌𝜖,𝑎𝑏 (𝑦) [ (𝑦𝑡𝑝 −𝜽𝒙𝑇 −𝑚)2 ∑ − 1 × (𝑦𝑡𝑝−𝜽𝒙𝑇 −𝑚)2 e 2(1+𝜎 2 )‖𝜽‖2 (1 + 𝜎 2 )‖𝜽‖2 𝑁

as follows according to the total differential theorem: 𝐷𝑘𝑙 , FAR

d𝑦,

𝑁 𝒙∈𝒳𝑎𝑏

in the PDF estimation. In this way, the sensitivities over the mean and

𝜕FAR d𝑚 + 𝜕FAR d𝜎 𝜕𝑚 𝜕𝜎 𝜕𝐷𝑘𝑙 d𝑚→0,d𝜎→0 𝜕𝐷𝑘𝑙 d𝑚 + 𝜕𝜎 d𝜎 d𝑚=𝜆d𝜎 𝜕𝑚

(𝑦𝑡𝑝 −𝜽𝒙𝑇 −𝑚)2 2(1+𝜎 2 )‖𝜽‖2

where the partial derivative terms are shown in (28) to (31). ( )𝑁 𝜕MAR 1 1 =− √ 𝜕𝑚 𝑙𝑎𝑏 2𝜋(1 + 𝜎 2 )‖𝜽‖

where erf(⋅) is the error function.

lim

−

𝜕MAR d𝑚 + 𝜕MAR d𝜎 𝜕𝑚 𝜕𝜎 𝜕𝐷𝑘𝑙 𝜕𝐷𝑘𝑙 d𝑚→0,d𝜎→0 d𝑚 + 𝜕𝜎 d𝜎 d𝑚=𝜆d𝜎 𝜕𝑚

×

𝐷

e

lim

𝒙∈𝒳𝑛

𝑘𝑙 𝑆FAR =

(24)

,

Via tedious but straightforward calculations, each term in (22) can be obtained in (23) to (26) using (5) and (21). The details are omitted due to page limit. Similarly, the sensitivity of the MAR over the KLD can be calculated as follows:

+∞

1 𝑙𝑛

(𝑦𝑡𝑝 −𝜽𝒙𝑇 −𝑚)2 2(1+𝜎 2 )‖𝜽‖2

𝒙∈𝒳𝑛𝑁

as:

(

2

∑ − (𝑦𝑡𝑝 −𝜽𝒙 −𝑚) − e 2(1+𝜎 2 )‖𝜽‖2

where 𝜓𝑒 (⋅) is the characteristic function of the Gaussian noise distribu-

FAR =

−

+∞ ( )𝑁 𝑓𝑌𝑛 (𝑦) 𝜕𝐷𝑘𝑙 1 𝜎 = √ ∫−∞ 𝜕𝜎 𝑙𝑛 𝑓 2 2 2𝜋(1 + 𝜎 )(1 + 𝜎 )‖𝜽‖ 𝑌𝜖,𝑛 (𝑦) [ (𝑦𝑡𝑝 −𝜽𝒙𝑇 −𝑚)2 ∑ − 1 (𝑦𝑡𝑝− 𝜽𝒙𝑇 −𝑚)2 e 2(1+𝜎 2 )‖𝜽‖2 × (1 + 𝜎 2 )‖𝜽‖2 𝑁

𝒙∈𝒳𝑛𝑁

1 e−𝑗𝑦𝑡 𝜓𝑌𝜖 ,𝑛 (𝑡)d𝑡 𝑓𝑌 𝜖,𝑛 (𝑦) = 2𝜋 ∫−∞ ( )𝑁 (𝑦−𝜽𝒙𝑇 −𝑚)2 ∑ − 1 1 = e 2(1+𝜎 2 )‖𝜽‖2 . √ 𝑙𝑛 2𝜋(1 + 𝜎 2 )‖𝜽‖ 𝑁

(23)

𝜎 1 √ (1 + 𝜎 2 ) 2𝜋(1 + 𝜎 2 )‖𝜽‖

𝒙∈𝒳𝑛𝑁

𝑖=0

∑

)𝑁

∑

×

which can be calculated as: (𝑁−1 ) (𝑁−1 ) ∏ ∏ 𝜓𝑌 𝜖,𝑛 (𝑡) = 𝜓𝑋,𝑛 (𝜃𝑖 𝑡) 𝜓𝑒 (𝜃𝑖 𝑡) )𝑁

2

+∞ 𝑓 (𝑦) ( )𝑁 𝜕𝐷𝑘𝑙 𝑌𝑛 1 1 = √ ∫−∞ 𝑓𝑌 (𝑦) 𝑙𝑛 𝜕𝑚 2𝜋(1 + 𝜎 2 )‖𝜽‖ 𝜖,𝑛

Thus, by following similar steps in the previous subsection, the

𝑖=0

𝑇

∑ − (𝑦𝑡𝑝 −𝜽𝒙 −𝑚) 1 e 2(1+𝜎 2 )‖𝜽‖2 , √ 2𝜋(1 + 𝜎 2 )‖𝜽‖ 𝒙∈𝒳𝑛𝑁

𝒙∈𝒳𝑛𝑁

(18)

characteristic function of 𝑦𝑛 with Gaussian noises is denoted by 𝜓𝑌 𝜖,𝑛 (⋅),

1 𝑙𝑛

1 𝑙𝑛

×

𝑖=0

(

1 𝑙𝑛

)𝑁

(22)

𝒙∈𝒳𝑎𝑏

where the condition d𝑚 = 𝜆d𝜎 is to specify a linear relationship between −

d𝑚 and d𝜎. This condition is necessary to guarantee that the limit is well

∑

e

(𝑦𝑡𝑝 −𝜽𝒙𝑇 −𝑚)2 − 2(1+𝜎 2 )‖𝜽‖2

] d𝑦.

(31)

𝑁 𝒙∈𝒳𝑎𝑏

defined. The coefficient 𝜆 can be seen as a weighting factor to quantify the step size ratio of the mean over that of the standard deviation. A larger 𝜆 means that a heavier weight is put on the standard deviation

3.3. Sensitivity-based linear filter design

change than the mean change. In applications, 𝜆 should be tuned larger in situations where the system is more vulnerable to the variance of

Linear alarm filter design should consider not only the detection errors, but also the tolerance to system uncertainties, namely, the sensitivities to both the trip point changes and noise. In this paper, a new

noise than to the mean. In this paper, 𝜆 = 1 is chosen for simplicity, i.e., the mean and standard deviation changes are treated equally. 32

Y. Xiong et al.

Control Engineering Practice 70 (2018) 29–39

linear filter design is proposed to achieve the least detection errors given sensitivity requirements. Thus, the weighted sum of performance errors shown in (2) is taken as the objective function, while the sensitivity measures are used as the constraints in the filter design problem. There are two constraints involved. One is on the weighted sum of sensitivity measures of the FAR and MAR over the trip point. Since the trip point is one of the optimization variables, the formulation of this constraint is straightforward. Another is on the weighted sum of the sensitivities of the FAR and MAR over the KLD. But sensitivity measures with respect to the KLD need to be further modeled, because the mean and standard deviation of noise are system parameters that can take any value in a continuous range. In this work, the average sensitivities over certain mean and variance ranges are used. Assume that the mean and standard deviation of noise have uniformly distributions respectively under the normal mode, i.e., 𝑚𝑛 ∼ 𝑈 (𝑚𝑛1 , 𝑚𝑛2 ) and 𝜎𝑛 ∼ 𝑈 (𝜎𝑛1 , 𝜎𝑛2 ). The average sensitivity of the FAR over the KLD can be obtained as: 𝐷𝑘𝑙 𝑆̄𝐹 𝐴𝑅 =

(𝑚𝑛2

𝑚

𝜎

1

1

𝑛2 𝑛2 1 𝐷 𝑆 𝑘𝑙 (𝑚, 𝜎)d𝑚d𝜎. − 𝑚𝑛1)(𝜎𝑛2 − 𝜎𝑛1∫ ) 𝑚𝑛 ∫𝜎𝑛 𝐹 𝐴𝑅

reach an appropriate balance between them, firstly, the ranges of 𝑐3 and 𝑐4 should be calculated. Define 𝑐3max =

𝑐4max =

𝑚

(32)

𝜎

𝑎𝑏2 𝑎𝑏2 1 𝐷 𝑆 𝑘𝑙 (𝑚,𝜎)d𝑚d𝜎. (𝑚𝑎𝑏2−𝑚𝑎𝑏1)(𝜎𝑎𝑏2 −𝜎𝑎𝑏1)∫𝑚𝑎𝑏 ∫𝜎𝑎𝑏 𝐹 𝐴𝑅 1

(34)

𝜃0 ,…,𝜃𝑁−1 ,𝑦𝑡𝑝

𝑦

𝑦

𝐷

𝐷

𝑘𝑙 𝑘𝑙 𝑐1 |𝑆̄𝐹 𝐴𝑅 | + 𝑐2 |𝑆̄𝑀𝐴𝑅 | ≤ 𝑈4 ,

(35) (36)

where 𝑈3 and 𝑈4 are upper bound values of the sensitivities constraints over the trip point and over the KLD, respectively. They are set as: ( ) 𝑦 𝑦𝑡𝑝 𝑈3 = 𝑐 3 min 𝑐1 |𝑆𝐹𝑡𝑝𝐴𝑅 | + 𝑐2 |𝑆𝑀𝐴𝑅 | , 𝜃0 ,…,𝜃𝑁−1 ,𝑦𝑡𝑝

𝑈4 = 𝑐 4

min

𝜃0 ,…,𝜃𝑁−1 ,𝑦𝑡𝑝

𝑦

𝐷

𝐷

𝐷

𝐷

𝑘𝑙 𝑘𝑙 max𝜃0 ,…,𝜃𝑁−1 ,𝑦𝑡𝑝 (𝑐1 𝑆̄𝐹 𝐴𝑅 + 𝑐2 𝑆̄𝑀𝐴𝑅 ) 𝑘𝑙 𝑘𝑙 min𝜃0 ,…,𝜃𝑁−1 ,𝑦𝑡𝑝 (𝑐1 𝑆̄𝐹 𝐴𝑅 + 𝑐2 𝑆̄𝑀𝐴𝑅 )

,

,

To solve the constrained optimization problem in (34), a grid search method shown in Algorithm 1 is used. In the grid generation, 𝑁 − 1 dimensions of the hypercube [−1, 1]𝑁−1 are considered for the filter coefficients, since the summation of the filter coefficients is normalized to be 1. Let 𝑠 be the step size of the grid search in each coefficient dimension. There are (𝑝𝑁−1 ) candidate vectors for the filter coefficients, where 𝑝 = 2∕𝑠. Another dimension of the grid is the trip point. Let its range be [𝑦1 , 𝑦2 ] and the step size be 𝑦𝑠 . There are (𝑦2 − 𝑦1 )∕𝑦𝑠 candidate values for the trip point. Firstly, for every possible filter vector, the characteristic functions and the PDFs for both normal and abnormal data sets are calculated. Then, for every possible trip point value, the four sensitivity measures are calculated. In this step, sensitivities of the FAR and MAR over the KLD are calculated separately. Precisely, the mean sensitivity of the FAR over the KLD is calculated with (32) and that of the MAR is calculated with (33). Following that, two pairs of weighted summations are compared with the corresponding sensitivity requirements in (35) and (36). If both requirements are satisfied, the weighted detection error performance is calculated, the filter coefficient vector is saved, together with its corresponding feasible trip point range. For every filter vector, the trip point within the feasible range that with the minimum weighted sum of the detection errors and the corresponding optimal trip point are found and saved. By searching over all saved filter coefficients, one with the minimum weighted summation of detection errors is discovered.

1

𝑡𝑝 subject to 𝑐1 |𝑆𝐹𝑡𝑝𝐴𝑅 | + 𝑐2 |𝑆𝑀𝐴𝑅 | ≤ 𝑈3 ,

𝑦

𝑡𝑝 min𝜃0 ,…,𝜃𝑁−1 ,𝑦𝑡𝑝 (𝑐1 𝑆𝐹𝑡𝑝𝐴𝑅 + 𝑐2 𝑆𝑀𝐴𝑅 )

𝜃0 ,…,𝜃𝑁−1 ,𝑦𝑡𝑝

(33)

In practice, the ranges of the mean and variance can be defined from the knowledge of the disturbance or noise. Thus, the new sensitivity-based filter design problem is formulated as:

argmin 𝑐1 FAR + 𝑐2 MAR

𝑦

which are the maximum values of 𝑐3 and 𝑐4 for the constraints to be functional. Thus, the ranges for 𝑐3 and 𝑐4 are 𝑐3 ∈ [1, 𝑐3max ] and 𝑐4 ∈ [1, 𝑐4max ] respectively. Secondly, different sensitivity constraint levels can be described by the corresponding percentages over the entire coefficient range denoted by 𝜇3 and 𝜇4 . Then, the two coefficients should be set as 𝑐3 = 1 + 𝜇3 (𝑐3max − 1) and 𝑐4 = 1 + 𝜇4 (𝑐4max − 1). The sensitivity levels over the trip point and the KLD can be adjusted by 𝜇3 and 𝜇4 . If the sensitivity requirement over the trip point is higher than that over the KLD, 𝜇3 should be set to be smaller than 𝜇4 , and vise versa. When the sensitivity requirement over the trip point is treated as equal to that over the KLD, then 𝜇3 should be set to be equal to 𝜇4 . Accordingly, given percentage values of 𝜇3 and 𝜇4 , 𝑈3 and 𝑈4 can be further calculated as: ( ) 𝑦 𝑦𝑡𝑝 𝑈3 =[1 + 𝜇3 (𝑐3max − 1)] min 𝑐1 |𝑆𝐹𝑡𝑝𝐴𝑅 | + 𝑐2 |𝑆𝑀𝐴𝑅 | , (37) 𝜃0 ,…,𝜃𝑁−1 ,𝑦𝑡𝑝 ( ) 𝐷𝑘𝑙 𝐷𝑘𝑙 𝑈4 =[1 + 𝜇4 (𝑐4max − 1)] min 𝑐1 |𝑆̄𝐹 𝐴𝑅 | + 𝑐2 |𝑆̄𝑀𝐴𝑅 | . (38)

Similarly, by denoting the uniform distributions for the noise mean and the standard deviation in the abnormal mode as 𝑚𝑎𝑏 ∼ 𝑈 (𝑚𝑎𝑏1 , 𝑚𝑎𝑏2 ) and 𝜎𝑎𝑏 ∼ 𝑈 (𝜎𝑎𝑏1 , 𝜎𝑎𝑏2 ), the average sensitivity of the MAR over the KLD can be expressed as: 𝐷𝑘𝑙 𝑆̄𝑀𝐴𝑅 =

𝑦

𝑡𝑝 max𝜃0 ,…,𝜃𝑁−1 ,𝑦𝑡𝑝 (𝑐1 𝑆𝐹𝑡𝑝𝐴𝑅 + 𝑐2 𝑆𝑀𝐴𝑅 )

( ) 𝐷𝑘𝑙 𝐷𝑘𝑙 𝑐1 |𝑆̄𝐹 𝐴𝑅 | + 𝑐2 |𝑆̄𝑀𝐴𝑅 | ,

where 𝑐3 and 𝑐4 are coefficients to adjust the constraints. In this model, the upper bounds depend on the minimum weighted sum of the sensitivity values. Notice that 𝑐3 ≥ 1 and 𝑐4 ≥ 1 are necessary for the feasibility of the optimization. In (34), the 𝑐𝑖 ’s are weighting coefficients. They should be set and adjusted to meet practical engineering needs. Specifically, the weight ratio 𝑐1 ∕𝑐2 should be larger if the cost for false alarms is higher than that for miss alarms. Since the FAR and MAR are scaled by 𝑐1 ad 𝑐2 , their corresponding sensitivities are added up with the same weights 𝑐1 and 𝑐2 . Besides, the coefficients 𝑐3 and 𝑐4 are directly related to the sensitivity requirements. Smaller coefficients 𝑐3 and 𝑐4 indicate stronger sensitivity requirement on the corresponding uncertainties. However, the ratio 𝑐3 ∕𝑐4 does not have straightforwardly physical meaning unlike 𝑐1 ∕𝑐2 , since the sensitivity over the trip point can take a totally different scale with that over the KLD. For example, when coefficients 𝑐3 and 𝑐4 are set to be equal, it is likely that one of the sensitivity constraint is not effective anymore, which means that an arbitrary linear filter can meet this sensitivity requirement; while another constraint is still too strong to allow any feasible solution. In order to quantify strength levels of two sensitivity constraints and to

4. Case study In this section, simulation results based on data sets from Gaussian distribution and from an industrial plant are shown to verify the analytical results and perform the proposed filter design procedures in Section 3. Two widely used linear filters, the MA filter and the LW filter (Tsay, 2005), are considered and compared with the proposed optimal filter. The MA filter has equal coefficients, shown as: 𝑦(𝑘) =

𝑁−1 1 ∑ 𝑥(𝑘 − 𝑖); 𝑁 𝑖=0

(39)

while the LW filter has ascending coefficients, expressed as: 𝑦(𝑘) = 33

𝑁−1 1 ∑ (𝑖 + 1)𝑥(𝑘 − 𝑖). 𝑁 2 𝑖=0

(40)

Y. Xiong et al.

Control Engineering Practice 70 (2018) 29–39 Table 1 Sensitivity values at the optimal trip points for MA and LW filters under Gaussian distributions.

Algorithm 1: Sensitivity-based filter design.

1 2

3 4 5 6 7 8 9 10 11 12

13 14 15 16 17

18

Input: Data collected under normal and abnormal condition: 𝒙𝑛 , 𝒙𝑎𝑏 Output: Solution of the optimization problem (34) (𝜽∗ and 𝑦∗𝑡𝑝 Initialization: Generate (𝑁 − 1) × 𝑝 candidate vectors 𝜽 with 𝜃𝑖 ∈ [−1, 1]; 𝑖 = 1; ← Save the number of feasible filters. for 𝑗 ← 1 to (𝑁 − 1) × 𝑝 do For the 𝑗 th candidate of the filter coefficient vector 𝜽𝑗 , calculate 𝑓𝑌𝑛 (𝑦) and 𝑓𝑌 𝜖,𝑛 (𝑦) by (14) and (19); Similarly, obtain 𝑓𝑌𝑎𝑏 (𝑦) and 𝑓𝑌 𝜖,𝑎𝑏 (𝑦); for 𝑘 ← 1 to (𝑦2 − 𝑦1 ) × 𝑦𝑠 do Set 𝑦𝑡𝑝 to be the 𝑘 th candidate of the trip point;; 𝑦 𝐷 Calculate 𝑆 𝑡𝑝 and 𝑆̄ 𝑘𝑙 by (16) and (32); 𝐹 𝐴𝑅 𝑦

𝑦

𝑦

𝑦

𝑦

Filter

𝑡𝑝 𝑆FAR

𝑡𝑝 𝑆MAR

𝑡𝑝 𝑡𝑝 | | + |𝑆MAR |𝑆FAR

MA Filter|𝑦𝑡𝑝 = 7.89 LW Filter|𝑦𝑡𝑝 = 4.76

-10.95 -9.48

7.20 5.55

18.15 15.03

Via straightforward calculations according to (4), the sensitivity of the FAR and MAR over the trip point can be calculated as 𝑦

𝑡𝑝 𝑆FAR =−

𝐹 𝐴𝑅 𝐷

𝑡𝑝 𝑘𝑙 Calculate 𝑆𝑀𝐴𝑅 and 𝑆̄𝑀𝐴𝑅 by (17) and (33); if (35) and (36) hold then 𝜣(𝑖, ∶) = 𝜽𝑗 ; 𝒀 𝑡𝑝 (𝑖, 𝑘) = 𝑦𝑡𝑝 ; 𝐽 (𝑖, 𝑘) = 𝑐1 FAR + 𝑐2 MAR ; [𝐽𝑚𝑖𝑛 (𝑖), 𝑚] = min(𝐽 (𝑖, ∶)); ← Save the minimum cost function for each filter. 𝑀(𝑖) = 𝑚; ← Save the optimal trip point for each filter. 𝑖 = 𝑖 + 1; close; end end [𝐽𝑚𝑖𝑛 , 𝑛] = min(𝐽𝑚𝑖𝑛 (∶)); ← Find the minimum cost function in all feasible filters; Return 𝜽∗ = 𝜣(𝑛, ∶) and 𝑦∗𝑡𝑝 = 𝒀 𝑡𝑝 (𝑛, 𝑀(𝑛)).

𝑦

𝑡𝑝 𝑆MAR =

− 𝑦𝑡𝑝 e √ FAR 2𝜋𝜎𝑛 ‖𝜽‖

(𝑦𝑡𝑝 −𝑚𝑛 )2 2‖𝜽‖2 𝜎𝑛2

− 𝑦𝑡𝑝 e √ MAR 2𝜋𝜎𝑎𝑏 ‖𝜽‖

,

(𝑦𝑡𝑝 −𝑚𝑎𝑏 )2 2‖𝜽‖2 𝜎 2 𝑎𝑏

.

Firstly, it can be observed that the derived sensitivity formulas tightly match the theoretical ones for most of the parameter range. The small differences in Fig. 1 are caused by errors in the PDF estimation with the Gaussian kernel method. The gap can be further reduced by careful tuning of the bandwidth in the kernel function (8) (Hall, Sheather, Jones, & Marron, 1991). This verifies the effectiveness of the proposed sensitivity analysis under unknown distribution cases. Secondly, properties of sensitivity measures in Section 2 can be used to understand sensitivities of the FAR and MAR. According to Fig. 1(a), the sensitivity of the FAR has a downward trend with a negative sign. This means that the FAR decreases when the trip point grows. And the increasing magnitude means that the FAR gets more sensitive as the trip point rises. Similar analysis is applicable for the MAR in Fig. 1(b). The magnitudes of the sensitivities show that the MAR is less sensitive than the FAR with respect to the trip point change. In addition, the sensitivities of the two filters are compared. When the objective function in (2) has equal weights, i.e., 𝑐1 = 𝑐2 , the optimal trip points for the MA filter and the LW filter are 𝑦MA = 7.89 and 𝑡𝑝 𝑦LW = 4.76 respectively. The sensitivity values of the FAR and MAR 𝑡𝑝 for the two filters at their optimal trip points are listed in Table 1. The derived results in (16) and (17) are used for the calculations. It can be seen that the LW filter is less sensitive to the trip point variation. Moreover, the trend of sensitivity changing with different filter orders can be observed from Fig. 2. Specifically, in Figs. 2(a) and 2(b), the weighted summation of the sensitivity magnitude of the FAR and MAR shown in (35) is calculated with 𝑐1 = 𝑐2 = 1 for the MA and LW filters, respectively. Both figures show that the magnitude of the sensitivity measures become larger with higher filter orders. This means both the MA and LW filters tend to be more sensitive to the trip point change when the order increases.

The LW filter reduces both the mean and variance of data, while the MA filter only changes the variance and keeps the mean unchanged. The variance reduction effect of the LW filter is stronger than the MA filter. In this sense, the LW filter can achieve lower FAR and MAR than the MA filter for some asymmetric distributions, making the LW filter a competition to the MA filter. Besides, there is no constraints in terms of the filter length; however, longer length induces greater alarm delay, and how to analytically achieve the trade-off between the detection errors and the alarm delay is out of the scope of this paper. Similar to the order selection for delaytimers in Adnan et al. (2011), a simple way is to choose the least order of the MA filter that fits the detection error requirements based on industrial needs. In the simulation, the detection error requirements are set as FAR ≤ 10% and MAR ≤ 10% for both Gaussian and industrial data sets. With the help of the receiver operating characteristic (ROC) curve (details can be referred to Adnan et al., 2011), the order of 𝑁 = 5 has been chosen and considered to verify the effectiveness of the proposed method. Also, the sensitivity measures with different orders (𝑁 = 4, 𝑁 = 5 and 𝑁 = 6) are evaluated and compared.

4.1.2. Sensitivity over noise In this subsection, the sensitivity over noise is studied, where both the proposed analysis and the theoretical analysis are used. In the simulations, the ranges of the mean and standard deviation of the noise are chosen as 𝑚 ∈ [−0.2𝑚, ̄ 0.2𝑚] ̄ and 𝜎 ∈ [0.8𝜎, ̄ 1.2𝜎], ̄ where 𝑚̄ and 𝜎̄ represent the mean and standard deviation of the data set respectively. Thus, the noise mean and standard deviation are within 20% of the counterparts of data. The step sizes for the mean and standard deviation are set as d𝑚 = d𝜎 = 0.1 to calculate the sensitivity values over the KLD. Considering Gaussian noise 𝜖𝑛 ∼  (𝑚, 𝜎 2 ), the filtered data is also Gaussian distributed with 𝑦𝜖,𝑛 ∼  (𝑚 + 𝑚𝑛 , 𝜎 2 ‖𝜽‖ + 𝜎𝑛2 ‖𝜽‖). With the help of the error function, the partial differentiations of the FAR over the noise mean and variance in (22) are obtained as

4.1. Simulation results on sensitivity measures In order to check the effectiveness of the sensitivity analysis in Section 3, Gaussian distributions are considered firstly. The data sets under the normal and abnormal modes are generated from Gaussian distributions, where 𝑥𝑛 ∼  (6, 9) and 𝑥𝑎𝑏 ∼  (10, 16). The size of each set is 1000, i.e., 𝑙𝑛 = 𝑙𝑎𝑏 = 1000.

𝜕FAR = √ 𝜕𝑚

4.1.1. Sensitivity over trip point In Fig. 1, both proposed and theoretical sensitivity results over the trip point using the MA filter are shown. The proposed results are obtained from (16) and (17), where the distribution knowledge is unknown. For the theoretical curves, the PDFs of the data are assumed 2 ). to be known and denoted with 𝑥𝑛 ∼  (𝑚𝑛 , 𝜎𝑛2 ) and 𝑥𝑎𝑏 ∼  (𝑚𝑎𝑏 , 𝜎𝑎𝑏

−

1

e

(𝑦𝑡𝑝 −𝑚𝑛 −𝑚)2 2‖𝜽‖2 (𝜎𝑛2 +𝜎 2 )

,

2𝜋(𝜎𝑛2 + 𝜎 2 )‖𝜽‖ (𝑦 −𝑚 −𝑚)2

(𝑦𝑡𝑝 − 𝑚𝑛 − 𝑚) − 𝑡𝑝 2 𝑛2 2 𝜕FAR 𝜎 = e 2‖𝜽‖ (𝜎𝑛 +𝜎 ) . √ 2 2 𝜕𝜎 (𝜎𝑛 + 𝜎 ) 2𝜋(𝜎 2 + 𝜎 2 ) 𝑛

34

Y. Xiong et al.

Control Engineering Practice 70 (2018) 29–39

(a).

(a).

(b).

(b).

Fig. 1. Sensitivity over the trip point for MA and LW filters: (a) FAR; (b) MAR.

Fig. 2. Sensitivity over the trip point for MA and LW filters with different orders: (a) MA; (b) LW.

Besides, according to (Do & Vetterli, 2002), the closed form expression of the KLD with the filtered Gaussian data are 𝐷𝐾𝐿 (𝑓𝑌 ,𝑛 (𝑦) ∥ 𝑓𝑌 𝜖,𝑛 (𝑦)) =

1 ln 2

𝜎𝑛2 + 𝜎 2 𝜎𝑛2

+

‖𝜽‖2 𝜎𝑛2 + 𝑚2 2‖𝜽‖2 (𝜎𝑛2 + 𝜎 2 )

the theoretical one. Specifically, peaks of both methods appear synchronously. The differences in the magnitudes of the peaks are caused by the limited size of the data set. It has been tested that the differences decrease with larger data size. Secondly, the effectiveness of the proposed method is verified by the average sensitivity curves shown in Figs. 5 and 6. It is obvious that for both the MA and the LW filters, the proposed sensitivity results with respect to the FAR and MAR can tightly match the ones derived with the PDF information. The differences are also caused by the limited size of the data set. Additionally, from Figs. 5 and 6, it can be concluded that the average sensitivity values of the LW filter over the KLD are smaller than those of the MA filter in most of the entire trip point range in terms of the sensitivity of the FAR over the KLD. As for the sensitivity of the MAR, the magnitude of the LW filter is greater than that of the MA filter for the trip point ranging from 2 to 7, and it becomes smaller than that of the MA filter between 7 to 10. This reveals that the LW filter is less sensitive in terms of the FAR but more sensitive in the MAR case. Meanwhile, the sensitivity trends of the MA and LW filters with different orders are shown Fig. 7 by setting 𝑐1 = 𝑐2 = 1 in (36). It can be seen from Fig. 7(a) that higher order MA filters gets less sensitive to noise in the collected data, while the LW filter becomes more sensitive to noise according to Fig. 7(b). One explanation for this phenomenon

.

Then it is straightforward to get the partial differentiations of the KLD over the noise mean and standard deviation as 𝜕𝐷𝑘𝑙 𝑚 = , 𝜕𝑚 ‖𝜽‖2 (𝜎𝑛2 + 𝜎 2 ) 𝜎(‖𝜽‖2 𝜎𝑛2 + 𝑚2 ) 𝜕𝐷𝑘𝑙 𝜎 = − . 2 2 𝜕𝜎 𝜎𝑛 + 𝜎 ‖𝜽‖2 (𝜎𝑛2 + 𝜎 2 )2 With the same filter and trip point, the results with respect to the MAR in (27) can be calculated by similar steps. For the MA filter at the trip point 𝑦MA 𝑡𝑝 = 7.89, the sensitivities of the FAR and MAR over the KLD are shown in Figs. 3 and 4 respectively. The proposed results are shown in Fig. 3(a) and 4(a); while he theoretical results are depicted in Figs. 3(b) and 4(b). Additionally, when the trip point ranges from 2 to 10 with the step size 0.01, the average sensitivity of the FAR over the KLD calculated with (32) using both the MA and the LW filters are plotted in Fig. 5. Similarly, the average sensitivity of the MAR over the KLD calculated with (33) are depicted in Fig. 6. Firstly, from Figs. 3 and 4, it can be observed that the sensitivities derived with the proposed method can precisely catch the trend of 35

Y. Xiong et al.

Control Engineering Practice 70 (2018) 29–39

(a).

(a).

(b).

(b).

Fig. 3. Sensitivity of the FAR over the KLD for the MA filter: (a) proposed results without PDF information; (b) theoretical results with PDF information.

Fig. 4. Sensitivity of the MAR over the KLD for the MA filter: (a) proposed results without PDF information; (b) theoretical results with PDF information.

is that the noise reduction effect of the MA filters gets more powerful when the order increases. 4.2. Sensitivity-based filter design under Gaussian distribution In this subsection, the sensitivity-based design procedure introduced in Section 3 is applied for the same Gaussian data. In the simulation, the step size in each dimension of the filter coefficient is 0.1, i.e., 𝑠 = 0.1. Other parameters (d𝑚, d𝜎 and the ranges of the noise mean and standard deviation) are set to be the same with those in the previous subsection. The sensitivity measures of the proposed filter, the MA filter and the LW filter are listed in Table 2. Also, the sensitivity over the trip point is treated equally to that over the KLD, i.e., 𝜇3 = 𝜇4 . Three constraint levels are simulated, which are denoted with different sensitivity constraint levels: 𝜇3 = 𝜇4 = 5%, 𝜇3 = 𝜇4 = 10% and 𝜇3 = 𝜇4 = 15%. The corresponding upper-bounds of the sensitivities in each constraint (35) and (36) are listed in the first row in Table 2, denoted with 𝑈3 and 𝑈4 and calculated with (37) and (38), respectively. Besides, taking the case when 𝜇3 = 𝜇4 = 15% as an example, the sensitivities of the proposed filter, the MA and the LW filters over the trip point are shown in Fig. 8; and the sensitivity results over the KLD using the three filters are plotted in Fig. 9. Firstly, from Table 2, the trade-off between detection errors and sensitivity values can be observed. On one hand, the MA filter outperforms the other two linear filters in terms of the detection error, but its

Fig. 5. The average sensitivity of the FAR over the KLD for the MA and LW filters.

sensitivity is larger than others. On the other hand, when the constraint level changes from 𝜇3 = 𝜇4 = 5% to 𝜇3 = 𝜇4 = 15%, the summation of detection errors decreases, while the sensitivity of the proposed filter becomes larger. Also, it has been confirmed by simulations that when 36

Y. Xiong et al.

Control Engineering Practice 70 (2018) 29–39

Table 2 Sensitivity values with different constraint levels with Gaussian data. Constraint Levels [𝜇3 , 𝜇4 ] Constraints [𝑈3 , 𝑈4 ]

[5%, 5%] [4.43, 0.69]

[10%, 10%] [8.27, 1.28]

[15%, 15%] [12.34, 1.88]

MA

LW

Proposed 𝜽 𝑦𝑡𝑝

[ − 0.7, 0.1, 0.4, 0.8, 0.4 ] 8.78

[ − 0.1, 0.3, 0.4, 0.5, −0.1 ] 8.21

[ − 0.1, 0.2, 0.3, 0.4, 0.2 ] 7.97

[ 0.20, 0.20, 0.20 0.20, 0.20 ] 7.89

[ 0.04, 0.08, 0.12, 0.16, 0.20 ] 4.76

𝐹 𝐴𝑅 + 𝑀𝐴𝑅 𝑦 𝑦𝑡𝑝 |𝑆𝐹𝑡𝑝𝐴𝑅 | + |𝑆𝑀𝐴𝑅 | |𝑆̄ 𝐾𝐿𝐷 | + |𝑆̄ 𝐾𝐿𝐷 |

0.67 4.36 0.64

0.42 8.21 1.13

0.32 12.23 1.56

0.19 17.96 1.73

0.24 15.97 2.27

𝐹 𝐴𝑅

𝑀𝐴𝑅

(a). Fig. 6. The average sensitivity of the MAR over the KLD for the MA filter.

the constraint level is 𝜇3 = 𝜇4 = 20%, i.e., 𝑈3 = 16.46 and 𝑈4 = 2.50, the MA filter becomes the sensitivity-based optimal filter. Secondly, how the sensitivity constraints work in the filter design is studied when the sensitivity level is 𝜇3 = 𝜇4 = 15%, i.e., 𝑈3 = 12.34, 𝑈4 = 1.88. The constraint on the trip point is considered first in Fig. 8, where the summation of sensitivity magnitudes over the FAR and MAR is drawn with different trip points. The constraint is labeled with dotted green line for reference. It can be seen that because of the constraint, the feasible trip point ranges of the proposed filter, the LW filter and the MA filter are truncated to be [2, 4.02], [2, 5.93] and [2, 7.91] respectively. Similarly, the second sensitivity constraint about the sensitivity over the KLD in (34) is shown in Fig. 9. It can be seen that the proposed filter and the MA filter can satisfy the constraint with the feasible ranges for the trip point denoted as [3.02, 9.98] and [7.69, 8.81], respectively; while the LW filter cannot meet the second sensitivity constraint, thus it is infeasible. Then, combining the two constraints, there is no intersections for neither of the LW and MA filters. Finally, it can be concluded that the proposed filter has a better balance on the detection error and sensitivity than the other two linear filters.

(b).

Fig. 7. Sensitivity over the KLD for MA and LW filters with different orders: (a) MA; (b) LW.

4.3. Sensitivity-based filter design with industrial data

simulated. The corresponding upper bounds are listed in the first row of Table 3. By taking the third case as an example, the effects of the constraint with respect to the sensitivities over the trip point on the proposed filter, the MA and the LW filters are shown in Fig. 11, and those with respect to the sensitivity over the KLD are shown in Fig. 12. Similar to the observations in the Gaussian case, the trade-offs between the performance error and the sensitivity measure also exist with industrial data. Specifically, the sensitivity over the trip point and over the KLD both increase with a lower sensitivity constraint, while the weighted sum of the detection errors experiences a downward trend. Still, the MA filter is more sensitive to both the trip point and the KLD changes than the proposed filter, even though the MA filter has lower detection errors.

In this section, industrial data collected from a Chinese power plant in 2013 is used for the sensitivity analysis and design. The data plot is shown in Fig. 10, where the data points are sampled every 0.1 seconds and the data sizes are 𝑙𝑛 = 𝑙𝑎𝑏 = 1000. Same as the Gaussian case, in the simulation, the step size of the trip point is set to be 0.01, and the step size of the linear filter is 0.1. In the calculation of the sensitivity over the KLD, the range of the noise mean under the normal and abnormal modes are 𝑚𝑛 ∈ [−0.45, 0.45] and 𝑚𝑎𝑏 ∈ [−0.49, 0.49] with the step size 0.01; and the ranges of the noise standard deviation are 𝜎𝑛 ∈ [0.52, 0.78] and 𝜎𝑎𝑏 ∈ [0.53, 0.79] with a step size 0.01. By following the sensitivity-based filter design procedure in Section III, the sensitivity values with different constraint levels 𝜇3 = 𝜇4 = 5%, 𝜇3 = 𝜇4 = 10% and 𝜇3 = 𝜇4 = 20% are 37

Y. Xiong et al.

Control Engineering Practice 70 (2018) 29–39 Table 3 Sensitivity values with different constraints levels with the industrial data. Constraint Levels [𝜇3 , 𝜇4 ] Constraint [𝑈3 , 𝑈4 ]

[5%, 5%] [25.60, 5.33]

[10%, 10%] [34.70, 10.25]

[20%, 20%] [52.91, 22.34]

MA

LW

Proposed 𝜽 𝑦𝑡𝑝

[−1, 0.3, 1.7] 1.0, −1.6 ] 23.59

[ 0.4, 0.5, 0.7 0.5, −0.1 ] 23.60

[ 0.1, 0.2, 0.4, 1, −0.7 ] 23.61

[ 0.20, 0.20, 0.20, 0.20, 0.20 ] 23.61

[ 0.04, 0.08, 0.12 0.16, 0.20 ] 14.23

𝐹 𝐴𝑅 + 𝑀𝐴𝑅 𝑦 𝑦𝑡𝑝 |𝑆𝐹𝑡𝑝𝐴𝑅 | + |𝑆𝑀𝐴𝑅 | |𝑆̄ 𝐾𝐿𝐷 | + |𝑆̄ 𝐾𝐿𝐷 |

0.30 23.81 4.79

0.26 33.46 8.68

0.23 51.14 15.61

0.18 60.72 28.67

0.41 32.86 2.49

𝐹 𝐴𝑅

𝑀𝐴𝑅

Fig. 10. Industrial data.

Fig. 8. Three filters comparison in terms of the average sensitivity over the trip point with Gaussian data.

Fig. 11. Three filters comparison in terms of the average sensitivity over the trip point with industrial data.

Fig. 9. Three filters comparison in terms of the average sensitivity over the KLD with Gaussian data.

as the inputs; while the FAR and MAR are modeled as the outputs. By relating the definition of the elasticity, sensitivity measures are defined to quantify the performance change caused by variations in the inputs. With the help of the Gaussian Kernel based method, analytical results of sensitivity values are derived, which do not require knowledge on data distributions. Based on these sensitivity measures, a linear filter design problem is formulated in a constrained minimization structure, where sensitivity requirements are incorporated as constraints and the summation of the FAR and MAR is the objective function. The grid search is used to find the optimal sensitivity-based linear filter. Simulation results with the proposed filter and the MA and LW filters are conducted for the Gaussian and industrial data respectively. By comparing proposed results derived without any distribution information and the theoretical ones obtained with the known PDFs, the proposed

Also, from Fig. 11, it can be noticed that neither the MA filter nor the LW filter can meet the sensitivity requirement with respect to the trip point. Thus, the proposed filter still outperforms the other two filters in terms of the sensitivity taking into account both detection errors and sensitivity. 5. Conclusion In this paper, both sensitivity analysis and sensitivity-based design are studied for linear alarm filters. To quantify the sensitivity, a local sensitivity measure, namely, elasticity, is adopted and customized. In the sensitivity model, uncertainties in the trip point and in the KLD between the PDF of the noise-free data and that of the noisy data are modeled 38

Y. Xiong et al.

Control Engineering Practice 70 (2018) 29–39 Ahmed, K., Izadi, I., Chen, T., Joe, D., & Burton, T. (2013). Similarity analysis of industrial alarm flood data. IEEE Transactions on Automation Science and Engineering , 10(2), 452– 457. Borgonovo, E., & Apostolakis, G. E. (2001). A new importance measure for risk-informed decision making. Reliability Engineering & System Safety, 72(2), 193–212. Borgonovo, E., & Plischke, E. (2016). Sensitivity analysis: A review of recent advances. European Journal of Operational Research, 248(3), 869–887. Cheng, Y., Izadi, I., & Chen, T. (2011). On optimal alarm filter design. In Proc. of IEEE international symposium on advanced control of industrial processes, ADCONIP (pp. 139– 145). Cheng, Y., Izadi, I., & Chen, T. (2013). Optimal alarm signal processing: Filter design and performance analysis. IEEE Transactions on Automation Science and Engineering , 10(2), 446–451. Do, M. N., & Vetterli, M. (2002). Wavelet-based texture retrieval using generalized Gaussian density and Kullback-Leibler distance. IEEE Transactions on Image Processing , 11(2), 146–158. Gunawan, S., & Azarm, S. (2004). Non-gradient based parameter sensitivity estimation for single objective robust design optimization. Journal of Mechanical Design, 126(3), 395–402. Hall, P., Sheather, S. J., Jones, M., & Marron, J. S. (1991). On optimal data-based bandwidth selection in kernel density estimation. Biometrika, 78(2), 263–269. Izadi, I., Shah, S. L., Shook, D. S., Kondaveeti, S. R., & Chen, T. (2009). A framework for optimal design of alarm systems. In Proceedings of 7th IFAC symposium on fault detection, supervision and safety of technical processes (pp. 651–656). Karnavas, W. J., Sanchez, P. J., & Bahill, A. T. (1993). Sensitivity analyses of continuous and discrete systems in the time and frequency domains. IEEE Transactions on Systems, Man and Cybernetics, 23(2), 488–501. Kullback, S., & Leibler, R. A. (1951). On information and sufficiency. The Annals of Mathematical Statistics, 22(1), 79–86. Lai, S., & Chen, T. (2015). A method for pattern mining in multiple alarm flood sequences. Chemical Engineering Research and Design. Naghoosi, E., Izadi, I., & Chen, T. (2011). Estimation of alarm chattering. Journal of Process Control, 21(9), 1243–1249. Parzen, E. (1962). On estimation of a probability density function and mode. The Annals of Mathematical Statistics, 33(3), 1065–1076. Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., et al. (2008). Global sensitivity analysis: The primer. John Wiley & Sons. Tsay, R. S. (2005). Analysis of financial time series. New York: John Wiley & Sons. Walck, C. (2007). Handbook on statistical distributions for experimentalists, International report, SUF-PFY/96-01. University of Stockholm. Xiong, Y., Jing, Y., & Chen, T. (2017). Performance analysis on linear alarm filters. In Proc. American control conference, ACC (pp. 4424–4429). Yang, F., Shah, S. L., Xiao, D., & Chen, T. (2012). Improved correlation analysis and visualization of industrial alarm data. ISA Transactions, 51(4), 499–506. Zames, G., & Francis, B. (1983). Feedback, minimax sensitivity, and optimal robustness. IEEE Transactions on Automatic Control, 28(5), 585–601. Zang, H., Yang, F., & Huang, D. (2015). Design and analysis of improved alarm delaytimers. IFAC-PapersOnLine, 48(8), 669–674.

Fig. 12. Three filters comparison in terms of the average sensitivity over the KLD with industrial data.

performance sensitivity analysis is verified to be effective. Furthermore, trade-offs between sensitivities and performance errors are observed. The proposed linear filter is shown to outperform the MA and LW filters when the sensitivity is considered. Acknowledgments The authors would like to thank the Associate Editor and reviewers for their constructive comments with have improved this work, they are also grateful to the Natural Sciences and Engineering Research Council (NSERC) of Canada and the Chinese Scholarship Council (CSC) for the financial support. References Abdelaziz, T. H. (2015). Robust pole assignment using velocity–acceleration feedback for second-order dynamical systems with singular mass matrix. ISA Transactions, 57 , 71– 84. Adnan, N. A., Izadi, I., & Chen, T. (2011). On expected detection delays for alarm systems with deadbands and delay-timers. Journal of Process Control, 21(9), 1318–1331.

39