Real time algorithm based on time series data abstraction and hybrid bond graph model for diagnosis of switched system

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

Contents lists available at ScienceDirect

Engineering Applications of Artificial Intelligence journal homepage: www.elsevier.com/locate/engappai

Real time algorithm based on time series data abstraction and hybrid bond graph model for diagnosis of switched system

MARK

⁎

Mariam Taktak , Slim Triki, Anas Kamoun Research Unit on Renewable Energies and Electric Vehicles (RELEV), University of Sfax, Sfax Engineering School (ENIS), Tunisia

A R T I C L E I N F O

A BS T RAC T

Keywords: Switched System Hybrid Bond Graph Qualitative Diagnosis Parametrized Temporal Causal Graph Piecewise Aggregate Approximation Page-Hinkley Test

In this paper, we propose a real time algorithm to realize a diagnosis of switched systems for abrupt parametric faults. This algorithm is based on interaction between a Qualitative Diagnosis (QD) and a monitoring component that performs a Qualitative Trend Analysis (QTA) of residual signals generated from Bond Graph (BG) elements called residual sinks. The QTA is applied in order to on-line detect change in the mean of residual signal based on combination of Piecewise Aggregate Approximation (PAA) with Page-Hinkley Test (PHT). The QD procedure is performed in two stages. In the first off-line stage, Symbolic Fault Signature Matrix (SFSM) is generated from a Parameterized Temporal Causal Graph (PTCG). The PTCG is valid for all system modes and deduced from a unified Hybrid Bond Graph (HBG) model by converting its elements into node and labeled edge. Each entry in the SFSM matrix gives the residual symbolic signature which corresponds to the lower-order signature predicted using the PTCG model by propagating initial deviation from the fault parameter in the label of edge to the residual node. In the second on-line stage, trend extraction by linear regression is triggered after change detection in order to estimate the lower order time-derivative symbol for each residual sinks. Subsequently, we propose a stepwise similarity measure for fault isolation task. The functioning of this approach is illustrated in simulating with a switched quarter-car active suspension system.

1. Introduction Physical systems with switching are termed hybrid dynamic systems (Kypuros, 2001). Moreover, switches are physical mechanisms that discontinuously redirect the flow of power in a system. Semiconductor switches, hydraulic valves and mechanical clutches are example of switched elements mechanisms. When switching occurs, the system may change its mode of operation and system with n switching states has 2n possible operating modes. Consequently, systems that employ physically switched elements are difficult to model. As the complexity of systems increases, Fault Detection and Isolation (FDI) become more important since it is an essential means to maintain system safety and reliability. Typically, fault diagnosis includes three tasks: fault detection, fault isolation and fault identification (Gertler, 1998; Chiang et al., 2001; Precup et al., 2015; Isermann, 2009; Ding, 2008; Korbicz et al., 2004). For fault detection task, residual-based approach has attracted much attention (Niu et al., 2015; Serdio et al., 2014a; Borutzky, 2009a; Staroswiecki and Comtet-Verga, 2001). It consists in the generation and evaluation of diagnostic signals called residuals. Depending upon the knowledge (model-based or datadriven) used and the nature of information processing (quantitative or

⁎

qualitative) various techniques can be applied to generate residuals and evaluate them. However, the design of a residual-based approach should be adapted when we deal with a switching behavior within system dynamic. In fact, due to change in the operating mode of the switched system, time series of measurement exhibits singularities features even when no faults are present in the system. Recently, a data-driven fault detection (FD) tool is presented in (Serdio et al., 2014a, 2014b) and extended in (Serdio et al., 2015) by integrating a fault isolation (FI) level based on gradient and quality information. In the FD framework, residuals are obtained through data-oriented system identification (SysId) models. The particularity of this framework is its ability to resolve the singularities, in measurements, by transforming the original time-amplitude space into a reduced dimensional product-amplitude space (also called “residual space”) over a different set of measurements. While the SysId models are given under several forms (including linear, nonlinear and fuzzy models), the proposed framework is only well suited for sensors and actuators FDI. In order to isolate parametric faults, as it is our purpose in this work, more a priori knowledge is needed to have properly structured residuals. Among several residual generation methods that have been developed, Analytical Redundancy Relation (ARR) provides the possi-

Corresponding author. E-mail addresses: [email protected] (M. Taktak), [email protected] (S. Triki), [email protected] (A. Kamoun).

http://dx.doi.org/10.1016/j.engappai.2016.12.009 Received 14 December 2015; Received in revised form 17 August 2016; Accepted 5 December 2016 0952-1976/ © 2016 Elsevier Ltd. All rights reserved.

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

M. Taktak et al.

expected behavior of the system, given by model, with its actual behavior. The first step of a FDI procedure consists in generating a set of residuals which reflect the discrepancy between the two behaviors.

bility to have such structured property. That is, ARR correspond to relations between system variables which can be explained from the physical model and defined under symbolic or numeric form (Staroswiecki and Comtet-Verga, 2001; Ding, 2008). Moreover, ARR can be easily generated from graphical models that are commonly used to describe complex systems such as causal graphs, bipartite graphs and bond graphs (Yang et al., 2014). Compared to other graphical model, the causal properties of the Bond Graph (BG) provide more detailed information that allows the FDI algorithm to isolate a fault at component level easily (Ould Bouamama et al., 2014; Borutzky, 2009b). BG is based on analyzing the flow of energy, which is the product of a flow variable and an effort variable, through the interconnections of components called ports (Karnopp et al., 2000). The flow of energy between ports is shown by a line, which is called bond. When two elements are interconnected, there are power interactions between them. The power can be in different forms: mechanical, electrical, hydraulic or thermal. Power is product of two variables called power variables; e. g. mechanical power is the product of force (effort) and velocity (flow), whereas hydraulic power is the product of pressure (effort) and flow rate (flow). In BG model, any physical system (or subsystem) can be described by elementary components which include source elements Se and Sf, dissipative element R, storage elements C and I, two junctions 0 and 1, and two transformers TF and GY (Karnopp et al., 2000; Kypuros, 2013). As BG modeling is based on the exchange of energy between components, the methodology was initially used to capture continuous time phenomena. Extensions that also cover hybrid system models have been reported in the literature (Dauphin-Tanguy and Rombaut, 1989; Strömberg et al., 1993; Mosterman and Biswas, 2000; Umarikar and Umanand, 2005). Modeling of hybrid systems using a BG is one of the topics of research and various models have been proposed (Borutzky, 2009a; Pulido and Alonso-Gonzlez, 2004; Mosterman and Biswas, 1999; Manders et al., 2000; Narasimhan, 2002). To include discrete transition and modeling switching phenomena, additional mechanisms are introduced into the continuous BG language. In this work we adopt the concept of Hybrid Bond Graph (HBG) proposed by (Mosterman and Biswas, 2000), which extends the ability of BG to model hybrid systems using controlled junctions. Each controlled junction has two possible states, ON or OFF, which corresponds to an active or an inactive junction. Therefore, a controlled 1-junction is used to inhibit flow and a controlled 0-junction is used to inhibit effort. The concept of controlled junctions is an intuitive representation for structural discontinuities because they show clearly where elements connect and disconnect and break the path of power flow (Rebecca et al., 2013). Those junction switching function are implemented as a finite state automaton control specification. The finite state automaton may have several states, and each state maps to either the OFF mode or the ON mode of the junction. Modes transition expressed only by external controller signals define controlled switching, and those expressed by internal variables crossing boundary values define autonomous switching. The overall mode of the system is determined by a parallel composition of modes of the individual switched junctions.

2.1. Quantitative approaches for FDI of switched systems In the quantitative FDI approaches, ARR methods are classically used for residual generation (Borutzky, 2009a; Low et al., 2010a). ARRs are static or dynamic constraints which link the time evolution of known (i.e. parameters, inputs and measurements) variables when the system operates according to its normal operation model. ARRs have to be sensitive to faults and insensitive to perturbations. ARRs are derived from the set of over-determined equations obtained from the structural system model. In the BG framework, this translates to first generating equations that correspond to the conservation laws at each 0− and 1− junction, and manipulating these equations till only known variables remains. A concept of Global Analytical Redundancy Relation (GARR) has been proposed to extend the ARR based fault for switched dynamic systems (Low et al., 2010b). Aiming at a set of ARRs that holds for all system modes, authors in (Low et al., 2008) study the computational causality effects controlled junctions have on adjacent parts of the BG under alternative causality assignments and propose a modification of the Sequential Causality Assignment Procedure (SCAP). The main goal of the so-called Sequential Causality Assignment Procedure for Hybrid System (SCAPH) was to assign the HBG a desirable causality that facilitates analysis and FDI designs by introducing a Boolean variable to model the eliminations of flow or effort variables due to the OFF states of the controlled junctions. Fig. 1 shows the flowchart of the SCAPH procedure and the summary of the algorithm to derive the Diagnostic Hybrid Bond Graph. In model-based FDI for switched dynamic systems, GARRs are HBG-based method for symbolic residual generation. However, as for continuous dynamic systems, this method requires high computational costs for equation derivation and can not be applied when unknown variables can not be eliminated because of the presence of algebraic loops and non-linear non-invertible constraints (Borutzky, 2009a, 2012; Samantaray et al., 2006). As solution to these problems, author in (Borutzky, 2009a, 2012) proposed additional BG element termed as “residual sinks” used for a numerical computation of the residuals without having to derive ARRs in symbolic form. In contrast to another proposed method (Samantaray et al., 2006) where sensor variables are represented as sub-graph coupled to a Diagnosis Bond Graph model with preferred derivative causality, residual bond graph sinks use integral causality as the preferred computational causality in BG model. Fig. 2 gives a representation of residual effort sink and residual flow sink denoted as rSe and rSf, respectively. For both, the difference (Δe and Δf) between the measured variable obtained from sensor in real process and the output of the normal behavior BG-model is fed into modulated sinks that deliver the residual signal res. As shown in Fig. 2, these outputs from the residual sinks, denoted res, become additional inputs to the BG-model that alter its normal behavior until the difference (Δe and Δf) vanishes.

2. Switched systems modeling and diagnosis 2.2. Qualitative approaches for FDI of switched systems HBG is a BG-based modeling approach which provides an effective tool not only for dynamic modeling but also for FDI of switching systems (Borutzky, 2012; Low et al., 2010a). BG has been proven useful for FDI for continuous systems (Borutzky, 2009a; Pulido and AlonsoGonzlez, 2004; Mosterman and Biswas, 1999; Manders et al., 2000; Samantaray et al., 2006). In these model-based approaches, modeling has been used for both qualitative and quantitative FDI. The general principle of these BG model-based FDI approaches is to compare the

In qualitative FDI approaches, methods such as possible conflicts (Pulido and Alonso-Gonzlez, 2004), and analysis of Temporal Causal Graphs (TCG) (Mosterman and Biswas, 1999) has developed for diagnosis of systems. These methods are based on the structural analysis of dynamic models. For example, authors in (Mosterman and Biswas, 1999) developed a diagnosis schema, which called TRANSCEND, based on the qualitative analysis of fault transient

52

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

M. Taktak et al.

Fig. 1. Flowchart of the SCAPH procedure and the generation of the Diagnosis Hybrid Bond Graph.

to Hybrid TRNASCEND where the diagnosis algorithms use a combination of qualitative and quantitative reasoning mechanisms and the fault detection is based on a hybrid observer. Subsequently in (Daigle, 2008), author propose a qualitative event-based approach to fault diagnosis of hybrid systems that extends the TRANSCEND and Hybrid TRANSCEND methodologies to incorporate discrete faults. Other attention to hybrid diagnosis concentrates on efficiently processing a

information in the observed variables for diagnosis of continuous systems. In TRANSCEND, the fault isolation scheme is implemented in three step process: (1) hypothesis generation after fault detection, (2) fault signature generation for all hypothesized faults, and (3) hypothesis refinement through progressive monitoring. Since hybrid systems cannot be described by a single continuous or discrete event model, author in (Narasimhan, 2002) extend the TRANSCEND system

Fig. 2. Representation of residual effort sink and residual flow sink (Borutzky, 2012).

53

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

M. Taktak et al.

Fig. 3. Temporal Causal Graph transformation from BG model (Daigle, 2008).

constructed for a particular mode because the change of mode requires the re-assignement of causality using classical SCAP procedure. But, in case of many locally acting switches, the combinatorial explosion quickly leads to an intractable problem.

set of TCGs in (Narasimhan et al., 2000). In all proposed frameworks, system uses a Temporal Causal Graph model in order to analysis fault transients. Then, the diagnosis is based on this analysis, where observed deviations from nominal behavior expressed in qualitative form are compared against qualitative predictions of faulty behavior, i.e., fault signature, to isolate faults. Note that TCG is a linear causal graph which represents the model structure by linking various nodes representing variables in the model (Mosterman and Biswas, 1999). The directed edges between nodes represent constitutive relations of elements in a BG model and additive constraints are represented at the nodes. In addition, the TCG can be traversed in both forward and backward directions from an observed or hypothesized fault (Mosterman and Biswas, 1999). The backward propagation is used to construct a list of fault candidates, whereas the forward propagation derives predictions for posteriori temporal evolution for hypothesized faults. In the continuous system model, the BG is with fixed topology and the causality can be automatically assigned to the bonds model using four types of causal constraints defined by SCAP procedure (Karnopp et al., 2000): (1) fixed causality of source elements, (2) constrained causality for 0- and 1- junctions, (3) preferred integral causality for energy storage elements, and (4) arbitrary causality for dissipative elements. Then, the TCG may be directly derived from a BG by including temporal constraints that are important for analysis of dynamic systems. Fig. 3 gives the TCG transformation from BG model as defined by (Daigle, 2008). However, due to mode-varying causality properties of HBG, the TCG of switched system must be derived for each possible system configuration or mode. Hence, the TCG model is

2.3. Proposed framework for abrupt FDI of switched systems In this paper we propose a combination of quantitative and qualitative FDI approach that use a unified HBG model to generate and evaluate residual signals. For the residual generation, “residual sinks” are coupled with the sensor from the real process and added to a faultless unified HBG-model in order to numerically evaluate the GARR. Outputs signals delivred by these residual sinks are subsequently feeding into QTA, as shown in Fig. 4, in order to initially transform time series data into a sequence of coefficient that represent a constant approximation of the original data. We will call these coefficients “local means” and the proposed on-line change detection algorithm uses these coefficients in a modified Page-Hinkely Test. When change is detected, estimated symbolic signature of the residual signals are evaluated against predicted symbolic signature generated by forward propagation of deviation change through a Parametrized TCG, deduced from the unified HBG model. In (Mosterman, 2001), author describes how a hybrid model can be made amenable to the diagnosis algorithms that were developed in (Mosterman and Biswas, 1999) by systematically generating one parametrized TCG. The latter parametrized TCG is not derived from HBG, as we propose in this work, but from an implicit model formulation that includes mode selection

54

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

M. Taktak et al.

Fig. 4. Architecture of the proposed QTA and QD for abrupt fault diagnosis of switched system.

is solved numerically from a HBG model with preferred integral causality assignment to have a consistent and coherence results with QD when the PTCG will be used for symbolic signature prediction. The QD procedure is performed in two stages. In the first off-line stage, Symbolic Fault Signature Matrix (SFSM) is generated from PTCG model. This matrix is defined by the product of the set of residual sinks R={Resj | j=1,2,…,J}, and the set of parametric faults F={Fi| i=1,2,..,I} where parametric fault magnitude is abstracted symbolically as ‘+’ or ‘−’ symbol, representing above or below nominal value. The matrix element in the i-th row and the j-th column has the symbolic signature which correspond to the lower-order signature predicted using the PTCG model by propagating initial deviation from the labeled edge Fi to the residual node Resj (see Section 2.2 and Section 3.1). Table 1 j shows a SFSM matrix where Si represent the predicted symbolic signature of Resj to parametric fault Fi. In the second on-line stage, trend extraction by linear regression is triggered after change detection in order to estimate the lower order time-derivative symbol for each residual sinks. These estimated values are also abstracted symbolically as ‘+’, ‘0′ and ‘−’ symbols. Subsequently, we propose a stepwise similarity measure for fault isolation task (see Fig. 4). In step-1, using results of on-line modified PHT and trend extraction algorithm, a local similarity match between estimated and predicted symbolic signature for the first detected residual sink is used to eliminate faults from SFSM that have predicted pattern distant from the estimated pattern. In step-2, a global similarity measure (by considering all estimated residual sinks signatures) is then applied to a small set of fault candidates in order to identify the most likely fault scenario. Because

Table 1 Symbolic Fault Signature Matrix. F/R

Res1

F1

S1

1

S1

F2

1 S2

… FI

…

ResJ

2

…

S1

2 S2

…

…

…

1

SI

… …

SI

Res2

2

J J

S2 …

J

SI

parameters to switch between equations. By using the implicit model formulation, causalities relations between variables vanishes and conflicting relation may arise. To resolve this problem, author in (Mosterman, 2001) propose the introduction of additional parameters to enforce mutually exclusion between different causal relations. This proposed solution may be well adapted for simple case and become ambiguous for more than ternary relations. In this work we focus on the unified HBG model to build a PTCG which require less computational resource. As the SCAPH deal with the change of causality in the HBG model, we propose to apply the SCAPH procedure on the unified HBG model. Then, the conversion of the causal unified HBG model to a PTCG will be valid for the overall system configuration where the causal paths are parameterized by the state of the controlled junction. However, storage components of the unified HBG model must be in integral causalities instead of in derivative as proposed in (Low et al., 2010b). In fact, it is crucial that residual signal

55

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

M. Taktak et al.

Table 2 Pattern Similarity Matrix. –

S0

0

+

S1

–

0

+

?

–

0

+

?

–

0

+

?

–

– 0 + ?

1 0.4 0.3 0.6

0.4 1 0.4 0.6

0.3 0.4 1 0.6

0.6 0.6 0.6 1

0.2 0 0 0.1

0 0 0 0

0 0 0.2 0.1

0.1 0 0.1 0.2

0 0 0 0

0 0 0 0

0 0 0 0

0 0 0 0

0

– 0 + ?

0.2 0 0 0.1

0 0 0 0

0 0 0.2 0.1

0.1 0 0.1 0.2

1 0.4 0 0.6

0.4 1 0.4 0.6

0 0.4 1 0.6

0.6 0.6 0.6 1

0.2 0 0 0.1

0 0 0 0

0 0 0.2 0.1

0.1 0 0.1 0.2

+

– 0 + ?

0 0 0 0

0 0 0 0

0 0 0 0

0 0 0 0

0.2 0 0 0.1

0 0 0 0

0 0 0.2 0.1

0.1 0 0.1 0.2

1 0.4 0.3 0.6

0.4 1 0.4 0.6

0.3 0.4 1 0.6

0.6 0.6 0.6 1 Fig. 5. Time dependency of faults a) abrupt, b) incipient and c) intermittent.

S0: symbol of zeroth order time-derivatire S1: symbol of first order time-derivative

3. Sygnal to symbol generation for diagnosis of switched system off-line symbolic signature prediction may leads to an uncertain value, ‘? ’, (see case study section for details) a fuzzy-match rather than crispmatch between symbolic signature is considered. Hence, local and global similarity indices are defined according to a user-defined fuzzy similarity-match matrix as proposed in (Dash et al., 2003) and (Maurya et al., 2007). For our case study application, we define the pattern similarity matrix shown in Table 2 in order to quantify local similarity between patterns in a fuzzy manner where the degree of match vary from 0 to 1. Each entry in Table 2 gives a local similarity index, LSI(Sym_Sig1,Sym_Sig2), between patterns Sym_Sig1 and Sym_Sig2 where each pattern is defined as a pair of symbols {S0,S1} (see Section 3.1). As we can see from Table 1, the signture of the i-th fault corresponds to the SFSM i-th row. Thus, the Fault Signature (FS) of parametric fault Fi is defined as the vector of residual sinks symbols corresponding to that fault: 1

2

J

In previous work (Triki and Kamoun, 2014), a trend extraction algorithm have been used to extract useful features of biological signal, generated from an ocular command device, in order to distinguish between ocular command movements and then assist a disabled person in navigation task. In this work we reveal the importance of the time series abstraction and the trend extraction approach; (i) to design a online abrupt fault detection algorithm based on a sliding window technique and a modified Page-Hinckley Test that enhance the robustness by evaluating the difference between what we call local and global means of the time series data, and (ii) to estimate the symbolic signature of 0th and 1st order time-derivative of the transient in the dynamic of residual sinks after the fault occurrence. 3.1. Related works for symbolic generation of fault transient after abrupt change

(1)

FS (Fi ) = [Si , Si , …, Si ]

Significant non-zero residual indicate faults which is defined as an undesired deviation of at least one characteristic property of a variable from a normal behavior. Therefore, the fault is a state that may lead to a malfunction or failure of the system. The time dependency of faults can be distinguished as shown in Fig. 5 (Isermann, 1984; Luo, 2006). In this work we focus on the abrupt change profile. For continuous systems, it is reported in (Mosterman and Biswas, 1999) that abrupt faults create transients in the dynamic plant behavior. This result is extended to the switched system in (Narasimhan et al., 2000), where fault occurrence or discrete mode switching cause instantaneous change in some system variables. These discontinuities are perceptible as transient behavior in the residuals which are defined as the difference between the actual plant output and the nominal plant output (r = y − yˆ ). Recall that in contrast to the GARR deduced from the Hybrid Bond Graph model, the nominal plant behavior in (Manders et al., 2000; Narasimhan, 2002; Daigle, 2008) is tracked by the mean of a hybrid observer. In TRANSCEND and HybridTRANSCNED, diagnosis framework is based on the qualitative abstraction of the transient observed in the residual signal. For this purpose, the residual after the occurrence of the fault is approximated by the kth order Taylor series expansion given by:

Consequently, the i-th fault can be written in the following rule form: 1

J

If (S1 → Si ) ∧ … ∧ (S J → Si ), thenFi

(2)

where

• • •

S j is the estimated symbolic signature of Resj j Si is the predicted symbolic signature of Resj to parametric fault Fi, and ∧ denotes a conjunction operator

The i-th fault isolation consists in determining global similarity index (GSI) of rule (2). Note that GSI is a measure of the degree of match between observed fault signature and that stored in SFSM. As stated in (Maurya et al., 2007), the GSI of the rule (2) whose premise is a conjunction of simple premises is calculated using the MIN operator because classification requires discrimination between different faults. Thus, the GSI of the rule (2) is calculated according to: 1

J

GSI (Fi ) = min(LSI (S1, Si ), …, LSI (S J , Si ))

(3)

where the j-th rule's simple premise fulfillment degree is obtained from j pattern similarity matrix and represent the LSI between S j and Si .

r (t ) = r (t f ) +

56

dr (t ) dt

(t − t f ) t=t f

1!

+

d 2r (t ) dt 2

(t − t f ) 2 t=t f

2!

+⋯+Rk (4)

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

M. Taktak et al.

Fig. 6. Transient evolution that may be observed in residual signal and their symbolic signature which include the 0th and 1st order time derivative.

respectively. As illustrated in Fig. 6, the set of possible lower order fault signatures, given by the combination of the three symbols {−,0, +}, is SymSig = {00, 0 + ,0 − ,+0, −0, + +,+ −,− −,−+}. Each symbolic signature in the SymSig set defines a pattern. Given that the HybridTRANSCEND framework in (Narasimhan et al., 2000) addresses only continuous fault isolation for switched system, authors in (Daigle, 2008) develop an event-based framework to addresses the discrete fault diagnosis and extend the symbol generation method by introducing three additional symbols, {N , X , Z}, in order to distinguish continuous from discrete faults. These three additional symbols are deduced from measured variable, that exhibit a discrete change, and not from residual as for defined symbols of continuous faults. The methodology proposed in (Daigle, 2008; Manders et al., 1999) assumes that discontinuities of the form {++} or {−−} will not occur because they imply positive feedback and then an unstable system. That is, the proposed symbol generation based on the Z-test will generate {0 + } for the case where a + is observed in both the magnitude and the slope. Whereas a − in the magnitude and a − in the slope will be reported as {0 − }. However, author in (Gelso et al., 2007) demonstrate that signatures of the form {++} and {−−} can exist for some continuous parametric faults in stable systems.

where t f is the time of abrupt fault detection and Rk is a remainder term based on kth order derivatives of r (t ). Thus, the fault signature, denoted as S k (r (t )), is the set of feature values consisting of the magnitude and the 1st through kth order derivative value computed from the residual at instant time of abrupt fault detection:

⎛ dr (t ) S k (r (t )) = ⎜⎜r (t f ), dt ⎝

, t=t f

d 2r (t ) dt 2

, …, t=t f

d k r (t ) dt k

⎞ ⎟ ⎟ t=t f ⎠

(5)

The Taylor series is a good approximation of the noise-free signal r (t ) when t is close to t f . Since measurements are always contaminate by the noise, the estimation of the higher order derivative of the transient observed in the residual signal becomes a challenging task. A common ways to remove noise are filtering methods such as a low-pass filter or a moving average. From a FD point of view, authors in (Serdio et al., 2014) demonstrate that the smoothing effect of filters on residuals signal when it is compared against adaptive threshold decrease the false alarms and improve the accuracy of the fault detection task. However, from a FDI point of view based on the shape of the transient in residual, post-processing filtering may tend to distort the residual signal during the transient and consequently smooth the effect of the abrupt fault or the discrete change mode. To overcome this problem, the lower order feature that only capture the zeroth and first order derivative can be estimated using the methodology explained in (Daigle, 2008; Manders et al., 1999) where the detection of abrupt fault is posed as a hypothesis test using Z-test and a set of three sliding windows with varying size. When a fault is detected, the qualitative magnitude (i. e. 0th order) is provided at the time of fault detection and a +(res. − ) symbol is taken if the residual was positive (res. negative). Then, the 1st order of the deviation in the residual, after fault detection, is calculated using also Z-test over several small windows at different position from the point of fault detection. Hence, the magnitude and derivative value are abstracted symbolically as ‘+’, ‘0′ and ‘− ’ symbols, representing above, at, and below nominal value

3.2. The proposed on-line algorithm for abrupt change detection and symbolic signature estimation We propose here an on-line algorithm based on a sliding window technique in order to (i) detect abrupt parametric fault from residual bond graph sinks and (ii) to sequentially generate the lower order of symbolic fault signature based on trend extraction and classification approach. Indeed, the proposed algorithm will independently estimate the zeroth and first order change such that a signature of {++} (res. {−−}) will not be reported as {0 + }(res. {0 − }). For the detection process, a combination of the Page-Hinkley Test (Mouss et al., 2004) with a time series data abstraction technique is proposed to enhance

57

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

M. Taktak et al.

Fig. 7. abrupt change detection based on PHT applied on two similar signals with different AWGN.

where a forgetting mechanism was added to weight data in order to favorite recent samples in the detection process and consequently reduces the detection time delay. The proposed mechanism was tested in a real-time algorithm for changes detection in depth of anesthesia signal. Recently, authors in (Shaker and Lughofer, 2014) extend the application of forgetting mechanism on both mT and MT and also integrate a dynamic forgetting factor which adapt the degree of forgetting independency of the change amplitude. Page-Hinckley algorithm show good performance for signal with high SNR and require optimal setting parameters to deal with low SNR. In the latter case a multiple trial-and-error test runs on various signals may leads to an optimal choice of these parameters. Given that the parameter of admissible change, δ , is highly dependent of the characteristic of the signal under study, it is quite sensitive to the SNR and require retuning stage even when signal shows slightly variation in SNR. The upper right and left plots on the Fig. 7 shows two synthetic signals, yσ1 and yσ2 , with Additive zero-mean White Gaussian Noise (AWGN) having standard deviation σ1 = 1.4 and σ2 = 1.5, respectively. Each signal has 1000 samples where a positive change in the mean occur between samples 400 and 700. Before and after the mean change, the signal without noise has a magnitude equal to 0.1 which mean that the noise is larger than the signal. The following is the equation used for generation the synthetic signal:

the performance of the PHT for signal with low SNR (Signal to Noise Ratio), i. e. for noisy signal. For the qualitative signature estimation, a trend extraction technique based on feature extraction about the data tendency and a classification process is proposed. Two main approaches, for abrupt change detection, can be found in the literature; statistical signal processing (Basseville and Nikiforov, 1993) and timefrequency analysis (Mallat and Hwang, 1992). As stated in (Mouss et al., 2004; Sebastiao et al., 2013) and (Shaker and Lughofer, 2014), PHT is one of a simple algorithm, easy to implement and commonly used to on-line detect changes in a signal processing. It represents a sequential adaptation for detecting abrupt change in the mean within Gaussian distributed observations. For example, to detect change in signal x (t ), PHT consists of running two tests simultaneously. For this purpose, two cumulative variables denotes as UT , for the increase direction, and L T , for the decreasing one, are defined as the cumulative difference between the observed values and their mean, T 1 xT = T ∑t =1 x (t ), till the current moment T : T ⎧ (x (t ) − xT − δ ) ⎪UT = ∑ t =1 ⎨ T ⎪ ⎩ L T = ∑t =1 (x (t ) − xT + δ )

(6)

The parameter δ , in Eq. (6), correspond to the admissible change that is allowed to avoid false detection due to noise. Letting mT = min(Ut , t = 1 … T ) and MT = max(L t , t = 1 … T ) the minimum and maximum value of cumulated variables in (6), respectively. The detection of change in the mean of the signal x (t ) is, then, determined using the following tests which are run in parallel:

⎧ if (PHU = UT − mT ) ≥ λ then increase change detected ⎪ ⎨ if (PHL = MT − L T ) ≥ λ then decrease change detected ⎪ otherwise no change detected ⎩

⎧ μ0 + σ⋅n (k )for0 ≤ k < 400 ⎪ yσ (k ) = ⎨ μ0 + σ⋅n (k ) + d for400 ≤ k ≤ 700 ⎪ μ + σ⋅n (k )for700 < k ≤ 1000 ⎩ 0 where:

• • • • • •

(7)

where λ is a predefined threshold that is chosen considering a tradeoff between delay time detection and admissible false alarm rate. Note that the mean x in Eq. (6) is incrementally updated according to the following equation:

x (t ) =

T−1 1 x (t − 1) + x (t ) T T

(9)

(8)

A novel concept for PHT was proposed in (Sebastiao et al., 2013),

58

μ0 : mean value of the discrete signal under normal behavior (fixed to 0.1), σ : standard deviation of the noise signal, n (∙): noise function that generates random normally distributed numbers between 0 and 1, d : magnitude deviation of the change, k : sample index, yσ (k ): sample value for a chosen standard deviation of the noise.

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

M. Taktak et al.

Fig. 8. Influence of the magnitude deviation and variance of AWGN in the PHT performance.

corresponds to a correctly detected change. In order to enhance the robustness of the PHT algorithm, we will propose some modification inspired from the methodology proposed in (Daigle, 2008). In fact, the mechanism of multiple windows with varying size used to estimate the mean of the residual in different position within a sliding window inspired us to combine the PHT with a technique of Piecewise Aggregate Approximation (PAA) (Keogh and Pazzani, 2000) widely used in time series data mining community. Basically, PAA is used to achieve dimensionality reduction (or abstraction) of time series data and build index in order to speed up the similarity search in large time series database. For illustration, let P = { p1, p2 , …, pn} be an initial time series of dimension n , the compressed time series data after dimensionality reduction using PAA is P = { p1, p2 , …, pN } where N is the dimensionality of the transformed space (1 ≤ N ≤ n ) and each element in P can be obtained by:

The results of PHT algorithm (as proposed in (Shaker and Lughofer, 2014) for predicting the change in the mean) applied on these signals are depicted on the bottom plot of the Fig. 7. In this example, the initialization parameters are fixed as follow; the forgetting factor α = 0.999 , the admissible change δ = 0.237 and the threshold to detect abrupt change λ = 0.645. These parameters are fixed after some trial applied on signal with AWGN having σ1 = 1.4 . We use the predefined parameters for both case and we remark that a slightly increase of the Gaussian white noise standard deviation leads to some few false alarm. Obviously, a re-tuning stage will lead to adapted values of parameters that reduce false alarm. However, depending on the magnitude deviation of the change and the SNR, false alarm may be either decrease or increase. For the simple case of the synthetic signal, we can see on the plots of Fig. 8 that when magnitude deviation increase and variance of the AWGN decrease, false alarms are considerably reduced. Otherwise, a false alarm of type I, know as a False Positive, arise especially when magnitude deviation is small and the variance of the AWGN is large. Note that, a False Positive (FP) corresponds to a wrongly detected change; a False Negative (FN) is the number of change that is not detected and a True Positive (TP)

pk =

N n

nk N

∑ i = n (k −1)+1 N

pi (10)

As we can remark from Eq. (10), the PAA consist of splitting the

59

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

M. Taktak et al.

The functional form Y j (k ), therefore, can be used to refer to the time series vector Y j . Indeed, Y j will be used to denote the jth onedimensional time series vector or data set. The proposed on-line algorithm for simultaneously detection and estimation of the symbolic signature of the transient caused by an abrupt fault can be described in the following steps: INPUT: a sliding window size N1 that contains n data points of the set Y , a fixed slide move distance N2 , a predefined dimensionality N for reduction of the time series in the sliding window (for convenience choose N a factor of n), a predefined threshold for change detection λ , a predefined threshold to distinguish change “with” or “without” discontinuity Λ and a predefined threshold to classify linear segment ths . OUTPUT: estimated symbolic signature for each time series data within the sliding window where the first change is detected. STEP1: normalize each time series within the sliding window to have a mean of zero and variance of one. That is,

Fig. 9. Illustrative example of the data reduction using PAA technique (n = 300 ,N = 10 ).

Yj ←

X2

a) Sequentially transform the set of time series in Y that's fall within the sliding window into a PAA form, denoted as Y . That is, Y = [Y 1, Y 2, …, Y m]T where Y j = {y1j , y2j , …, y Nj}. At each sequence, the ith element of Y j is calculated using Eq. (10):

yhj =

d 22 d32

d 23 d33

⋯

=

dkj

N

n

n

2: if Uhj < MinU j then MinU j ← Uhj 3: 4: endif 5: if Uhj > MaxU j then MaxU j ← Uhj 6: 7: endif 8: PHPaaLhj ← MaxU j − Uhj 9: ΔPHPaaLhj ← PHPaaLhj − PHPaaLhj−1 10: if ΔPHPaaLhj < 0 then MaxUpVariation ← PHPaaLhj−1 11: 12: IndexMaxUpVariation ← h − 1

(11)

13: ΔUp ← ΔPHPaaLhj MaxU j ← 0 14: Uhj ← 0 15: 16: endif 17: PHPaaUhj ← Uhj − MinU j

Note that X is an ordered sequence of n vectors each having dimension m such that m = r + d , while Y is an unordered collection of m one-dimensional time series vectors each of length n . Since each value dkj is associated with a particular sample time index k , Y j can be rewritten in the equivalent functional form:

Yj

nh

∑zN= n (h −1)+1 dzj where yhj is the local mean of the jth measure-

1: Uhj ← Uhj−1 ∪ {yhj − GM (Y j )}

d 2m d3m

⋯ X3 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ Xn dn1 dn2 dn3 ⋯ dnm

N n

ment for the interval [ N (k − 1) + 1, N k ]. b) collect, simultaneously: (i) the cumulative difference between local and global mean stored in U j , (ii) the minimum and maximum of the set U j , (iii) the difference (Uhj − MinU j ) stored in PHPaaUhj and the difference (MaxU j − U j ) stored in PHPaaLhj and (iv) the increment of PHPaaUhj and PHPaaLhj stored in ΔPHPaaUhj and ΔPHPaaLhj , respectively. That is,

Y1 Y 2 Y 3 ⋯ Y m d11 d12 d13 ⋯ d1m d 21 d31

for j = {1, …, m}

STEP2: for each measurement, determine the vector n 1 GM (Y ) = [GM (Y1), GM (Y 2 ), …, GM (Y m )]T where GM (Y j ) = n ∑i =1 di j is the global mean of the jth measurement within the sliding window. STEP3: initialize the following sets: U0j ← 0 , MinU j ← 0 , MaxU j ← 0 , ΔPHPaaU0j ← 0 , PHPaaU0j ← 0 , PHPaaL 0j ← 0 , ΔPHPaaL 0j ← 0 for j = {1, 2, …, m}, and initialize h ← 1 STEP 4: do the following subsets while h < N1:

initial time series data into N segments (or interval) and compute the average (mean) value of each segment. Fig. 9 illustrates the abstraction of time series using PAA. The proposed combination of PAA with PHT, to enhance the robustness for abrupt change detection of the mean in a signal, is implemented in three steps. The first step consists of collecting data of length n within a sliding window and calculates their mean value. This mean value is called a global mean. In second step, data are subdivided into N equi-spaced segments and the mean value of each segment is obtained by Eq. (10). The vector of N means values is called the vector of the local mean. In the third step, the PHT is applied through segments in the sliding window by substituting x (t ) and xT in Eq. (6) by the corresponding local and global mean, respectively. More formally, a time series data set X , in our application, is an ordered sequence {X1, X2, …, Xn} where Xk is a vector of samples data collected at the discrete sampling time index k such that the samples in Xk were recorded previous to those in Xk+1 for all 1 ≤ k ≤ n . Each vector Xk = {dk1, dk2, …, dkm} is composed of m measurements collected at the discrete sampling time index k . A transposition of the time series data set X results in a corresponding data set Y which is a collection of vectors {Y1, Y 2, …, Y m} where Y j contains the time series {d1j , d 2j, …, dnj} for all 1 ≤ j ≤ m . Here, the set Y is divided into two mutually exclusive set {Yr , Yd} such that Y = Yr ∪ Yd . Although discrete change is out of the scope of this work, the set Yd is defined, here, for a future extension of our work. Therefore, Yr ⊂ Y is the set of the r residual signal generated from the “residual sinks” coupled with the sensor of the real switched system and Yd ⊂ Y represents the set of d measured variables that exhibit the discrete change behavior in the switched system. Given this definition, a time series data set has the matrix form:

X1

Y j − mean (Y j ) standard deviation (Y j )

18: ΔPHPaaUhj ← PHPaaUhj − PHPaaUhj−1

(12)

60

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

M. Taktak et al.

Fig. 10. Rules of slope classification for the first order change estimation.

Fig. 11. Results of PHT with PAA applied on two synthetic signals with different AWGN standard deviation.

61

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

M. Taktak et al.

Fig. 12. Scheme of the quarter-car active suspension system.

Fig. 13. Hybrid Bond Graph model of the quarter-car system. 19: if ΔPHPaaUhj < 0 then

a) check whether the MaxUpVariation is greater than the predefined threshold λ :

20: MaxDownVariation ← PHPaaUhj−1 21: IndexMaxDownVariation ← h − 1 22: ΔDown ← ΔPHPaaUhj 23: 24:

1. if the above test is TRUE, then a positive change in the mean occurs in the jth measurement at sample which belong to the segment n n sF = [ N (h − 1), N h]. This change is qualified as “without discontinu-

MinU j ← 0 Uhj ← 0

25: endif

62

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

M. Taktak et al.

Fig. 14. Coupling of measured variable from real system and a unified HBG with nominal parameters for the quarter-car active suspension system.

tree (see Fig. 10) using the slope quantitative parameter and a predefined threshold ths . ths determines whether the slope is increasing {+}, decreasing {−} or steady {0}. That is, the qualitative signature is automatically generated as follow:

Table 3 Nominal physical parameters of the quarter-car system. Parameter

Value

Unit

Ms Mu Ks Kt Cs A1 A2 V1 V2 Ps Ctm ρ β Cd ω

290 59 16800 190000 1000 3.8 × 10−4 3.5 × 10−4 1.2 × 10−4 1.1 × 10−4 1.5 × 107 0.3 × 10–12 900 1.2 × 109 0.6 1.5 × 10−3

Kg Kg N/m N/m N s/m m2 m2 m3 m3 Pa m3 s−1/Pa Kg m−3 Pa – m

1: if Y j ⊂ {Yr} then 2:

S (Y j ) ← (0, +) 3: 4: else S 0 (Y j ) ← (+) 5: 6: if slope j ≤ ths then S1 (Y j ) ← (0) 7: elseif sign (slope j ) > 0 then S1 (Y j ) ← (+) 8: else S1 (Y j ) ← (−) 9: endif 10 : endif 11 : for k ← 1 to r do 12 : ifk ≠ j do S 0 (Y k ) ← (0) 13 : 14 : if slope k ≤ ths then S1 (Y k ) ← (0) 15 : elseif sign (slope k ) > 0 then S1 (Y k ) ← (+) 16 : else S1 (Y k ) ← (−) 17 : endif 18: endif 19 : endfor 20 : returnS (Y ) 21 : else if Y j ⊂ {Yd} 22 : // the class of discrete fault is reported to future work 23 : end if

ity” if the incremental value ΔUp is less than the predefined threshold Λ. Hence, a symbolic signature of 0th and 1st order is immediately deduced for the jth measurement as: S (Y j ) ← (0, +). Otherwise, the change is qualified as “with discontinuity” (or abrupt) and only the qualitative signature for the 0th order of the jth measurement is deduced as: S 0 (Y j ) ← (+). In order to estimate the qualitative signature of the 1st order derivative, a linear regression is computed from the original sampled data for an interval I of L data point collected just after the segment sF , n n I = [ N h , N h + L ]. Thus, the least square optimization of the linear slope leads to: slope j =

n h+L ⎛ N⎞ ∑ N n ⎜ i − ⎟ (d i j − d j ) i= h ⎝ n⎠ N n h+L ⎛ 2 N⎞ ∑ N n ⎜i − ⎟ i= h ⎝ n⎠ N

if ΔPHPaaLhj < Λ

where d j represent the 2. elsewhere, no change is detected and increment h . That is, h ← h + 1

average value of d j in the interval I . Hence, the qualitative value of the 1st order derivative time is obtained thanks to a simple decision

a) Check whether the MaxDownVariation is greater than the predefined 63

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

M. Taktak et al.

STEP 5: slide the window a fixed distance N2 and GO TO STEP 3. The robustness of the proposed PHT with PAA algorithm is illustrated through the simple example of synthetic signal generated from Eq. (9). Let us consider the worst case for the PHT algorithm where the false alarm of type I arise, i. e. the case where magnitude deviation is small and the standard deviation of the additive white Gaussian noise is high. So, for a small magnitude deviation of 0.4 , two

threshold λ : 1. if the above test is TRUE, then do same thing as in the first subsubstep c) except for line 3 and 5 where the sign (+ ) must be replaced by sign (− ) 2. elsewhere, no change is detected and increment h . That is, h ← h + 1

Fig. 15. PTCG derived from the coupled HBG model of the quarter-car system.

Fig. 16. road profile – valve displacement and the discrete state [α1, α2] of the switched system.

64

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

M. Taktak et al.

Fig. 17. Adjacency matrix of the quarter-car model operating at mode 1.

Fig. 18. Adjacency matrix of the quarter-car model operating at mode 2.

original signal as well as the interval of the fifteenth segment, [700, 750], include a maximum negative change in the mean of the original signal. The correct change detection point belong to the identified segment in both cases, therefore detection robustness is improved compared to the PHT without PAA. In addition, we can remark that detection accuracy depends to the segment size and consequently to the chosen dimension reduction value. However, in this work, we look for a simple and robust way to estimate the qualitative signature during the transient dynamic in the residual signal when an abrupt parametric fault occurs. Hence, rather than use PHT for a one to one sample test which may lead to false alarm especially when a slightly variation of the SNR happen, we simply use PHT through segments where each segment is a constant line that aggregate more than one sample. The proposed approach to estimate the qualitative signature is illustrated, in subsequent section, through an example of quarter-car active suspension switched system.

signals with additive white Gaussian noise of standard deviation σ1 = 1.5 and σ2 = 1.8 are used in this illustration. Note that, for visibility purpose, the proposed algorithm is applied off-line with dimensionality reduction N = 20 . Fig. 11 shows on the upper plot graph, the normalized noisy signal (where the zero mean represent the global mean) and the signal-reduced representation under piecewise aggregate approximation. In this graph, original signal is divided into twenty equi-sized n segments where each segment contain fifty sampled values, ( N = 50), and the mean value of these data represent the PAA coefficient for each segment. Hence, the PAA coefficients represent the local mean as defined in the proposed algorithm at subset (a) of STEP 4. For both signals, the PHPaa_L and PHPaa_U indicators, as defined in lines 8 and 17 of the subset (b) in the STEP 4, are depicted on the bottom plot of the Fig. 11. As we can remark, these indicators shows almost the same range for both cases. Indeed, the value of MaxUpVariation (res. MaxDownVariation) when the incremental value ΔPHPaaL (res. ΔPHPaaU ) becomes negative occurs at the IndexMaxUpVariation=8 (res. IndexMaxDownVariation=14) for both cases. The latter results mean that the interval of the ninth segment, [400, 450], include a maximum positive change in the mean of the

4. Case study: simulation results Both QTA and QD are tested on a quarter car active suspension 65

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

M. Taktak et al.

Fig. 19. Reachability matrix of the quarter-car model operating at mode 1.

Fig. 20. Reachability matrix of the quarter-car model operating at mode 2.

quarter-car active suspension system shown in Fig. 12. In hydraulic part, we start drwing the HBG model by assigning 0-junction to represent each points with different pressures. In hydraulic actuator ofFig. 12, the distinct pressure points in the two cylinder chambers of the rod piston are identified and 0-junctions (02 and 03) are used to represent them. In addition, a 0-junction (01) is used to represent the input pressure in the spool valve (see Fig. 12). In next step, we insert the component models between appropriate 0-junction pairs using 1junction. For this reason, a relative pressure between the two cylinder chambers is formed using a 1-junction (11). A resistive element, R, is

system shown in Fig. 12. This system is a standard vehicle dynamic model used in automotive engineering and consists of the sprung mass Ms (car body) and the unsprung mass Mu that accounts for the wheel axle masses supported by the tire. The active suspension is modeld as a spring and a dumper with hydraulic actuator in parallel, which connect the unsprung to the sprung mass. The tire is modeled as a spring and represents the transfert of the road profile to the unsprung mass. A more extended description of the suspension systems and their corresponding continuous BG model can be found in Chapter 6 of this reference (Marzouki et al., 2013). Fig. 13 depicts the HBG model of the

66

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

M. Taktak et al.

− − Fig. 21. Forward tracing for parametric faults:{MC2+, MC1−, IMs , IMu , MR2−, MR1+}.

state variable is αi = 1 and when the junction is OFF, the state variable is zero. Thus, the mode of the hydraulic actuator is given by the discrete state vector [α1 α2]. For example, the defining relationships of controlled junctions X11 and X13 are given by following equations:

connected to the 1-junction (11) and represents the internal leakage of the cylinder. Notice that external leakage flow is assumed to be negligible. In this study, switched phenomena is the autonomous mode change in the four way spool valve. In fact, the spool valve flow rate dynamic equation has two modes depending on the valve position as given in the following equation:

⎧ ⎡ ⎪Q1 = Cd ω x v ⎢sw1 (x v ) ⎪ ⎣ ⎨ ⎡ ⎪ ⎪Q2 = Cd ω x v ⎢⎣sw1 (x v ) ⎩

⎤ − Pr ) ⎥ ⎦ ⎤ 2 2 (P − Pr ) + sw2 (x v ) ρ (Ps − P2 ) ⎥ ρ 2 ⎦

2 (P ρ s

− P1) + sw2 (x v )

2 (P ρ 1

(13)

⎧1ifx v (t ) ≥ 0 , sw2 (x v ) = 1 − sw1 (x v ), x v (t ) is the displawhere sw1 (x v ) = ⎨ ⎩ 0otherwise cement of spool valve, Cd is the discharge coefficient, ω is the spool valve area gradient and ρ is the fluid density. For each mode, two controlled 1-junction ( X11,4 and X12,3) are used to model the flow rate direction in the spool valve. The controlled signals for turning these junctions ON and OFF are generated by the two finite state automata shown in the bottom of the Fig. 13. A binary state variable αi(i = {1, 2}) is associated with each controlled junction, when the junction is ON its

⎧ α1 f4 = α1 f2 = α1 f5 ⎪ ⎨ α1 (e4 − e2 + e5) = 0 ⎪α f = α f = α f = 0 1 4 15 ⎩ 13

(14)

⎧ α2 f14 = α2 f3 = α2 f13 ⎪ ⎨ α2 (e14 − e3 + e13) = 0 ⎪α f = α f = α f = 0 2 11 2 12 ⎩ 2 10

(15)

The second-order transfer function between the spool valve displacement xv and supplied voltage u is given in bloc “valve” as:

x v (s ) = u (s )

Kv s2 ωv2

+

2ξ s ωv

+1

(16)

where Kv is the gain of the servo-valve, ωv is the natural frequency of the servo-valve and ξ is the damping ratio of the servo-valve. Derived from Eq. (13), the modulated resistor elements MR1, MR2, MR3 and 67

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

M. Taktak et al.

⎧ ⎛ ⎪ α2⋅f10 = ΦR2 (x v , e10 , α2 ) = α2⋅⎜Cd⋅ω⋅x v (t )⋅ ⎨ ⎝ ⎪ ⎩ α2 ⋅f10 = 0

Table 4 Symbolic Fault Signature Matrix of the quarter-car active suspension system. F/R

Res1: e16

Res2: e7

Res3: f22

Res4: f26

F1: MC1–

α1 {+,?} α2 {+,?}

{0, +}

{0, +}

{0, +}

F2: MC2+

{0, −}

α1 {−,?} α2 {−,?}

{0, +}

{0, +}

F3: IMs–

{0, +}

{0, −}

{−,?}

{0, +}

F4: IMu

–

F5: IMu

+

⎞ e10 ⎟ = α2⋅MR2⋅ e10 ⎠

(18)

⎧ ⎛ ⎪ α1⋅f14 = ΦR3 (x v , e14 , α2 ) = α2⋅⎜ −Cd⋅ω⋅x v (t )⋅ ⎨ ⎝ ⎪ ⎩ α2 ⋅f14 = 0

2 ρ

⎞ e14 ⎟ = −α2⋅MR3⋅ e14 ⎠ (19)

⎧ ⎛ ⎪ α1⋅f18 = ΦR 4 (x v , e18, α1) = α1⋅⎜Cd⋅ω⋅x v (t )⋅ ⎨ ⎝ ⎪ ⎩ α1 ⋅f18 = 0

{0, −}

{0, +}

{+,−}

{+,?}

{0, +}

{0, −}

{−,+}

{−,?}

F6: MR1+

α1 {0, −}

α1 {0, −} α1 α2 {0, −,+}

α1 {0, +}

α1 {0, +}

F7: MR2–

α2 {0, −}

α2 {0, −} α1 α2 {0, −,+}

α2 {0, +}

α2 {0, +}

F8: MR3+

α2 {0, −,+} α1 α2 {0, −,+}

α2 {0, −}

α2 {0, −}

α2 {0, −}

F9: MR4–

α1 {0, −} α1 α2 {0, −,+}

α1 {0, −}

α1 {0, −}

α1 {0, −}

2 ρ

2 ρ

⎞ e18 ⎟ = α1⋅MR 4⋅ e18 ⎠ (20)

From the above equations, it is clear that resistor elements are modulated by the displacement of spool valve xv(t) as shown in hydraulic part of the HBG model with additional information bond from the output of the valve. In addition, the Boolean parameters α1 and α2 selects the set of equations that are valid given the discrete state vector [α1 α2]. Furthermore, the parameter representing C-element of the cylinder is modulated by the velocity of the quarter car body, vs (t ), because the fluid compliance of the two chambers of the cylinder is not constant and writing as:

⎧ V1 − A1 ∫ xṡ dt ⎪ MC1 = β ⎨ ⎪ MC = V2 + A2 ∫ xṡ dt 2 β ⎩

MR4 (see Fig. 13) represent the non linear relationship between flow rate and chamber pressure P1 and P2. The constitutive relation of these modulated resistors are given as:

⎧ ⎛ ⎪ α1⋅f4 = ΦR1 (x v , e4 , α1) = α1⋅⎜ −Cd⋅ω⋅x v (t )⋅ ⎨ ⎝ ⎪ ⎩ α1 ⋅f4 = 0

2 ρ

⎞ e4 ⎟ = −α1⋅MR1⋅ e4 ⎠

(21)

where V1 and V2 are the original total control volumes of the two cylinder chambers, A1 and A2 are the ram area of the two chambers and β is the effective bulk modulus of the hydraulic fluid. Hence, depending on the valve position, the flow into the cylinder moves the piston according to the mass conservation law on each chamber. Constitutive

(17)

− − Fig. 22. Time evolution of residual signals for simulated parametric faults:{MC2+, MC1−, IMs , IMu , MR2−, MR1+}.

68

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

M. Taktak et al.

Table 5 Predefined parameters algorithm used in case study runs. Window size N1

Sliding size N2

Dimensionality Reduction N

Detection thresholdλ

λ Discontinuity threshold Λ

L

Classification threshold ths

100

20

25

7

0.3

20

7e−4

fault set. Hence, the set of parametric faults take the following form: {MCi+, MCi−, MRj+ , MRj−, Ms+, Ms−, Mu+, Mu−} where i = {1, 2} and j = {1, 2, 3, 4}. In the case of the modulated parameters, the abrupt change value is applied in a constant positive parameter. For example, an abrupt increasing change in the modulated capacitance value MC2 , denoted as MC2+ , can be expressed, from equation (21.2), as

relations of 0-juntions, 02 and 03, represent this conservation of mass on each chamber. The 2-ports transformer element, TF1 and TF2 , are needed to couple the mechanical and hydraulic energy domains. They represent the relationship between velocity of the single rod piston and the flow rate changes of two chambers. The effort source, Se:Ps, which represents the supply pressure for the hydraulic actuator, is attached to the 0-junction, 01. Finally, De:P1 and De:P2 represent the pressure sensors mounted in the two cylinder chambers of the rod piston. The mechanical part can be treated in a manner quite similar to the hydraulic part (Karnopp et al., 2000). However, for mechanical part, 1 1 see Fig. 13, elements C: K and C: K represent the compliance of two t s R : C springs, s denotes the damper coefficient, I : Ms and I : Mu represent the inertias of the sprung mass and unsprung mass, respectively. In addition, elements Df : xṡ and Df : xu̇ represent velocity sensor mounted in sprung and unspring mass, respectively. Fig. 13 depicts the HBG model of the quarter car active suspension system where all controlled 1-junction are assigned with their preferred causality as required by the HSCAP procedure and all storage elements are in preferred integral causality. The coupling of the measured variable obtained from sensor in real process and the output of the unified HBG-model with nominal parameters is depicted in Fig. 14 where residual sinks are denoted by rSf and rSe . Note that nominal parameters values, considered in this investigation, are given in Table 3. From the model given in Fig. 14, the PTCG of the quarter-car active suspension model is systematically build according to transformation rules given in Fig. 3 in which the arcs with dotted and highlighted solid line represent directed edges between flow and effort nodes respectively. However, constitutive relations of elements in HBG are shown, in Fig. 15, with directed arcs in solid line. Circled vertices, in the PTCG, represent the measured variables and labels on the directed edges capture the causal relation in effect between the variables. Note that the numbered bonds, in the HBG model of Fig. 13, are converted to corresponding variables with subscript indicating the bond number. The edges of controlled junctions are parametrized by their related binary state variables (for visibility purpose, binary state variable are displayed inside a square box in the Fig. 15). MATLAB is used to simulate the HBG model of the switched quarter car active suspension system which is coupled with a faulty model that emulate the real system and generate measured variables. To simulate reality, different level of AWGN is added to sensor outputs from the faulty HBG model. This model was simulated using a fixed-step simulation with sample period of ts = 0.001s where input 2π valve signal is u (t ) = 2 sin( 0.4 ts ) and the vehicle is assumed to be initially traveling a “smooth” curb before reaching a flat road profile. The road profile, the input valve signal and the state of the two switches sw1 and sw2 are shown in Fig. 16. Hence, the switching signal is known and switched system operates alternatively between mode 1 and 2, where mode1=[α1 = 1 , α2 = 0] and mode2=[α1 = 0 , α2 = 1]. The simulation model was run for 0.8s of simulation time with solver “ode4″. Let us consider, in this case study, all parameters of the hydraulic part, i. e. {MC1, MC2, MR1, MR2 , MR3, MR 4}, and the two inertia parameter, i. e. {Ms, Mu}, in the mechanical part of the quarter-car active suspension system as the set of parametric faults. So, in the faulty HBG model, abrupt increase change, denoted as (+), for each parameter value in the parametric fault set, is added to simulate a single abrupt fault occurrence which is injected at 0.35s in a first run and at 0.45s in a second run. A similar work is made to simulate an abrupt decrease change, denoted as (−), for each parameter value in the parametric

⎛ ΔV ⎞ V2 ⎜1 + 2 ⎟ + A2 ∫ xṡ dt V ⎠ ⎝

2 MC2+ = where ΔV2 represents the magnitude of the β fault. Therefore, from the coupled HBG model, the numerical values of the residual sinks are obtained for each parametric fault set in both case of abrupt increase and decrease applied in mode 1 as well as in mode 2. If we try to follow only a quantitative diagnosis approach in order to study the parametric fault isolability of the quarter-car active suspension system, then a so-called all-mode Fault Signature Matrix (FSM) should be available from the Diagnosis HBG model of the switched system. Note that all-mode FSM is an extension of a FSM (Tagina et al., 1995), used for Diagnosis BG model of continuous system, that contains discrete switch state variables. However, when the residuals are in numerical form, the all-mode FSM can only obtained from the analysis of the causal paths to each residual sinks in the faultless HBG model (Borutzky, 2012). Instead of inspecting the causal path from the HBG model, the proposed PTCG provides a useful tool to speed up the analysis of the causal paths because the PTCG is a parametrized labeled directed graph that explicitly represent the causal and temporal relations between system variables. In graph theory, it is possible to represent a directed graph by an adjacency matrix where the rows and columns both represent the effort and flow variables (Yang et al., 2014; Jiang et al., 2008). The adjacency matrix of the quarter car PTCG at mode 1, denoted A[α1=1, α2=0], and at mode 2, denoted A[α1=0, α2=1], are shown in Figs. 17 and 18 respectively. Note that the (i , j )th entry in adjacency matrix is assigned a value of “1” if there is a directed edge from vertice i to vertice j , otherwise it is assigned a value of “0”. For a PTCG with V vertices and an adjacency matrix A[α1, α2] at a given operating mode, the following matrix:

R[#α1, α2] = (A[α1, α2] + A[2α1, α2] + A[3α1, α2] +…+A[Vα1, α2] )

#

(22)

R[#α1, α2]

is defined as the reachability matrix, where the (i , j )th element of indicates whether there exists any directed path from vertice i to vertice j . The superscript # in Eq. (22) stands for the Boolean equivalent of the matrix R defined by the following relationship:

⎧ 0ifR (i , j ) = 0 R # (i , j ) = ⎨ ⎩1ifR (i , j ) ≠ 0

(23)

The reachability matrix of the adjacency matrix A[α1=1, α2=0] and A[α1=0, α2=1] is shown in Figs. 19 and 20 respectively. It is clear from reachability matrix R[#α1=1, α2=0] or R[#α1=0, α2=1] that a component (displayed in the left hand side of the first column) can reach measured vertices (e7, e16 , f22 andf26 ) if the corresponding entry has a value of “1” which means that a fault in this component contributes to the corresponding residual sink and, then, can be detected. Moreover, we can see from reachability matrix R[#α1=1, α2=0] and R[#α1=0, α2=1] , that the components involved in the measured vertices are identical for the four residual sinks which mean that the component fault signature is similar and, then, the fault cannot be isolated. For more isolation ability and without modifying the number of sensors, the QTA with QD approach is followed as we propose in this paper. Hence, the PTCG is used to predict qualitative effects of parametric faults in residual sinks. To this end, a forward-propagation 69

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

M. Taktak et al.

− − Fig. 23. Results of QTA and QD algorithm for parametric faults: {MC2+, MC1−, IMs , IMu , MR2−, MR1+}.

70

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

M. Taktak et al.

isolated. The QD algorithm confirms this result and returns the single fault candidate {F1: MC1−} which has the GSI value of 0.6. Fig. 23.(c) shows the − . We can see that an abrupt negative QTA-QD results of fault scenario, IMs change of Res3 was detected before the slow increase change in Res4 when the slice reaches the sample 360. Estimated symbolic signature return a {−,+} for Res3, {0, +} for Res4 and {0, 0} for both Res1 and Res2. In step1, the QD algorithm returns the following candidate faults: [{F3}(LSI = 0.6), {F1, F2, F6}(LSI = 0.2), {F4, F8, F9}(LSI = 0)]T where faults are rank ordered in the decreasing order of LSI. In step-2 and for each fault candidate, determined in step-1, GSI is calculated using Eq. (3). Final result of the QD algorithm returns the following candidate faults: [{F3}(GSI = 0.4), {F7}(GSI = 0.2), {F5}(GSI = 0.1), {F1, F2}(GSI = 0)]T , − which mean that for the fault scenario IMs , F3 is the most likely occurred − with GSI value of 0.4. Unlike the fault scenario IMs , the simulated abrupt − decrease in unsprung inertia IMu shown in Fig. 23.(d) return a positive abrupt increase in magnitude of Res4 which is firstly detected at slice position 480. In our proposed algorithm it is assumed that “abrupt” change in residual sinks cannot occur simultaneously. However, the predicted − qualitative signature for F4: IMu given in Table 4 as well as the time evolution of residual signal given in Fig. 22 shows a simultaneous “+” abrupt increase in the magnitude of Res3 and Res4. From results shown in table of the Fig. 23.(d), we can observe that ΔUp of Res3 is quite higher than the predefined threshold Λ and MaxUpVariation of Res3 is slightly less than the predefined threshold λ . Although estimated qualitative − signature of Res3 has “no change” in the zeroth order change, IMu fault − } scenario is isolated but QD algorithm return single fault candidate {F4: IMu with low GSI index value of 0.2. Finally, we observe that estimated qualitative signature for MR2− and MR1+, shown in tables of Fig. 23.(e) and (f) respectively, are perfect similar where no “abrupt” change was detected. However, on-line algorithm results shows that a positive change was firstly detected in Res3 at slice position 660 for MR2− and in Res4 at slice position 480 for MR1+. Because the abrupt decrease in MR2 was injected at 0.45 s, i.e. when the system is at operating mode 1, the binary state variable α2 = 0 and consequently the controlled junction X12 is OFF. That is, no change was expected in residual sinks until system switch from mode 1 to mode 2 at 0.6 s (see Fig. 16). Likewise, abrupt increase in MR1+ was injected at 0.35 s i.e. when the system is at mode 2; thus the binary state variable α1 = 0 and controlled junction X11 is OFF. Hence, no change was expected in residual sinks until system switch from mode 2 to mode 1 at 0.4 s as depicted in Fig. 16. Although estimated qualitative signature are similar for MR2− and MR1+ , but they occurs at opposite operating mode. Therefore, MR2− is detectable and isolable (GSI=1 for fault candidate F7) when system is at operating mode 2 and MR1+ is detectable and isolable (GSI=1 for fault candidate F6) when system is at operating mode 1.

algorithm (Mosterman and Biswas, 1999; Narasimhan, 2002; Daigle, 2008) generates the fault signatures by propagating the qualitative effect of a fault," + " or " − ", to all residual sinks vertices in the PTCG. All deviation propagation start off as zeroth order effects and when a temporal integrating edge is traversed order of change increment, that is zeroth order becomes a first order change and so on. In addition, deviation of change will be dependently of the switch state when a parametrized edge is traversed. Fig. 21 shows the forward propagation {MC2+, MC1−, corresponding to the parametric faults: − − IMs , IMu , MR2−, MR1+} . Note that symbol “? ”, in signature of residual sinks, mean that change deviation is qualitatively indeterminate. Let's consider the forward propagation of the fault MC2+ along the sub path marked by A in the PTCG of Fig. 15, 1 dt C2

− MC2+ → e7− → α2 f10 →α2 f 7+ ⎯⎯⎯⎯→α2 e7−? . On the one hand, the first change of the measured variable e7 is dependently of the switch state α2 , and on the other hand this change is qualitatively indeterminate, "? ", because f7 → e7 propagate a + effect, while 1 (the inverse of MC2+ ) propagate a C2 − effect. Table 4 shows the SFSM which summarize the all-mode fault signatures of the fault components set generated by the forward propagation along the PTCG. Note that a predicted signature of the form {0, 0, +} (res. {0, 0, −}) will be reported as {0, +} (res. {0, −}). For our proposed on-line algorithm developed to detect abrupt faulty change and estimate the lower-order symbolic signature, a fixed window size of length N1 is defined as a slice and each subsequent slice is created by moving the window in time a fixed distance of N2 such that N2 < N1. The algorithm was tested on data set, Y = {Res1, Res2, Res3, Res 4}, of simulation results obtained by numerical evaluation of the coupled quarted-car model where single parametric fault is added in faulty behavior model as above mentioned. Fig. 22 shows the simulation results of residual sinks, for the − − , IMu , MR2−, MR1+}, corresponding parametric faulty set: {MC2+, MC1−, IMs where in all the plot, the time is expressed in samples, and recall that the sample time is 0.001s . Using the parameters given in Table 5, each data was converted into individual slice of 100 samples by using a sliding window that moved in 20 samples and data within slice were subdivided into 25 equi-sized interval. Fig. 23 shows the time evolution of all residual sinks with PHPaa_U and PHPaa_L indicators when the first change is detected. Particularly, Fig. 23.(a) show the results obtained for the fault scenario where abrupt increase in the modulated capacitance, MC2+ , is injected at 0.35 s (i.e. where the system is at mode 2). As we can see in table of QTA-QD results, our algorithm detect a first change in residual sink Res2 when the slice reach the sample 380 (the PHT with PAA indicators are marked in bold case) and change point detection fall within the interval [348352]. The symbolic signature of the first change in Res2 is estimated as: S 2 = {−,0}*. Furthermore, the estimated symbolic signature of all remaining residual sinks is given as described in substep (d) of STEP 4 of our proposed QTA algorithm. From this estimated fault signature, FS = [{0, −}, {−,0}* , {0, +}, {0, +}], the QD component return the single fault candidate {F2: MC2+} from the SFSM with LSI=0.6 (step-1) and GSI=0.6 (step-2). Thus, the fault scenario MC2+ is detected and isolated. Fig. 23.(b) depicts the results of run 2 where the faulty decrease in modulated capacitance, MC1−, is simulated at 0.45 s (i.e. where the system is at mode 1). Likewise in above case, table below the plot graph of Fig. 23. (b) shows results of on-line QTA-QD algorithm when a first change is detected. Therefore, an abrupt (i.e. with discontinuity) positive increase of Res1 was firstly detected over the interval [448452] when the slice reaches the sample 480. The estimated symbolic signature shows a “+” first order change in Res3 and Res4 and no change, “0”, in the slope of Res1 and Res2. As the first order change in the predicted signature of Res1 is “? ” for both mode 1 and 2, as we can see in Table 4 for the faulty case F1: MC1−, then this uncertainly change can be “+”, “− ” or “0”. The resulting estimated signature of Res1, Res2 and Res3 has a most likely match with MC1− in the predicted signature Table 4. Hence, this fault can be detected and, also, be

5. Conclusion In this work, we present a QTA with QD procedure for abrupt parametric faults diagnosis of switched system. QTA is a Qualitative Trend Analysis applied on residual signals which are numerically evaluated by mean of residual sinks added to a Hybrid Bond Graphmodel. QTA analysis is used for on-line detection of abrupt change based on combination of PAA and PHT. The QD is a Qualitative Diagnosis procedure which is performed in two stages. In the first offline stage, we exploit the causal and temporal relations between process parameters and the residual sinks provided by the PTCG in order to generate a Symbolic Fault Signature Matrix. In the second online stage, trend extraction by linear regression is triggered after change detection in order to estimate the lower order time-derivative symbol for each residual sinks. Finally, a stepwise similarity measure between estimated and predicted symbolic fault signature is computed to isolate the most likely occurred fault. The proposed approach was tested on a simulated switched quarter-car active suspension system. 71

M. Taktak et al.

Engineering Applications of Artificial Intelligence 59 (2017) 51–72

References

Maurya, M.R., Rengaswamy, R., Venkatasubramanian, V., 2007. Fault diagnosis using dynamic trend analysis: a review and recent developments. Eng. Appl. Artif. Intell. 20, 133–146. Mosterman, P.J., 2001. Diagnosis of physical systems with hybrid models using parametrized causality. Institute of robotics and mechatronics, Germany. Mosterman, P.J., Biswas, G., 1999. Diagnosis of continuous valued systems in transient operating regions. IEEE Trans. Syst., Man Cybern. 29 (6), 554–565. Mosterman, P.J., Biswas, G., 2000. A comprehensive methodology for building hybrid models of physical systems. Artif. Intell. 121, 171–209. Mouss, H., Mouss, D., Mouss, N., Sefouhi,L., 2004. Test of Page-Hinkley, an approach for fault detection in an agro-alimentary production system, pp. 815–818, Proceedings of the Asian Control Conference, Volume 2 Narasimhan, S., 2002. Model-based diagnosis of hybrid systems (Ph. D. Dissertation). Vanderbilt University. Narasimhan, S., Zhao, F., Biswas, G., Hung, E., 2000. An integrated framework for combining global and local analyses in diagnosing hybrid systems. In Eleventh International Workshop on Principles of Diagnosis, pages 163–170 Niu, G., Zheo, Y., Tung Tran, V., 2015. Fault detection and isolation based on bond graph modeling and empirical residual evaluation. J. Mech. Eng. Sci. 229 (3), 417–428. Ould Bouamama, B., Biswas, G., Loureiro, R., Merzouki, R., 2014. Graphical methods for diagnosis of dynamic systems: review. Annu. Rev. Control 38 (2), 199–219. Precup, R.E., Angelov, P. B., Costa, S. J., Sayed-Mouchaweh, M., 2015. An overview on fault diagnosis and nature-inspired optimal control of industrial process applications, Computers in Industry Pulido, B., Alonso-Gonzlez, C., 2004. Possible conflicts: a compilation technique for consistency-based diagnosis. IEEE Trans. Syst., Man, Cybern. (Part B) 34 (5), 2192–2206. Rebecca, M., Roger F, N., da Cruz Marcelin, F., 2013. Construction and analysis of causally dynamic hybrid bond graphs. Proc. Inst. Mech. Eng. Part I – J. Syst. Control Eng. 227 (3), 329–346. Samantaray, A.K., Medjaher, K., Bouamama, B.O., Staroswiecki, M., Dauphin-Tanguy, G., 2006. Diagnostic bond graphs for online fault detection and isolation. Simul. Model. Pract. Theory 14 (3), 237–262. Sebastiao, R., Silva, M., Rabico, R., Gama, J., Mendonca, T., 2013. Realtime algorithm for changes detection in depth of anesthesia signals. Evol. Syst. 4 (1), 3–12. Serdio, F., Lughofer, E., Pichler, K., Buchegger, T., Efendic, H., 2014a. Residual-based fault detection using soft computing techniques for condition monitoring at rolling mills. Inf. Sci. 259, 304–320. Serdio, F., Lughofer, E., Pichler, K., Pichler, M., Buchegger, T., Efendic, H., 2014b. Fault detection in multi-sensor networks based on multivariate time-series models and orthogonal transformations. Inf. Fusion 20, 272–291. Serdio, F., Lughofer, E., Pichler, K., Pichler, M., Buchegger, T., Efendic, H., 2015. Fuzzy fault isolation using gradient information and quality criteria from system identification models. Inf. Sci. 316, 18–39. Serdio, F., Lughofer, E., Pichler, K., Buchegger, T., Pichler, M., Efendic, H., 2014. Reducing False Positives for Residual-Based On-line Fault Detection by Means of Adaptive Filters, pp. 2803–2808, Proceedings of the IEEE SMC 2014 Conference, San Diego, CA, U.S.A. Shaker, A., Lughofer, E., 2014. Self-adaptive and local strategies for a smooth treatment of drifts in data streams. Evol. Syst. 5 (4), 239–257. Staroswiecki, M., Comtet-Verga, G., 2001. Analytical redundancy relations for fault detection and isolation in algebraic dynamic systems. Automatica 37, 687–699. Strömberg, J. E., Top, J., Södermann, U. 1993. Variable causality in bond graphs by discrete effects. International Conference on Bond Graph modelling, ICBGM’93, Proc. Of 1993 Western Simulation Multi-conference, pages 115-119 Tagina, M., Cassar, J.P., Dauphin-Tanguy, G., Staroswiecki, M., 1995. Monitoring of systems modelled by bond-graphs. International Conference On Bond Graph Modeling And Simulation, ICBGM’95.(vol. 27(1), of Simulation Series, pp.275–280). Las Vegas, Nevada, USA: SCS Publishing Triki, S., Kamoun, A., 2014. EOG classification based on Fuzzyfied symbolic representation. J. Intell. Fuzzy Syst. 26 (6), 2841–2851. Umarikar, A.C., Umanand, L., 2005. Modelling of switching systems in bond graphs using the concept of switched power junctions. J. Frankl. Inst. 342 (2), 131–147. Yang, F., Duan, P., Shah, S.L., Chen, T., 2014. Capturing connectivity and causality in complex industrial process, Springer, SpringerBriefs in Applied Science and Technology

Basseville, M., Nikiforov, I., 1993. Detection of abrupt changes: theory and applications. Prentice-Hall. Borutzky, W., 2009a. Bond graph model-based fault detection using residual sinks. Proc. Inst. Mech. Eng., Part I: J. Syst. Control Eng. 223, 337–352. Borutzky, W., 2012. Bond-graph-based fault detection and isolation for hybrid system models. Proc. Inst. Mech. Eng., Part I: J. Syst. Control Eng. 226 (6), 742–760. Borutzky, W. 2009. Bond Graph Modelling and Simulation of Multidisciplinary Systems An Introduction Simulation Modelling Practice and Theory, Volume/Issue 17/1, pp. 3-21 Chiang, L.H., Russell, E.L., Braatz, R.D., 2001. Fault Detection and Diagnosis in Industrial Systems. Springer, London Berlin Heidelberg. Daigle, M.J., 2008. A qualitative event-based approach to fault diagnosis of hybrid systems (Ph. D. Dissertation). Vanderbilt University. Dash, S., Rengaswamy, R., Venkatasubramanian, V., 2003. Fuzzy-logic based trend classification for fault diagnosis of chemical processes. Comput. Chem. Eng. 27 (3), 347–362. Dauphin-Tanguy, C.S., Rombaut, C., 1989. Bond graph approach of commutating phenomenon. In: Husson, R. (Ed.), Proceedings of Advanced Info. Processing in Automatic Control. Nancy, IFAC, 339–343. Ding, S.X., 2008. Model-based Fault Diagnosis Techniques: Design Schemes, Algorithms, and Tools. Springer, Berlin Heidelberg. Gelso, E.R., Castillo, S.M., Armengol, J., 2007. Diagnosis based on interval analytical redundancy relations and signs of the symptoms, In IFAC International Workshop on Intelligent Manufacturing Systems IMS’07, Alicante, Spain Gertler, J.J., 1998. Fault detection and diagnosis in Engineering Systems. Marcel Deckker, New York, NY, USA. Isermann, R., 1984. Process fault detection based on modelling and estimation methods A survey. Automatica 20–4, 384–404. Isermann, R., 2009. Fault Diagnosis Systems: An Introduction from Fault Detection to Fault Tolerance. Springer, Berlin Heidelberg. Jiang, H., Patwardhan, R., Shah, S. L., 2008. Root Cause Diagnosis of Plant-wide Oscillations using the Adjacency Matrix, In International Workshop on Automatic Control (IFAC), Seoul, Korea Karnopp, D.C., Rosenberg, R.C., Margolis, D.L., 2000. System Dynamics: A Unified Approach 3rd ed. John Wiley, New York, NY, USA. Keogh, E., Pazzani, M., 2000. A Simple dimensionality reduction technique for fast similarity search in large time series databases, In: Proceedings of the 4th Pasific-Asia Conference on Knowledge Discovery and Data Mining, pp. 122-133 Korbicz, J., Koscielny, J.M., Kowalczuk, Z., Cholewa, W., 2004. Fault Diagnosis–Models, Artificial Intelligence and Applications. Springer Verlag, Berlin Heidelberg. Kypuros, J.A., 2001. Variable Structure Model Synthesis for Switched Systems (Ph.D. thesis). University of Texas at Austin. Kypuros, J.A., 2013. System Dynamics and Control with Bond Graph Modeling. Taylor & Francis Group. Low, C.B., Wang, D., Arogeti, S., Luo, M., 2010a. Quantitative hybrid bond graph-based fault detection and isolation. IEEE Trans. Autom. Sci. Eng. 7 (3), 558–569. Low, C.B., Wang, D., Arogeti, S., Zhang, J.B., 2010b. Causality assignment and model approximation for hybrid bond graph. IEEE Trans. Autom. Sci. Eng. 7 (3), 570–580. Low, C. B., Wang, D., Arogeti, S., Zhang, J. B., 2008. Monitoring ability and quantitative fault diagnosis using hybrid bond graph, in International Federation of Automatic Control (IFAC), pp. 10516-10521, Seoul Korea Luo, M., 2006. Data-Driven Fault Detection using Trending Analysis (PhD thesis). Department of Electrical and Computer Engineering, Louisiana State University and Agricultural and Mechanical College, Lousiana. Mallat, S., Hwang, W.L., 1992. Singularity detection and processing with wavelets. IEEE Trans. on Information Inf. Theory 38, 617–643. Manders, E. J., Narasimhan, S., Biswas, G., Mosterman, P. J., 2000. A combined qualitative/quantitative approach for fault isolation in continuous dynamic systems, Proc. Of- the 4th IFAC safeprocess symposium, Budapest, Hungary, pages 1074-1079 Manders, E. J., Mosterman, P. J., Biswas, G., 1999. Signal to Symbol Transformation Techniques for Robust Diagnosis in TRANSCEND, in Proceedings of the Tenth International Workshop on Principles of Diagnosis (DX99), pp. 155-165, Lock Awe Hotel, Scotland, June 8-11, Marzouki, R., Samantary, A.K., Pathak, P.M., Ould Bouamama, B., 2013. Intelligent Mechatronic Systems. Springer, Berlin Heidelberg.

72