Sequential fault diagnosis for mechatronics system using diagnostic hybrid bond graph and composite harmony search

Abstract

This article proposes a sequential fault diagnosis method to handle asynchronous distinct faults using diagnostic hybrid bond graph and composite harmony search. The faults under consideration include fault mode, abrupt fault, and intermittent fault. The faults can occur in different time instances, which add to the difficulty of decision making for fault diagnosis. This is because the earlier occurred fault can exhibit fault symptom which masks the fault symptom of latter occurred fault. In order to solve this problem, a sequential identification algorithm is developed in which the identification task is reactivated based on two conditions. The first condition is that the latter occurred fault has at least one inconsistent coherence vector element which is consistent in coherence vector of the earlier occurred fault, and the second condition is that the existing fault coherence vector has the ability to hide other faults and the second-level residual exceeds the threshold. A new composite harmony search which is capable of handling continuous variables and binary variables simultaneously is proposed for identification purpose. Experiments on a mobile robot system are conducted to assess the proposed sequential fault diagnosis algorithm.

Keywords

Sequential fault diagnosis diagnostic hybrid bond graph composite harmony search intermittent fault coherence vector mobile robot

Introduction

The growing demand for reliability and safety in modern industrial systems has led to the development of fault diagnosis methodologies. Faults need to be detected close to their occurrence time, so that corrective actions can be taken in a timely manner to avoid catastrophic consequences. In general, fault diagnosis methods can be classified into two categories: model-based and data-driven methods. Model-based method usually tries to compare the behavior of the mathematical model representing the monitored system and the behavior of the actual system. Any discrepancy between these two behaviors may indicate the fault occurrence. Several commonly used model-based methods are observer, analytical redundancy relation (ARR), and parameter estimation–based methods.

Model-based method is able to incorporate a physical understanding of the system for diagnosis purpose.¹ As the understanding of the monitored system improves, the model can be modified to increase its accuracy and address subtle performance problems. Unlike model-based method which requires accurate mathematical model, data-driven method typically attempts to obtain feature information from large amount of historical process data. The main advantage of data-driven method is that it can capture complicated phenomenon without a priori knowledge.²

For modeling complex system, bond graph (BG) is an efficient tool which has exhibited its superiority in fault diagnosis in continuous systems. In Djeziri et al.,¹ a BG-based robust fault detection and isolation (FDI) is developed for the vehicle traction system. This method is able to deal with structured and unstructured uncertainties using the model in the linear fractional transformation form. Robust ARRs with adaptive thresholds are generated to facilitate FDI. A method for diagnosability analysis and structural verification of fault recoverability using the BG tool is proposed in Loureiro et al.³ An approach that structurally states for which faults the system remains able to achieve its objectives without any calculations is developed. However, the works in Djeziri et al.¹ and Loureiro et al.³ only consider single fault scenario and do not concern fault modes. In Calderaro et al.,⁴ a method for failure identification and detection in smart grids using the Petri net is developed. The faults in the power distribution network and failures in the data transmission systems are considered. However, detection of parametric faults is not discussed.

In real-world applications, many systems, for example, high-speed printer, vehicle, and switched mode power converter, exhibit both continuous-time and discrete-event dynamics, and these systems are called hybrid systems. In the past few years, fault diagnosis of hybrid systems has received significant attention due to their wide applications.^5–10 In Narasimhan and Biswas,⁵ a methodology for FDI in hybrid systems, named Hybrid TRANSCEND, is developed. This method handles a single and abrupt parametric fault using qualitative and quantitative frameworks. The current mode of the systems is estimated by a hybrid observer. This observer combines a Kalman filter, which estimates state variables within a current mode and a mode change detector to decide which mode is the current mode. However, this method cannot deal with unknown fault modes.

A method for real-time mode tracking in the presence of fault is developed for hybrid systems in Arogeti et al.⁷ This method is based on the single fault assumption, meaning that once the system detects a fault, the certainty ARR is determined and will not be changed to uncertainty ARR during the monitoring process; this is not true for multiple faults condition. As a result, this approach is not able to handle mode tracking under multiple faults. In Arogeti et al.,⁸ a method for dealing with both parametric faults and fault modes is proposed. This method is able to single parametric fault with abrupt and incipient nature. However, no issue related to multiple distinct faults is addressed. A mode identification using observer for a switching linear system is developed in Li et al.¹⁰ The problem with this approach is that no parametric fault is considered.

To the best of our knowledge, no work related to fault diagnosis of asynchronous distinct faults, in which the earlier occurred fault may mask the latter occurred fault in terms of fault symptom (i.e. coherence vector (CV)), has been reported. In this article, a new sequential fault diagnosis method is developed to deal with the aforementioned issues. In the developed framework, the identification task is reactivated based on two conditions. In order to facilitate fault identification task after fault isolation, a new composite harmony search (CHS) which is capable of handling continuous variables and binary variables simultaneously is proposed.

The main contributions of this work are fourfold:

A fault detection method to deal with asynchronous distinct faults is developed by introducing the new concept of second-level residual (SLR). This method is able to detect the latter occurred fault whose signature is masked by the one of the earlier occurred faults.

A unified fault formulation for abrupt fault and intermittent fault, which enables estimations of fault characteristics (e.g. for intermittent fault, the frequency of repetitions, and magnitude), is proposed. This fault modeling technique allows the developments of the fault diagnosis method when the fault-type information cannot be obtained a priori.

A sequential fault identifier (SFI) scheme, which can identify the fault information as early as possible for future maintenance purpose, is developed.

A new harmony search (HS) algorithm called CHS is developed. The composite stems from its capability of finding real-valued variables together with binary-valued variables.

The rest of this article is organized as follows: in section “DHBG and AGARRs for fault diagnosis,” the concepts of diagnostic hybrid bond graph (DHBG) and augmented global analytical redundancy relations (AGARRs) for hybrid system diagnosis are introduced, section “Proposed sequential identification for asynchronous distinct faults” presents the proposed sequential fault diagnosis method for asynchronous distinct faults, section “Experimental study” describes the experiment process in a mobile robot system and also analyzes the experimental results, and finally, section “Conclusion” concludes this article.

DHBG and AGARRs for fault diagnosis

For model-based fault diagnosis, the modeling is a crucial step where the developed model should be able to sufficiently capture the physical phenomenon of the monitored system and at the same time easy to analyze. BG is a useful modeling tool and is originally geared towards modeling of continuous systems. To model a hybrid system using BG language, additional features are necessary to capture the discrete mode changes of a hybrid system. Hybrid bond graph (HBG) extends BG modeling by incorporating the controlled junctions to enable the hybrid system to be modeled using the BG components.⁵ With the aid of controlled junctions, the hybrid systems can be represented by HBG in a compact manner. In HBG, changes in configuration through the various operating modes of the system can result in a need to reassign the causality. In order to overcome this problem, a systematic procedure called the Sequential Causality Assignment Procedure for Hybrid Systems (SCAPH) is developed.⁶ Equipped with the SCAPH, the DHBG is developed where a set of dynamic constraints called AGARRs can be derived to describe the behaviors of a hybrid system at all operating modes in a unified manner. Without loss of generality, AGARR equations take the following form

AGAR R_{i} (θ, a, β, De, Df, u) = 0, i = 1, 2, \dots, m

(1)

where $θ \in R^{n_{θ}}$ is the nominal parameter vector of DHBG components, $a \in R^{n_{a}}$ is a binary vector representing the state of controlled junctions, $β \in R^{n_{β}}$ is the efficiency factor vector for sensors and actuators, $De \in R^{n_{De}}$ is the effort sensor vector, $Df \in R^{n_{Df}}$ is the flow sensor vector, $u \in R^{n_{u}}$ is the known input vector, and $m$ is the number of independent AGARR equations derived from the DHBG.

The introduction of efficiency factor vector $β$ in AGARRs enhance the capability of global analytical redundancy relations (GARRs) to identify the fault severity of system components (e.g. sensors and actuators) which cannot be described by physical parameters. As a result, $n_{β} = n_{De} + n_{Df} + n_{u}$ in equation (1). Since in the DHBG, sensor fault is modeled by multiplying the normal measurement by its corresponding efficiency factor, normal measurement of sensors is represented by dividing $De$ or $Df$ by its corresponding efficiency factor in AGARR equations described by equation (1). As for the actuator fault in DHBG, it keeps the same form as its fault modeling (i.e. actuator fault is modeled by multiplying the input by its corresponding efficiency factor) in equation (1) because AGARR equations require actual input signals.

According to the derived AGARR equations, fault signature matrix (FSM) at different operating modes, referred to as mode-dependent fault signature matrix (MD-FSM), can be established for monitorability analysis. Each entry of the MD-FSM table holds a Boolean value which is defined as

M_{i, j}^{a_{l}} = {\begin{matrix} 0 & \frac{\partial A G A R R_{i} (a_{l})}{\partial σ_{j}} = 0 \\ 1 & \frac{\partial A G A R R_{i} (a_{l})}{\partial σ_{j}} \neq 0 \end{matrix}, l = 1, \dots, 2^{n_{a}}, j = 1, \dots, n_{θ} + n_{β}

(2)

where $σ \in R^{n_{θ} + n_{β}}$ is the considered fault vector including both parametric fault and nonparametric fault, and $a_{l}$ is the lth mode.

For hybrid systems, residual (i.e. numerical evaluation of AGARRs) inconsistency is not necessarily caused by a component fault, but may also be caused by an unexpected change of mode. This unexpected change of mode is referred to as fault mode where the faulty state is known a priori and can be modeled by known parameters only. For example, short-circuit and open-circuit faults of power switches are commonly fault modes in power electronics. Similar to FSM, mode change signature matrix (MCSM) can be built from equation (1) which represents cause–effect relations between mode changes and residuals. The MCSM is used for mode change isolation.⁷ The component fault and mode change detection is carried out by evaluating the residuals, and a fault or a mode change is declared if one of the residuals exceeds the predefined threshold. A CV is defined as $C = [c_{1} c_{2} \dots c_{m}]$ to indicate the inconsistency of residuals and $c_{i} = 1, i = 1, 2, \dots, m$ ; if $g_{i}$ is inconsistent (i.e. exceeds the threshold) and zero otherwise. The binary CV will be zero when the system is fault free (or no mode change), and a nonzero CV means a faulty condition or an unknown mode change.

Proposed sequential identification for asynchronous distinct faults

After a nonzero CV is observed, the next step is to isolate the component fault or unknown mode change by comparing the CV with MD-FSM or MCSM. The isolation will lead to a set of fault candidates or suspected mode changes. In order to obtain the knowledge of fault severity, fault identification is carried out, and the identification result can also refine the isolation results. In this article, multiple distinct faults are considered, and these faults may occur in the monitored system in an asynchronous manner. The fault type could be abrupt or intermittent for component fault, and the knowledge of fault type is unknown beforehand.

Fault modeling

For fault identification, a certain form of fault modeling is required. In this article, a fault is characterized by the following form

\bar{F} (t, t_{0}) = {\begin{matrix} F \\ η_{F} \cdot S_{ab} + (1 - η_{F}) \cdot S_{int} \end{matrix} \begin{matrix} if t < t_{0} \\ if t \geq t_{0} \end{matrix}

(3)

with

\begin{matrix} S_{ab} = F_{ab} \\ S_{int} = \frac{F_{int} - F}{2} sign [\sin (\frac{2 π t - 2 π t_{0}}{T_{F}})] + \frac{F_{int} + F}{2} \end{matrix}

where $F$ is the nominal value of a component or an efficient factor, $η_{F}$ is the fault-type index, $F_{ab}$ is the abrupt fault amplitude, $F_{int}$ is the intermittent fault amplitude, $T_{F}$ is the intermittent fault period, and $t_{0}$ is the fault occurrence time.

In equation (3), $η_{F}$ is the fault-type index which is a Boolean variable. If $η_{F} = 1$ , the fault is an abrupt fault, and the fault is described by an intermittent process when $η_{F} = 0$ . The advantage of the fault profile defined in equation (3) lies in that it is able to describe different fault profiles for various faulty components in a compact manner.

Sequential fault diagnosis

For the case of asynchronous faults, the earlier occurred fault leads to fault symptom (i.e. a nonzero CV) which enables fault detection, isolation, and identification. After the identification process is finished, another fault occurs, and this latter occurred fault will have two possible effects on the observed CV. The first condition, referred to as A ₁, is that the latter occurred fault has at least one inconsistent CV element which is consistent in the previous CV. The second condition, referred to as A ₂, is that the latter occurred fault has a CV which is masked by the previous one. These two conditions can be expressed as

A_{1} : c_{1}^{e} - c_{1}^{l} \geq 0 \land c_{2}^{e} - c_{2}^{l} \geq 0 \land \dots \land c_{m}^{e} - c_{m}^{l} \geq 0 = 0

(4)

\begin{matrix} A_{2} : c_{1}^{e} - c_{1}^{l} \geq 0 \land c_{2}^{e} - c_{2}^{l} \geq 0 \land \dots \land c_{m}^{e} - c_{m}^{l} \geq 0 = 1 \\ with C^{l} \neq 0 \end{matrix}

(5)

where $c_{i}^{e}, i = 1, 2, \dots, m$ is the ith element of the earlier occurred fault CV (i.e. $C^{e}$ ), and $c_{i}^{l}, i = 1, 2, \dots, m$ is the ith element of the latter occurred fault CV (i.e. $C^{l}$ ).

It is worth noting that for the activation condition A ₂ in equation (5), the system cannot obtain the fault information $C^{l} \neq 0$ because the latter occurred fault is masked by the previous occurred fault. As a result, more information is required to obtain condition A ₂ for fault identification activation. Let us look at equation (1) more closely if a fault occurs in parametric component in $a$ or nonparametric component in $β$ , the corresponding AGARR equation which contains the fault component should exceed the threshold. Thus, its CV element is equal to one. Once a latter fault occurs which is masked by the previous fault, the accumulation effect of both faults will lead to a sudden increase in residual at the instance of latter fault occurrence. In order to capture this phenomenon, a SLR is defined as the derivative of AGARR equation which contains the latter fault as follows

SL R_{j} = \frac{d}{dt} AGAR R_{j}

(6)

where $AGAR R_{j}$ is the AGARR equation which contains the latter fault.

According to equations (5) and (6), the second activation condition can be reformulated as

\begin{matrix} A_{2} : c_{1}^{e} - c_{1}^{l} \geq 0 \land c_{2}^{e} - c_{2}^{l} \geq 0 \land \dots \land c_{m}^{e} - c_{m}^{l} \geq 0 = 1 \\ and | SL R_{j} | \geq ε \end{matrix}

(7)

where $ε$ is the threshold for SLR.

The complete structure of the proposed sequential fault diagnosis method is demonstrated in Figure 1. The initial inconsistency is detected by residuals of AGARRs. Since the mode information is unavailable, this inconsistency (i.e. $C \neq 0$ ) may be due to a component fault or an unknown mode change; thus a suspected mode change is examined first. The MCSM and CV are compared to find the suspected modes. After that, mode identification is enabled. If a new mode is identified, the system is considered as failure. If a new mode is not found, the component fault isolation is enabled to find the fault candidates set based on the observed CV.

Figure 1.

Schematic description of the proposed sequential fault diagnosis method.

The first SFI, denoted as SFI₁, is activated to identify the true faults if fault mode identification is unsuccessful and fault candidate set is obtained. After that, the monitoring process is carried on. If the identification activation condition, that is, A ₁ in equation (4) or A ₂ in equation (7), is satisfied, SFI₂ will be enabled based on the identification results from SFI₁. For each SFI except SFI₁, the required unknown parameters are reduced because it is based on the identification results from the previous SFI (i.e. SFI_i based on SFI_i−1 for $i \geq 2$ ). In addition, the fault information could be obtained as early as possible for further maintenance purpose using the proposed SFI scheme.

HS and CHS

HS is a newly developed meta-heuristic optimization method which is invented by Geem¹¹ and Geem et al.¹² It draws its inspiration from the improvisation process of musicians seeking a better state of harmony. The effort to find musically pleasing harmony determined by an esthetic standard is analogous to the effort to find the global optimal solution determined by an objective function in an optimization process. Many published works demonstrate that HS converges fast and is simple in implementation and requires only a few control parameters.^13–18 For most optimization problems, the design variables are continuous. However, many practical engineering problems require considering the design variables as integer or binary values.¹⁹ The presence of binary variables along with continuous variables adds to the complexity of the optimization problem. To the best of our knowledge, the majority of previous efforts on HS applications are focused on solving problems in discrete or continuous space, and few efforts are dedicated to problems including both binary and continuous variables.

This article presents a CHS approach for solving the fault identification problem formulated as a nonlinear optimization problem which contains binary variables as well as continuous variables. In the CHS, design variables are represented as a combination of binary strings for $η_{F}$ in equation (3) and continuous strings for $F_{ab}$ , $F_{int}$ , and $T_{F}$ in equation (3). Therefore, the number of design variables for fault identification is quadruple the number of fault candidates after component fault isolation. The purpose of fault identification is to find the true faults with relevant information. Thus, the identification problem can be reformulated as an optimization problem with the following objective function to be minimized

f_{obj} = (\sum_{l = 1}^{m} \sum_{n = 1}^{q} | G_{l}^{n} |)

(8)

where $G_{l}$ is the lth AGARR equation and $n$ is the discrete sampling index.

The first step of CHS is the initialization of the harmony memory as

HM = [\begin{matrix} x_{1}^{1} \\ x_{1}^{2} \\ ⋮ \\ x_{1}^{M} \end{matrix} \begin{matrix} \dots \\ \dots \\ \dots \\ \dots \end{matrix} \begin{matrix} x_{j}^{1} \\ x_{j}^{2} \\ ⋮ \\ x_{j}^{M} \end{matrix} \begin{matrix} \dots \\ \dots \\ \dots \\ \dots \end{matrix} \begin{matrix} x_{N}^{1} \\ x_{N}^{2} \\ ⋮ \\ x_{N}^{M} \end{matrix} | \begin{matrix} f_{obj} (X^{1}) \\ f_{obj} (X^{2}) \\ ⋮ \\ f_{obj} (X^{M}) \end{matrix}]

(9)

where $X^{i} = [x_{1}^{i}, \dots, x_{j}^{i}, \dots, x_{N}^{i}], i = 1, 2, \dots, M$ is the harmony vector where $x_{k}^{i}, k = 1, 2, \dots, j$ are continuous design variables and $x_{k}^{i}, k = j + 1, \dots, N$ are binary design variables; $N$ is the number of decision variables; and $M$ is the harmony memory size (HMS).

It is not hard to find that $j = 3 (N - j)$ in equation (9) because the number of continuous design variables triple the number of binary design variables according to equation (3). During initialization, each continuous design variable is randomly initialized within the feasible solution space, and each binary design variable is initialized with random binary number. After initialization, a new harmony vector $X' = [x_{1}', x_{2}', \dots, x_{j}', x_{j + 1}', \dots, x_{N}']$ is improvised by the following three mechanisms: memory consideration, random selection, and pitch adjustment.¹⁶

The harmony memory considering rate (HMCR) indicates whether the element of new candidates is generated from the harmony memory or random selection. In other words, the HMCR is the rate of choosing one value from the harmony memory. (1 − HMCR) is the rate of randomly selecting one value from the feasible solution space for continuous design variable and is the rate of randomly selecting to be 0 or 1 for binary design variable. Therefore, for continuous design variable

x'_{k} \leftarrow {\begin{matrix} {x'}_{k} \in {x_{k}^{1}, x_{k}^{2}, \dots, x_{k}^{M}} \\ {x'}_{k} \in [x_{k}^{L}, x_{k}^{U}] \end{matrix} \begin{matrix} w . p . HMCR \\ w . p . (1 - HMCR) \end{matrix}

(10)

where $x_{k}^{L}$ and $x_{k}^{U}$ are the lower and upper bounds for each continuous decision variable, respectively.

As for binary design variable, memory consideration and random selection can be expressed as¹⁹

x'_{k} \leftarrow {\begin{matrix} {x'}_{k} \in {x_{k}^{1}, x_{k}^{2}, \dots, x_{k}^{M}} \\ R \end{matrix} \begin{matrix} w . p . HMCR \\ w . p . (1 - HMCR) \end{matrix}

(11)

with

R = {\begin{matrix} 0 \\ 1 \end{matrix} \begin{matrix} if r_{k} < 0.5 \\ else \end{matrix}

(12)

where $r_{k}$ is a random number between 0 and 1.

Each decision variable obtained from the memory consideration is examined to determine whether it should be pitch adjusted with pitch adjusting rate (PAR). Thus, for continuous design variable, pitch adjustment can be formulated as¹⁷

x'_{k} \leftarrow {\begin{matrix} {x'}_{k} \pm U (- 1, 1) \times bw \\ {x'}_{k} \end{matrix} \begin{matrix} w . p . HMCR \times PAR \\ w . p . HMCR \times (1 - PAR) \end{matrix}

(13)

where $bw$ is bandwidth for pitch adjustment of continuous design variable and $U (- 1, 1)$ is a uniform distribution between 0 and 1.

As for binary design variable, the following pitch adjustment is used¹⁹

x'_{k} \leftarrow {\begin{matrix} {x'}_{k} = x_{k}^{best} \\ {x'}_{k} \end{matrix} \begin{matrix} w . p . HMCR \times PAR \\ w . p . HMCR \times (1 - PAR) \end{matrix}

(14)

where $x_{k}^{best}$ is the corresponding element of global optimal harmony vector in harmony memory.

During updating stage, if the new harmony vector is better than the worst harmony vector in harmony memory in terms of the objective function defined in equation (8), the new harmony is included in the harmony memory, and the existing worst harmony is excluded from the harmony memory. Finally, if the stopping criterion or maximum number of improvisation (NI) is satisfied, the computation is terminated. Note that for the suspected mode identification, only binary part of the CHS is utilized.

Experimental study

Mobile robot system description

In order to verify the proposed sequential fault diagnosis approach, a test on a mobile robot system is conducted. The robot system, as shown in Figure 2, consists of direct current (DC) motor, gear, belt, oil pump, hydraulic actuator, steering mechanism, and tires. The steering power is provided from the hydraulic actuator which generates pressure difference to push the oil piston inside the actuator. The oil piston is connected to the steering mechanism. The pump is driven by the DC motor through the belt and converts mechanical speed to proportional oil flow. Four sensors, that is, two pressure sensors, an absolute encoder, and an incremental encoder, are installed to measure various signals in the system. The absolute encoder is utilized to measure the steering angle, and the incremental encoder is adopted to measure the DC motor speed on the reducing output shaft.

Figure 2.

Test bed for mobile robot system.

The DHBG of the robot system is shown in Figure 3. A DHBG is a HBG that is assigned with proper set of causalities, such that the causality of every active BG component is valid and consistent at all modes. In this figure, source $M_{Sf} : u_{in}$ represents the input of motor driver. $0_{c 1}$ models the controlled junction which represents a burnt driver fault or a burnt motor fault. When a fault occurs in the driver or the motor, the state variable $a_{1}$ of junction $0_{c 1}$ is 0. Similarly, controlled junction $0_{c 2}$ is used to model the broken belt fault and $a_{2} = 0$ when the belt is broken. The steering mechanism is modeled by $F_{1} (θ_{2})$ and $F_{2} (θ_{2})$ which are derived from Ackerman’s geometry. In order to model the internal leakage of the piston, resistance element $R : R_{4}$ is used to represent the internal resistance to oil flow. When the piston is normal, $R_{4} \to \infty$ , and $R_{4} << \infty$ under leakage condition. The dynamics of the leakage can be described as

f_{19} = \frac{sign (e_{19}) \sqrt{| e_{19} |}}{R_{4}} = \frac{sign (Δ p) \sqrt{| Δ p |}}{R_{4}}

(15)

where $Δ p = p_{1} - p_{2}$ is the pressure difference of the two pressure sensors and $sign (\cdot)$ represents the sign function.

Figure 3.

DHBG model of the robot system.

$Df : {\overset{\cdot}{θ}}_{1}$ is the velocity ${\overset{\cdot}{θ}}_{1}$ at the gear output shaft, and $Df : {\overset{\cdot}{θ}}_{2}$ measures the steering angle $θ_{2}$ . Other parameters with their physical meanings are summarized in Table 1.

Table 1.

Physical parameters of the mobile robot system.

Parameter	Physical meaning	Nominal value	Unit
$k_{1}$	Voltage-to-current constant of motor driver	3	A/V
$k_{2}$	Current-to-torque ratio of DC motor	5.627e⁻⁴	N m/A
$J_{1}$	DC motor inertia	3.727e⁻⁵	kg m²
$k f_{2}$	Viscous friction coefficient in $R_{2}$ of DC motor	1.017e⁻⁵	N m/(rad/s)
$F u_{2}$	Coulomb friction coefficient in $R_{2}$ of DC motor	6.20e⁻²	N m
$k_{3}$	Reducing gear ratio	4/70	–
$k_{4}$	Belt ratio	1/2	–
$J_{2}$	Pump rotor inertia	5.627e⁻⁴	kg m²
$k f_{3}$	Viscous friction coefficient in $R_{3}$ of pump rotor	3.161e⁻¹	N m/(rad/s)
$F u_{3}$	Coulomb friction coefficient in $R_{3}$ of pump rotor	3.55e⁻²	N m
$k_{5}$	Ratio between pump angular velocity and oil flow	1.024e⁻⁴	m³
$A$	Effective cross area of the piston	4.46e⁻⁴	m²
$J_{3}$	Tire inertia	9.78e⁻²	kg m²
$k f_{5}$	Viscous friction coefficient in $R_{5}$ of tire	33.28	N m/(rad/s)
$F u_{5}$	Coulomb friction coefficient in $R_{5}$ of tire	22.11	N m

DC: direct current.

From the DHBG of the system, three independent AGARR equations are derived. First, consider the constitutive relation of junction $1_{3}$ attached with flow sensor $Df : {\overset{\cdot}{θ}}_{1}$

e_{11} = e_{10} - e_{12} = 0

(16)

Tracking back the causal paths on DHBG, the unknown variables $e_{10}$ and $e_{12}$ can be expressed as

\begin{matrix} e_{10} = k_{3}^{- 1} e_{9} = k_{3}^{- 1} (e_{6} - e_{7} - e_{8}) \\ = k_{3}^{- 1} (k_{2} f_{5} - J_{1} {\overset{\cdot}{f}}_{7} - k_{f_{2}} f_{8} - F u_{2} sign (f_{8})) \\ = k_{3}^{- 1} (a_{1} k_{1} k_{2} u_{in} - J_{1} k_{3}^{- 1} {\overset{\cdot\cdot}{θ}}_{1} - k_{f_{2}} k_{3}^{- 1} {\overset{\cdot}{θ}}_{1} - F u_{2} sign ({\overset{\cdot}{θ}}_{1})) \end{matrix}

(17)

\begin{matrix} e_{12} = a_{2} e_{13} = a_{2} k_{4} e_{14} = a_{2} k_{4} (e_{15} + e_{16} + e_{17}) \\ = a_{2} k_{4} (J_{2} {\overset{\cdot}{f}}_{15} + k_{f_{3}} f_{16} + F u_{3} sign (f_{16}) + k_{5} e_{18}) \\ = a_{2} k_{4} (J_{2} k_{4} {\overset{\cdot\cdot}{θ}}_{1} + k_{f_{3}} k_{4} {\overset{\cdot}{θ}}_{1} + F u_{3} sign ({\overset{\cdot}{θ}}_{1}) + k_{5} Δ p) \end{matrix}

(18)

Substituting equations (17) and (18) into equation (16) yields

A G A R R_{1} = k_{3}^{- 1} (a_{1} k_{1} k_{2} u_{i n} - J_{1} k_{3}^{- 1} \frac{d^{2}}{d t^{2}} (\frac{θ_{1}}{β_{θ_{1}}}) - k_{f_{2}} k_{3}^{- 1} \frac{d}{d t} (\frac{θ_{1}}{β_{θ_{1}}}) - F u_{2} s i g n (\frac{d}{d t} (\frac{θ_{1}}{β_{θ_{1}}}))) - a_{2} k_{4} (J_{2} k_{4} \frac{d^{2}}{d t^{2}} (\frac{θ_{1}}{β_{θ_{1}}}) + k_{f_{3}} k_{4} \frac{d}{d t} (\frac{θ_{1}}{β_{θ_{1}}}) + F u_{3} s i g n (\frac{d}{d t} (\frac{θ_{1}}{β_{θ_{1}}})) + k_{5} \frac{Δ p}{β_{Δ p}}) = 0

(19)

Next, consider junction 0₁ with sensor $De : Δ p$ , the constitutive relation of this junction can be formulated as

f_{20} = f_{18} - f_{19} - f_{21} = 0

(20)

The flow variables $f_{18}$ and $f_{21}$ can be represented as

f_{18} = k_{5} f_{17} = a_{2} k_{4} k_{5} {\overset{\cdot}{θ}}_{1}

(21)

f_{21} = A f_{22} = {AF}_{1}^{- 1} (θ_{2}) {\overset{\cdot}{θ}}_{2}

(22)

Combining equations (20)–(22) and equation (15), AGARR₂ can be expressed as

\begin{matrix} AGAR R_{2} = a_{2} k_{4} k_{5} \frac{d}{dt} (\frac{θ_{1}}{β_{θ_{1}}}) \\ - \frac{sign (\frac{Δ p}{β_{Δ p}}) \sqrt{| \frac{Δ p}{β_{Δ p}} |}}{R_{4}} - {AF}_{1}^{- 1} (\frac{θ_{2}}{β_{θ_{2}}}) \frac{d}{dt} (\frac{θ_{2}}{β_{θ_{2}}}) = 0 \end{matrix}

(23)

Finally, consider junction 1₅ attached with flow sensor $Df : {\overset{\cdot}{θ}}_{2}$

e_{24} = e_{23} - e_{25} = 0

(24)

Unknown variables $e_{23}$ and $e_{25}$ can be solved as

e_{23} = F_{1}^{- 1} (θ_{2}) e_{22} = {AF}_{1}^{- 1} (θ_{2}) Δ p

(25)

\begin{matrix} e_{25} = F_{2} (θ_{2}) e_{26} = F_{2} (θ_{2}) (e_{27} + e_{28}) \\ = F_{2} (θ_{2}) (J_{3} {\overset{\cdot}{f}}_{28} + k_{f_{5}} f_{27} + F u_{5} sign (f_{27})) \\ = F_{2} (θ_{2}) (J_{3} \frac{d}{dt} (F_{2} (θ_{2}) {\overset{\cdot}{θ}}_{2}) + k_{f_{5}} F_{2} (θ_{2}) {\overset{\cdot}{θ}}_{2} + F u_{5} sign (F_{2} (θ_{2}) {\overset{\cdot}{θ}}_{2})) \end{matrix}

(26)

The third AGARR can be obtained according to equations (24)–(26)

\begin{matrix} AGAR R_{3} = {AF}_{1}^{- 1} (\frac{θ_{2}}{β_{θ_{2}}}) \frac{Δ p}{β_{Δ p}} - J_{3} F_{2} (\frac{θ_{2}}{β_{θ_{2}}}) \frac{d}{dt} \\ [F_{2} (\frac{θ_{2}}{β_{θ_{2}}}) \frac{d}{dt} (\frac{θ_{2}}{β_{θ_{2}}})] - k_{f_{5}} F_{2}^{2} (\frac{θ_{2}}{β_{θ_{2}}}) \frac{d}{dt} (\frac{θ_{2}}{β_{θ_{2}}}) \\ - F u_{5} F_{2} (\frac{θ_{2}}{β_{θ_{2}}}) sign [F_{2} (\frac{θ_{2}}{β_{θ_{2}}}) \frac{d}{dt} (\frac{θ_{2}}{β_{θ_{2}}})] = 0 \end{matrix}

(27)

Table 2 is the MD-FSM which represents the cause–effect relations between component faults and residuals. Table 3 is the MCSM which represents cause–effect relations between mode changes and residuals. It is worth noting that more sophisticated model can be developed for the mobile robot system if effects of nonlinear phenomena such as saturation or friction are considered. At the same time, the model complexity is increased where many physical parameters related to nonlinear phenomena (e.g. the flat tire fault is represented by the change of friction coefficients) could be possible fault candidates. Therefore, fault isolation will lead to many fault candidates which adds to the complexity of fault identification. In order to achieve the tradeoff between the model accuracy and fault diagnosis complexity, the model with linear friction is preferred as long as the fault diagnosis performance is acceptable.

Table 2.

MD-FSM at mode [a ₁ a ₂] = [1 1].

	AGARR₁	AGARR₂	AGARR₃	$D_{b}$
$R_{4}$	0	1	0	1
$R_{5}$	0	0	1	1
$β_{θ_{1}}$	1	1	0	1
$β_{θ_{2}}$	0	1	1	1
$β_{Δ p}$	1	1	1	1

MD-FSM: mode-dependent fault signature matrix; AGARR: augmented global analytical redundancy relation.

Table 3.

MCSM of the robot system.

	AGARR₁	AGARR₂	AGARR₃	$D_{b}$	$I_{b}$
$a_{1}$	1	0	0	1	1
$a_{2}$	1	1	0	1	1

MCSM: mode change signature matrix; AGARR: augmented global analytical redundancy relation.

Experimental results and discussions

The robot system is an open-loop system. The motion of the DC motor of the system is generated by a current-modulated pulse width modulation (PWM) voltage from a current amplifier. The command input signal is defined as follows (in V)

u_{in} = 0.17 (5 \sin 2 t + 0.5 \sin t + 0.2 \sin 0.5 t + \sin 2.5 t)

(28)

During physical parameter identification of the mobile robot system under normal condition, a sine signal $u_{in} = 0.7 \sin t$ is used as the system input. In order to verify the fault diagnosis performance of the proposed method, another signal in equation (28) is used to test the algorithm in different working conditions.

Note that in closed-loop systems, controllers try to hide fault effects, and thus, it is usually impossible to find faults by simply observing the system outputs.²⁰ For example, a tank system with a controller to maintain a constant level, if a tank leakage fault occurs, the controller tries to hide fault effects by increasing its output to maintain a constant tank level. If one observes the system output, that is, the tank level, the fault cannot be detected. However, in DHBG-based fault diagnosis method using the AGARRs, the controller cannot hide the fault effects because AGARR residuals remain sensitive to fault hidden by the action of the closed-loop controller. Let us still consider the closed-loop tank system, if a leakage fault occurs, the AGARR residual will exceed the threshold because the residual of AGARR is persistently sensitive to the leakage fault. Thus, there is no difference between fault detection in closed-loop systems and the one in open-loop systems when DHBG-based methods are used.

Three fault scenarios are considered in the experiment. In the first fault scenario, an abrupt internal leakage fault at t = 32 s with faulty value ${\bar{R}}_{4_ab} = 2.5221 e^{7} k g^{- 1 / 2} m^{- 1 / 2}$ and an intermittent sensor fault in $Df : θ_{1}$ at t = 42 s with faulty values ${\bar{β}}_{θ_{1_int}} = 0.5$ and ${\bar{T}}_{β_{θ_{1}}} = 3 s$ are introduced in the system. When the cylinder is normal, oil cannot flow from one chamber to the other. The piston efficiency is decreased if the internal flow between the chambers exists. The injection of internal leakage fault in cylinder is realized by a manual valve as shown in Figure 4. The manual valve is used to control the oil flow in the bypass tube which is mounted to allow free oil flow from one side to the other side of the cylinder. When the valve is completely closed, the system is fault free. An open valve indicates abnormal condition where the opening level of the valve determines the fault severity.

Figure 4.

Injection of internal leakage fault by a manual valve.

Figure 5 illustrates the residual responses under the abrupt internal leakage fault and intermittent sensor fault where dashed lines are the thresholds. At 32 s, a CV = [0 1 0] is observed from residuals. It is known that the CV is not caused by fault mode from Table 3. As a result, the CV is caused by a component fault, that is, an internal leakage fault in $R_{4}$ , from Table 2.

Figure 5.

Residual responses under asynchronous faults: an abrupt internal leakage fault and an intermittent sensor fault in $Df : θ_{1}$ .

Since the knowledge of the fault type is unavailable in advance, the CHS algorithm is enabled for SFI₁ to identify the fault type with corresponding information. The parameters associated with CHS are chosen as: HMS = 100, PAR = 0.05, HMCR = 0.8, and NI = 200. Figures 6 and 7 show the identification results of SFI₁. From these figures, it is known that the internal leakage fault is abrupt because $η_{R_{4}} = 1$ . The identified fault magnitude is $R_{4_ab} = 2.5277 e^{7} k g^{- 1 / 2} m^{- 1 / 2}$ which is close to its true value $2.5221 e^{7} k g^{- 1 / 2} m^{- 1 / 2}$ . There is no hide effect since CV = [0 1 0] is unique, so condition A ₂ is not satisfied. The monitoring process is carried on until t = 42 s, a new CV = [1 1 0] which satisfies A ₁ is observed, and then SFI₂ is activated based on SFI₁ with fault candidate set $δ = [β_{θ_{1}}, R_{4}]$ . The identification results of SFI₂ are illustrated in Figures 8 –10. The fault type for $β_{θ_{1}}$ is intermittent because $η_{β_{θ_{1}}} = 0$ . The identified fault amplitude and period are $β_{θ_{1_int}} = 0.508$ and $T_{β_{θ_{1}}} = 3.05 s$ which match their designed values. The results successfully identify the intermittent sensor fault using the previous identified abrupt internal leakage fault.

Figure 6.

Identification for $R_{4_ab}$ .

Figure 7.

Identification for $η_{R_{4}}$ .

Figure 8.

Identification for $β_{θ_{1_int}}$ .

Figure 9.

Identification for $T_{β_{θ_{1}}}$ .

Figure 10.

Identification for $η_{β_{θ_{1}}}$ .

For the second fault scenario, an abrupt internal leakage fault with faulty value ${\bar{R}}_{4_ab} = 2.5221 e^{7} k g^{- 1 / 2} m^{- 1 / 2}$ at t = 32 s and an abrupt sensor fault in $Df : θ_{1}$ with designed value ${\bar{β}}_{θ_{1_ab}} = 0.5$ at t = 15 s are considered. Figure 11 shows the residual responses under the second fault scenario. It is observed that at t = 15 s, the CV = [1 1 0]. From Table 3, fault mode of broken belt is suspected, and the fault mode identification is enabled. The fault mode identification result is demonstrated in Figure 12 in which no new mode is identified. Thus, component fault isolation leads to fault candidate set $δ = [β_{θ_{1}}, R_{4}]$ . In order to obtain the true fault information, SFI₁ is activated with identification results as shown in Figures 13 –16. Note that for the fault in $R_{4}$ , the identified $η_{R_{4}} = 1$ and $R_{4_ab} = 8.9108 e^{9} k g^{- 1 / 2} m^{- 1 / 2}$ which is much bigger than the abrupt faulty value of internal leakage fault. Therefore, $R_{4}$ is considered as normal. The identified $η_{β_{θ_{1}}} = 1$ and $β_{θ_{1_ab}} = 0.5102$ which match designed values.

Figure 11.

Residual responses under asynchronous faults: an abrupt internal leakage fault and an abrupt sensor fault in $Df : θ_{1}$ .

Figure 12.

Identification for fault mode.

Figure 13.

Identification for $β_{θ_{1_ab}}$ .

Figure 14.

Identification for $η_{β_{θ_{1}}}$ .

Figure 15.

Identification for $R_{4_ab}$ .

Figure 16.

Identification for $η_{R_{4}}$ .

Since the observed CV = [1 1 0] has the ability to hide internal leakage fault with CV = [0 1 0] if sensor fault in $Df : θ_{1}$ occurs before the internal leakage fault, condition A ₂ could be satisfied if $| SL R_{2} | \geq ε$ is observed after 15 s. Around t = 32 s, an inconsistency in $SL R_{2}$ is detected as shown in Figure 17 which indicates a latter fault occurrence in $R_{4}$ ; thus SFI₂ is activated based on SFI₁, and the identified results are shown in Figures 18 and 19. The method successfully identifies the abrupt internal leakage fault with $η_{R_{4}} = 1$ and $R_{4_ab} = 2.6129 e^{7} k g^{- 1 / 2} m^{- 1 / 2}$ using the previous identified abrupt sensor fault from SFI₁.

Figure 17.

Response of SLR₂.

Figure 18.

Identification for $R_{4_ab}$ .

Figure 19.

Identification for $η_{R_{4}}$ .

In the third fault scenario, an actuator fault (i.e. a burnt DC motor fault) is considered. This fault is a fault mode where no power is delivered to the system upon the fault occurrence. The burnt motor fault is injected at t = 31 s by disabling the motor. Figure 20 shows the residual responses under the fault. Around t = 31 s, a CV of [1 0 0] is detected. Since the CV is unique according to Tables 1 and 2, it is concluded that this fault signature is caused by the burnt motor fault. Therefore, the burnt motor fault is detected and isolated.

Figure 20.

Residual responses under fault mode: burnt DC motor.

In order to demonstrate the superiority of the proposed CHS-based sequential fault identification method, comparison study with genetic algorithm (GA), adaptive genetic algorithm (AGA), and hybrid differential evolution (HDE) is conducted. All methods are tested on the same experimental data which are collected for SFI₁ of the first fault scenario. In GA, all possible solutions are encoded into binary chromosomes where $R_{4_ab}$ , $R_{4_int}$ , and $T_{R_{4}}$ use 3 bits each, and the binary variable $η_{R_{4}}$ uses 1 bit. During the search process, decoding is not required for the chromosomes representing the binary variable $η_{R_{4}}$ . To apply GA, individuals (i.e. possible solutions) strive for survival based on the Darwinian principle of the survival of the fittest. There are four major steps in GA: initialization, evaluation, selection and reproduction, and crossover and mutation. During initialization, a set of binary chromosomes which consist of all the aforementioned parameters for encoding form the individuals in the GA. In evaluation, each individual is evaluated by the fitness function in equation (8).

For selection and reproduction, individuals with the highest fitness values are kept in the next generation, while those with the lowest fitness values are abandoned. The high-quality individual has a greater chance to recombine with other individuals to reproduce offspring according to the genetic operators of crossover and mutation. The major structure of AGA is same with the one of GA. The only difference between them is that AGA is able to adaptively tune the crossover and mutation rates based on the performance of the current genetic operators.²¹ It can adjust the balance between the exploration and exploitation of the solution space.

In HDE, real-valued differential evolution (RDE) is adopted to find $R_{4_ab}$ , $R_{4_int}$ , and $T_{R_{4}}$ , and binary-valued differential evolution (BDE) is utilized to search the binary variable $η_{R_{4}}$ .⁹ Using strategy DE/best/1/bin in RDE, the mutation operator produces a new mutant individual by adding the weighted difference between two randomly chosen individuals to the best individual in population at current generation. In crossover, a trial vector is generated by mixing the parameters of the mutant individual with those of the target individual. Selection operates by comparing the individual fitness to generate the new population of next generation. Based on the selection process, all the individuals of the next generation are as good as or better than their counterparts in the current generation. In BDE, crossover operators can be directly used. Mutation is the only operator which should be modified where Boolean algebra “AND,”“OR,” and “XOR” are used to replace their counterparts in RDE.

To perform a fair comparison, the population size and the maximum iteration of all approaches are set equal to 100 and 200, respectively. For each algorithm, 50 independent runs are carried out. Table 4 summarizes the minimum objective function $f_{obj}^{\min}$ , the maximum objective function $f_{obj}^{\max}$ , and the average objective function $f_{obj}^{mean}$ over the 50 runs. In this table, the standard deviation values $σ$ of different approaches are also provided. It is obvious that CHS outperforms GA and AGA in terms of final solution and standard deviation. The performances of CHS and HDE are comparable, but CHS is more efficient because it is simple in implementation.

Table 4.

Performance comparison of different methods.

	$f_{obj}^{\min}$	$f_{obj}^{\max}$	$f_{obj}^{mean}$	$σ$
CHS	1.36e⁻³	1.72e⁻³	1.47e⁻³	5.26e⁻⁵
GA	1.85e⁻³	2.53e⁻³	2.03e⁻³	9.32e⁻⁵
AGA	1.63e⁻³	2.36e⁻³	1.77e⁻³	7.71e⁻⁵
HDE	1.46e⁻³	1.89e⁻³	1.55e⁻³	5.58e⁻⁵

CHS: composite harmony search; GA: genetic algorithm; AGA: adaptive genetic algorithm; HDE: hybrid differential evolution.

Conclusion

In this article, a new sequential fault diagnosis method to handle asynchronous distinct faults using DHBG and CHS is proposed. The monitored system has various fault types including fault mode, abrupt fault, and intermittent fault. The fault-type information is unknown in advance, and the SFI algorithm is reactivated during monitoring process based on two activation conditions. The proposed method is able to handle continuous variables and binary variables simultaneously. The developed method can identify the fault information as early as possible for future maintenance purpose. Three different fault scenarios have been experimentally tested on a mobile robot system to illustrate the effectiveness of the proposed method.

Footnotes

Academic Editor: Zhijun Li

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Huangshan Young Scholar Fund of Hefei University of Technology under grant no. 407037159.

References

Djeziri

Merzouki

Ould-Bouamama

. Robust monitoring of an electric vehicle with structured and unstructured uncertainties. IEEE T Veh Technol 2009; 58: 4710–4719.

Pan

. Machine health condition prediction via online dynamic fuzzy neural networks. Eng Appl Artif Intel 2014; 35: 105–113.

Loureiro

Merzouki

Ould-Bouamama

. Bond graph model based on structural diagnosability and recoverability analysis: application to intelligent autonomous vehicles. IEEE T Veh Technol 2012; 61: 986–997.

Calderaro

Hadjicostis

Piccolo

. Failure identification in smart grids based on petri net modeling. IEEE T Ind Electron 2011; 58: 4613–4623.

Narasimhan

Biswas

. Model-based diagnosis of hybrid systems. IEEE T Syst Man Cy A 2007; 37: 348–361.

Low

Wang

Arogeti

. Quantitative hybrid bond graph-based fault detection and isolation. IEEE T Autom Sci Eng 2010; 7: 558–569.

Arogeti

Wang

Low

. Mode identification of hybrid systems in the presence of fault. IEEE T Ind Electron 2010; 57: 1452–1467.

Arogeti

Wang

Low

. Fault detection isolation and estimation in a vehicle steering system. IEEE T Ind Electron 2012; 59: 4810–4820.

Wang

. Model based prognosis for hybrid systems with mode-dependent degradation behaviors. IEEE T Ind Electron 2014; 61: 546–554.

10.

Wen

Soh

. Observer-based stabilization of switching linear systems. Automatica 2003; 39: 517–524.

11.

Geem

. Music-inspired harmony search algorithm. Berlin: Springer, 2008.

12.

Geem

Kim

Loganathan

. A new heuristic optimization algorithm: harmony search. Simulation 2001; 76: 60–68.

13.

Manjarres

Landa-Torres

Gil-Lopez

. A survey on applications of the harmony search algorithm. Eng Appl Artif Intel 2013; 7: 1818–1831.

14.

Sarvari

Zamanifar

. Improvement of harmony search algorithm by using statistical analysis. Artif Intell Rev 2012; 37: 181–215.

15.

Omran

MGH

Mahdavi

. Global-best harmony search. Appl Math Comput 2008; 198: 643–656.

16.

Wang

Yuan

Wang

. Harmony search algorithm-based fuzzy-PID controller for electronic throttle valve. Neural Comput Appl 2013; 22: 329–336.

17.

Zou

Gao

. A novel global harmony search algorithm for reliability problems. Comput Ind Eng 2010; 58: 307–316.

18.

Sinsuphan

Leeton

Kulworawanichpong

. Optimal power flow solution using improved harmony search method. Appl Soft Comput 2013; 13: 2364–2374.

19.

Wang

Yang

. An improved adaptive binary harmony search algorithm. Inform Sciences 2013; 232: 58–87.

20.

Samantaray

Ould-Bouamama

. Model-based process supervision: a bond graph approach. London: Springer, 2008.

21.

Mak

Wong

Wang

. An adaptive genetic algorithm for manufacturing cell formation. Int J Adv Manuf Tech 2000; 16: 491–497.