Jha Mayank Shekhar

Particle Filter Based Hybrid Prognostics for Health Monitoring of Uncertain Systems in Bond Graph Framework

Keywords: Prognostics, Bond Graph, Intervals, Particle Filter, Remaining Useful Life, Robust Fault Detection

The paper"s main objective is to address the problem of health monitoring of system parameters in Bond Graph (BG) modeling framework, by exploiting its structural and causal properties. The system in feedback control loop is considered globally uncertain and parametric uncertainty is modeled in interval form. The system parameter is undergoing degradation (prognostic candidate) and its degradation model is assumed to be known a priori. The detection of degradation commencement is done in a passive manner which involves interval valued robust adaptive thresholds over the nominal part of the uncertain BG-derived Interval Valued Analytical Redundancy Relations (I-ARRs). The latter forms an efficient diagnostic module. The prognostics problem is cast as joint state-parameter estimation problem, a hybrid prognostic approach, wherein the fault model is constructed by considering the statistical degradation model of the system parameter (prognostic candidate). The observation equation is constructed from nominal part of the I-ARR. Using Particle Filter (PF) algorithms; the estimation of state of health (state of prognostic candidate) and associated hidden time-varying degradation progression parameters is achieved in probabilistic terms. A simplified variance adaptation scheme is proposed. Associated uncertainties which arise out of noisy measurements, parametric degradation process, environmental conditions etc. are effectively managed by PF. This allows the production of effective predictions of the remaining useful life of the prognostic candidate with suitable confidence bounds. The effectiveness of the novel methodology is demonstrated through simulations and experiments on a mechatronic system.

Introduction

Health Monitoring aims at ensuring system safety, reliability and efficient functionality and deals with fault detection and prediction of the Remaining Useful Life (RUL) of the system in a holistic way. While the former two is mainly dealt by using a diagnostic module, the latter is performed by a prognostic module. The primary focus lies in scheduling the maintenance actions according to progression of the system to a time where it may be considered beyond the limits of certified functionalities [START_REF] Jardine | A review on machinery diagnostics and prognostics implementing condition-based maintenance[END_REF]. Such a time-horizon of interest is termed as the End of Life (EOL) and the time remaining until that point is called RUL of the system [START_REF] Lee | Prognostics and health management design for rotary machinery systems-Reviews, methodology and applications[END_REF][START_REF] Sikorska | Prognostic modelling options for remaining useful life estimation by industry[END_REF]. Prognostics are focused on the study of fault (or damage) evolution and prediction of the RUL of the system/component. Accurate prediction of EOL/RUL enables efficient and optimal planning of the future maintenance actions, and renders the capability of assessing reliability of the system [START_REF] Heng | Rotating machinery prognostics: State of the art, challenges and opportunities[END_REF]. This leads to system/component"s life extension by modification of the system demand, operating conditions, workload etc. [START_REF] Sun | Application of a state space modeling technique to system prognostics based on a health index for condition-based maintenance[END_REF].

The failures of most systems can be attributed to the degradation of a given component, subsystem or material with time, environmental and operational conditions etc. Such system components/sub-systems can be identified as the poten-tial prognostic candidates through Failure Modes, Mechanisms and Effect Analysis or through other ways [START_REF] Louit | Condition-based spares ordering for critical components[END_REF]. The underlying physical degradation is usually captured by Degradation Model(s) (DM) that can be obtained based upon physics of degradation or statistical (experimental) modeling approach as described in Gebraeel et al. [START_REF] Gebraeel | Prognostic degradation models for computing and updating residual life distributions in a time-varying environment, Reliability[END_REF] and Guo et al. [START_REF] Guo | Practical Approaches for Reliability Evaluation Using Degradation Data[END_REF]. In cases where physics of degradation is not available or reliable, the respective DM can be obtained statistically by finding a mathematical model that best fits a given set of degradation data. In this context, commonly employed DMs to fit the data are of linear, logarithmic, power or exponential form [START_REF] Guo | Practical Approaches for Reliability Evaluation Using Degradation Data[END_REF]. For example, approximation of degradation model by a linear part and logarithmic/exponential part [START_REF] Jouin | Prognostics of PEM fuel cell in a particle filtering framework[END_REF], employment of exponential fit growth models [START_REF] Saha | Comparison of prognostic algorithms for estimating remaining useful life of batteries[END_REF], log-linear model for current drain degradation process [START_REF] Lu | Statistical inference of a time-to-failure distribution derived from linear degradation data[END_REF] and stochastic degradation model [START_REF] Si | A Wiener-process-based degradation model with a recursive filter algorithm for remaining useful life estimation[END_REF].

Prognostic approaches are broadly divided into three categories [START_REF] Sikorska | Prognostic modelling options for remaining useful life estimation by industry[END_REF][START_REF] Vachtsevanos | Intelligent Fault Diagnosis and Prognosis for Engineering Systems[END_REF]: model-based prognostics [START_REF] Chelidze | A dynamical systems approach to failure prognosis[END_REF], data-driven prognostics [START_REF] Kan | A review on prognostic techniques for non-stationary and non-linear rotating systems[END_REF][START_REF] Schwabacher | A survey of data-driven prognostics[END_REF] and hybrid prognostics [START_REF] Jouin | Prognostics of PEM fuel cell in a particle filtering framework[END_REF][START_REF] Neerukatti | Fatigue life prediction using hybrid prognosis for structural health monitoring[END_REF]. In model based approach, the degradation model is physics based and requires a detailed understanding of the underlying phenomenon [START_REF] Jardine | A review on machinery diagnostics and prognostics implementing condition-based maintenance[END_REF]. Inadequate modeling information, variation in behavioral physics or environmental conditions, un-modeled/unclassifiable sources of noise etc., result in limiting its adequacy.

Data-driven methods tend to learn the damage progression. However, they generalize damage progression over large sets of component population and remain unreliable in assessing the variability of damage progression trend from component to component in a population [START_REF] Jardine | A review on machinery diagnostics and prognostics implementing condition-based maintenance[END_REF]. As such, they provide inferior results especially in absence of complete data and large unit to unit variations. Hybrid approaches on the other hand, benefit from the fusion of the advantages of the former two [START_REF] Jouin | Prognostics of PEM fuel cell in a particle filtering framework[END_REF]. They employ physics or statistical based degradation models and use measured information to adapt the damage progression, accounting for un-modeled variations, environmental changes, external noise etc.

Prognostic approaches set as a joint state-parameter estimation problem [START_REF] Daigle | Model-based prognostics under limited sensing[END_REF], have been widely useful and may fall under hybrid approach wherein, the prediction of RUL is based on current estimate of damage state and state of damage propelling hidden parameters. Prediction of the RUL is obtained as probability distribution and accounts for the various uncertainties involved [START_REF] Daigle | Model-based prognostics under limited sensing[END_REF][START_REF] Daigle | A Model-Based Prognostics Approach Applied to Pneumatic Valves[END_REF][START_REF] Daigle | Model-based prognostics with concurrent damage progression processes[END_REF][START_REF] Roychoudhury | An integrated model-based diagnostic and prognostic framework[END_REF].

Choice of the filter for estimation and prediction process depends on the assumptions that can be made about the system, and desired performance [START_REF] Daigle | A comparison of filter-based approaches for model-based prognostics[END_REF]. Well-known Kalman filter, an optimal estimator for linear systems, has been used for prognostics in [START_REF] Celaya | A model-based prognostics methodology for electrolytic capacitors based on electrical overstress accelerated aging[END_REF][START_REF] Lim | Switching Kalman filter for failure prognostic[END_REF]. Extended Kalman Filter (EKF) [START_REF] Plett | Extended Kalman filtering for battery management systems of LiPB-based HEV battery packs: Part 3. State and parameter estimation[END_REF]or Unscented Kalman filter [START_REF] Daigle | Model-based prognostics of gear health using stochastic dynamical models[END_REF], may also be used for parameter estimation posing the problem as joint state-parameter estimation or as Expectation-Maximization problem [27] etc. However, they remain restricted to additive Gaussian noise. Also, EKF being sub-optimal diverges quickly if the initial estimate of state is significantly far from true value, or the model considered for estimation is not correct [START_REF] Saha | Comparison of prognostic algorithms for estimating remaining useful life of batteries[END_REF]. Comprehensive comparative studies of filters for prognostic purposes are found in [START_REF] Saha | Comparison of prognostic algorithms for estimating remaining useful life of batteries[END_REF][START_REF] Daigle | A comparison of filter-based approaches for model-based prognostics[END_REF][START_REF] An | Practical options for selecting data-driven or physics-based prognostics algorithms with reviews[END_REF].

Set in Monte-Carlo framework, Particle Filters (PF) or Sequential Monte Carlo methods [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF] form a suitable filter choice in this context, as it can be applied to non-linear systems corrupted with non-Gaussian noises for which, optimal solutions may be unavailable or intractable. Recently, particle filters have been exploited voraciously for prognostic methods [START_REF] Orchard | A particle filtering-based framework for on-line fault diagnosis and failure prognosis[END_REF][START_REF] Baraldi | Ensemble neural network-based particle filtering for prognostics[END_REF]. Significant works include prediction of end of discharge and EOL in lithium-ion batteries [START_REF] Saha | Modeling Li-ion battery capacity depletion in a particle filtering framework[END_REF], battery health monitoring [START_REF] Saha | Prognostics methods for battery health monitoring using a Bayesian framework, Instrumentation and Measurement[END_REF], prediction of battery grid corrosion [START_REF] Abbas | An intelligent diagnostic/prognostic framework for automotive electrical systems[END_REF], estimation and prediction of crack growth [START_REF] Bechhoefer | A method for generalized prognostics of a component using Paris law[END_REF][START_REF] Cadini | Monte Carlo-based filtering for fatigue crack growth estimation[END_REF][START_REF] Zio | Particle filtering prognostic estimation of the remaining useful life of nonlinear components[END_REF][START_REF] Zio | Prognostics and health management of industrial equipment[END_REF], fuel cell prognostics [START_REF] Jouin | Prognostics of PEM fuel cell in a particle filtering framework[END_REF], application to damage prognostics in pneumatic valve from the Space Shuttle cryogenic refueling system [START_REF] Daigle | Model-based prognostics under limited sensing[END_REF][START_REF] Daigle | A model-based prognostics approach applied to pneumatic valves[END_REF], estimation-prediction of wear as concurrent damage problem in centrifugal pumps with a variance control algorithm [START_REF] Daigle | Model-based prognostics with concurrent damage progression processes[END_REF], employment in distributed prognosis [START_REF] Roychoudhury | An integrated model-based diagnostic and prognostic framework[END_REF], exploring uncertainty management options for prognostics [START_REF] Baraldi | Investigation of uncertainty treatment capability of model-based and datadriven prognostic methods using simulated data[END_REF] etc.

Particle filters attract considerable attention [START_REF] An | Prognostics 101: A tutorial for particle filter-based prognostics algorithm using Matlab[END_REF], owing to the ever growing efforts being made for enhancement of performances and computational efficiency, such as the use of correction loops [START_REF] Orchard | Advances in uncertainty representation and management for particle filtering applied to prognostics, Prognostics and health management[END_REF], fixed-lag filters [START_REF] Daigle | Model-based prognostics with fixed-lag particle filters[END_REF] and the recently proposed adaption of the degradation model with a kernel smoothing method [START_REF] Hu | A particle filtering and kernel smoothing-based approach for new design component prognostics[END_REF]. Although a large amount of research exists in model based diagnostics and prognostics, very few promise the achievement/implementation of all key requirements in a common framework and the realization through a unified modeling paradigm [START_REF] Kunche | A review of PHM system's architectural frameworks[END_REF][START_REF] Moghaddass | An integrated framework for online diagnostic and prognostic health monitoring using a multistate deterioration process[END_REF]. In this context, Bond Graph (BG) modeling technique becomes a very suitable tool to deal with dynamic systems, especially those that belong to multi-energetic domains. A very brief and non-exhaustive introduction is given here. BG is a topological modeling language, where the exchange of energy between the different components of a dynamic system is captured in a graphical form. The energy exchange link is called a bond and there are two generic power variables named effort e and flow f, associated with every bond, such that e f Power  . The set of elements { , , } I C R , model the system parameters/component where I, C, and R are the inertial element, capacitance element and dissipation element respectively. The latter along with the elements {0,1, , } TF GY define the junction structure (global structure of the system) where TF and GY are the transformer element and gyrator element respectively. Junction 0 (or 1) implies that all the connected bonds have same effort (or flow) and the sum of flows (or efforts) equals zero. For efficient simulation of the physical behavior of the system, the computational order of the variables (e and f) must be decided systematically. For that purpose, the cause and effect decisions are described by the notion of causality. Causality in BG models is depicted by a perpendicular stroke on a bond. It determines whether the flow for a bond is computed from the effort or vice versa. The end of the bond which receives the effort is represented by the perpendicular stroke at that end of the bond. If all of the energy storage elements in a model are in integral form, the system is in integral causality. The constitutive equations of I and C respectively, in integral causality are 2 (linear case): ( ) (1/ ) ( )

f t I e t dt   and ( ) (1/ ) ( ) e t C f t dt   . For diagnosis task ,
where the initial conditions are unknown in real processes, the BG model is constructed in preferred derivative causality where the term "preferred" implies "wherever it is possible" . Derivative causality dictates the constitutive equation of I and C respectively, to be as (linear case): ( )

( ( )) / e t I d f t dt  and ( ) ( ( ))/ f t C d e t dt 
. For a detailed introduction from the ab initio and various related behavioral, structural and causal properties, the readers are referred to following works [START_REF] Karnopp | System Dynamics: Modeling, Simulation, and Control of Mechatronic Systems[END_REF][START_REF] Borutzky | Bond Graph Modelling of Engineering Systems[END_REF][START_REF] Mukherjee | Bond graph in modeling, simulation and fault identification[END_REF] and [START_REF] Thoma | Modelling and simulation in thermal and chemical engineering: A bond graph approach[END_REF].

For uncertain dynamic systems too, BG has been used extensively for modeling purposes and development of supervision techniques. This includes Fault Detection and Isolation (FDI) of complex systems [START_REF] Samantaray | Model-based process supervision[END_REF], highly non-linear and complex thermo-chemical systems [START_REF] Medjaher | Supervision of an industrial steam generator. Part II: Online implementation[END_REF], non-linear mechatronic systems [START_REF] Merzouki | Backlash fault detection in mechatronic system[END_REF], intelligent and autonomous systems [START_REF] Loureiro | Integration of Fault Diagnosis and Fault-Tolerant Control for Health Monitoring of a Class of MIMO Intelligent Autonomous Vehicles, Vehicular Technology[END_REF][START_REF] Benmoussa | Bond Graph Approach for Plant Fault Detection and Isolation: Application to Intelligent Autonomous Vehicle[END_REF][START_REF] Chatti | Functional and Behavior Models for the Supervision of an Intelligent and Autonomous System[END_REF],

industrial chemical reactors [START_REF] Harabi | Pseudo bond graph for fault detection and isolation of an industrial chemical reactor part I: bond graph modeling[END_REF], hybrid systems [START_REF] Ming | Prognosis of Hybrid Systems With Multiple Incipient Faults: Augmented Global Analytical Redundancy Relations Approach, Systems, Man and Cybernetics, Part A: Systems and Humans[END_REF] etc. In BG framework, the FDI is mainly based upon Analytical Redundancy Relations (ARRs) [START_REF] Bouamama | Derivation of constraint relations from bond graph models for fault detection and isolation[END_REF][START_REF] Samantaray | Model-based process supervision: a bond graph approach[END_REF] or by usage of the algebraic observers [START_REF] Benmoussa | Bond Graph Approach for Plant Fault Detection and Isolation: Application to Intelligent Autonomous Vehicle[END_REF][START_REF] Yang | New unknown input observer for control design: a bond graph approach[END_REF]. Specifically, for uncertain systems, BG in Linear Fractional Transformation (BG-LFT) [START_REF] Dauphin-Tanguy | How to model parameter uncertainties in a bond graph framework[END_REF][START_REF] Kam | Bond graph models of structured parameter uncertainties[END_REF] has been widely implemented for robust diagnosis by generating adaptive thresholds with respect to parametric uncertainties [START_REF] Djeziri | Robust fault diagnosis by using bond graph approach[END_REF][START_REF] Djeziri | Modelling and robust FDI of steam generator using uncertain bond graph model[END_REF][START_REF] Djeziri | LFT Bond Graph Model-Based Robust Fault Detection and Isolation[END_REF]. Very recently in Jha et al. [START_REF] Jha | Robust FDI based on LFT BG and relative activity at junction[END_REF], the authors have proposed modeling of parametric uncertainties in interval form. Unlike BG-LFT, where the threshold limits are simply the summation of the absolute values of each of the induced uncertain effort/flow at the junction [START_REF] Djeziri | Robust fault diagnosis by using bond graph approach[END_REF] leading to an overestimation of threshold bounds, the interval valued thresholds consider the sensitivity of each uncertain candidate to the respective residual. Even though there has been wide implementation of BG for robust diagnosis of complex processes [START_REF] Djeziri | Modelling and robust FDI of steam generator using uncertain bond graph model[END_REF], there have been very little efforts if none, towards the development/integration of prognostic techniques in BG framework.

The few motivations propelling the development of this work are:

 Initial steps towards system level prognostics in BG framework: There are many benefits of using BG, including but not limited to: systematic graphical representation of the governing differential equations, efficient decomposition of level etc. On the other hand, benefits of system level prognostics are many [START_REF] Sun | Benefits and challenges of system prognostics, Reliability[END_REF]. For almost all practical purposes, any plant (the dynamic system of interest) is a feedback closed loop system such that the system outputs follow a desired reference. As such, the system level prognostics present unique challenges in that incipient parametric degradation may progress unnoticed in presence of controller compensated system outputs, resulting in non-estimation of the same till the saturation limit of controller is reached. BG derived ARRs being sensitive to system parameters and control inputs can be exploited for the same at local component level while being in closed loop regime. Such a kind of BG enabled health monitoring, can be achieved in a unified framework at global system level.



Very few and inefficient existing residual based prognostic approaches: Most of the previous residual based attempts in BG framework consider damage progression deterministic in nature, incapable of adapting to the current damage progression and hence, reliability of prediction is minimal [START_REF] Medjaher | Hybrid prognostic method applied to mechatronic systems[END_REF][START_REF] Medjaher | Residual-based failure prognostic in dynamic systems, 7th[END_REF][START_REF] Djeziri | Data driven and model based fault prognosis applied to a mechatronic system[END_REF][START_REF] Djeziri | Fault prognosis based on fault reconstruction: Application to a mechatronic system[END_REF]. Moreover, uncertainties associated with measurements, operating conditions, process noise etc. have not been taken into account. This results in prediction of RUL without any associated confidence bounds, rendering it virtually useless for industrial certification and critical applications [START_REF] Uckun | Standardizing research methods for prognostics, Prognostics and Health Management[END_REF][START_REF] Saxena | Metrics for offline evaluation of prognostic performance[END_REF].

 Inclusion of the recently developed diagnostic methodology: Recently, [START_REF] Jha | Robust FDI based on LFT BG and relative activity at junction[END_REF][START_REF] Jha | Integrated Diagnosis and Prognosis of Uncertain Systems: A Bond Graph Approach[END_REF] proposed a methodology of modeling uncertainties in interval form and the subsequent diagnosis through interval valued thresholds. This generates a genuine interest in its inclusion in an appropriate health monitoring framework. This paper"s main objective is to address the problem of prognostics in BG modeling paradigm while the system is considered globally uncertain and parametric uncertainty is modeled in interval form. This is achieved by casting the problem as a joint state-parameter estimation problem, a hybrid prognostic approach, wherein the fault model is constructed by considering the statistical degradation model of the system parameter. The system parameter is known a priori to be undergoing degradation. Measurements are obtained from BG derived residuals (evaluation of ARRs). Using PF algorithms, estimation of state of the system parameter under degradation (prognostic candidate) along with the associated unknown hidden time varying Degradation Progression Parameters(s) (DPPs) is achieved and tracked to obtain the state of damage in probabilistic terms which is used for prediction of RUL of the system with respect to that parameter.

After this section, Section 2 details a fault detection algorithm where parametric uncertainty is modeled in interval form and Interval valued ARRs (I-ARRs) are derived systematically. The latter lead to the development of Interval valued robust thresholds over the nominal point-valued part of the I-ARRs. In Section 3, a novel methodology is proposed which includes construction of fault model and a novel way of obtaining the observation equation from the concerned nominal residual. In Section 4, the state of prognostic candidate and associated DPPs are estimated from the nominal part of the I-ARRs, sensitive to the latter and control inputs, in Monte-Carlo framework using PF algorithms. A novel variance control algorithm is proposed which ensures a suitable adaptation of random walk noise variance, once convergence is achieved in the estimation process. Prediction of the RUL is achieved in PF framework. Section 5 provides various evaluation metrics employed. In Section 6, the methodology is demonstrated through simulation and various issues are highlighted. In Section 7, the methodology is tested experimentally by variation of frictional torque on a mechatronic system and Section 8 draws the conclusions. The various novel contributions of the paper are listed as follows:

 Integration of BG modeling framework and Monte Carlo framework for estimation of state of health and prediction of RUL.

 Exploitation of nominal part of I-ARRs derived in (BG framework) for detection of degradation beginning and prognosis of incipient parametric degradation in Monte Carlo framework using PF.

 Obtaining the observation equation from the nominal part of I-ARRs and construction of local fault model such that state of the prognostic candidate and RUL prediction is obtained while system outputs are feedback controlled or otherwise.

 Accounting the various noises associated with degradation process and nominal residual output, for estimation and RUL prediction.

 Proportional control type variance adaptation algorithm with novel feedback condition that ensures a sustained convergence with low estimation variance (spread).

Major assumptions are:  Only the system parameters are considered uncertain. Sensors are considered non-faulty.

 The system parameter (prognostic candidate) that undergoes degradation is assumed to be known a priori. The issue of isolation or isolability of the faulty candidate is assumed resolved.

 DM of the prognostic candidate is assumed to be known a priori.

 Single fault (degradation) hypothesis is followed for an I-ARR considered.

 Noise associated with measurements (residuals) is assumed additive and normally distributed Gaussian in nature.

Bond Graph Based Uncertain System Modeling

Nominal model of any deterministic physical system may be modeled in BG form, in preferred integral causality, with nominal system parameters composed of basic elements

  n l n u ,,            (1)
    n n l u , , ,           (2)
Here, θ0 l  and Δθ u 0  are the additive uncertainty/deviation on the left and right sides, over the nominal value θ n .

[θ ,θ ]

nn is a degenerate interval with equal upper and lower bounds. For any additive uncertainty θ  over θ , the multipli- cative uncertainty is defined as:

n /      . Multiplicative interval uncertainty ,     is expressed as in (3) such that ,          .
Then, the uncertain θ may be expressed as shown in (4).

 

l n u n , / , /           + nn ,,                 (3) (4)
Interval valued parameters are represented on an uncertain BG closely following the BG-LFT representation details of which can be referred in [START_REF] Kam | Bond graph models of structured parameter uncertainties[END_REF]. For pedagogical illustration, a resistor element R in resistance (imposed flow) causality is considered.

 Nominal case (see Fig. 1a): The characteristic equation with parameter in nominal state (without any uncertainty) is expressed as:

RR e R. f  (5)
 Uncertain case (see Fig. 1b): With multiplicative interval uncertainty RR ,    , the characteristic law is expressed as:

  1+ R R R n R R R e ,e R,R . f R , . f                 (6)
    R ,unc R ,n R R n n R R e e
e ,e R ,R . f w     [START_REF] Gebraeel | Prognostic degradation models for computing and updating residual life distributions in a time-varying environment, Reliability[END_REF] where Similarly, interval uncertainty can be modeled and represented for the other BG elements I, C, GY, TF, RS etc. For diagnosis based on ARR generation, the detectors are dualized such that effort detector De becomes a source of effort signal SSe and imposes the effort signal at the 0-junction connected to the detector. Flow detector Df becomes a source of flow signal SSf and imposes flow at the 1-junction connected to the detector [START_REF] Samantaray | Diagnostic bond graphs for online fault detection and isolation[END_REF].

  R R R R R R n R w , .z , .R . f                 . Interval valued uncertain effort ,

Interval Valued ARR Generation

Classically, an ARR is a constraint relation derived from an over-constrained system/subsystem. It is expressed in terms of only known variables of the process [START_REF] Staroswiecki | Analytical redundancy relations for fault detection and isolation in algebraic dynamic systems[END_REF]. For any function f and set of known variables k, it has the form: ( ) 0  fk . In the context of BG modeling, an : ( ( ), ( ), ( ), ( ) ) 0

ARR f t t t t  θ SSe SSf
Se Sf , , where θ is vector of system parameters. For deterministic systems, the properties and ARR generation algorithm are detailed in [START_REF] Bouamama | Derivation of constraint relations from bond graph models for fault detection and isolation[END_REF]. BG model in preferred derivative causality with dualized sensors are utilized to avoid unknown initial condition problem. For the uncertain systems, robust FDI is achieved by generation of uncertain ARRs with perfectly separable nominal part and uncertain part [START_REF] Djeziri | Robust fault diagnosis by using bond graph approach[END_REF]. The approach of latter is described in [START_REF] Djeziri | Robust fault diagnosis by using bond graph approach[END_REF] and here, it is adapted to obtain interval valued ARRs (I-ARRs) in presence of interval valued uncertainties. Consider the uncertain parameter vector , m N

   θθ

and m NN   , following steps are taken to gen- erate I-ARRs:

Step 1: Preferred derivative causality is assigned to the nominal model.

Step 2: Parametric uncertainties are modeled in interval form and represented on the nominal BG, as explained in section 2.1.1 to obtain uncertain BG.

Step 3: The candidate ARRs are generated from "1" or "0" junction, where power conservation equation dictates that sum of efforts or flows, respectively, is equal to zero, as shown in ( 8) and ( 9) with s being the sign rendered to the bond due to energy convention.

 for 0-junction:s

  0 . , . : 0 m iN ii i s f f Sf s MSf w           (8)
 for 1-junction:

  0 . , . : 0 m iN ii i s e e Se s MSe w           (9)
Step 4: The unknown effort or flow variables are eliminated using covering causal paths from unknown variables to known (measured) variables (dualized detectors), to obtain the I-ARRs, , RR   consisting of known variables only as shown in [START_REF] Saha | Comparison of prognostic algorithms for estimating remaining useful life of batteries[END_REF]. The nominal part is characterized by point valued function 1  , with point valued nominal parameters as coefficients of point valued measured variables (cf.( 12)). It is separated from the interval valued part which is identified as an interval

function 2 Ψ (see Definition A.1 in Appendix A), sensitive to interval valued uncertainties (cf.( 13
)).     , : , , , , , ( ), (
)

n i R R w Se Sf t t          Ψ θ θ,θ SSe SSf (10)     1 2 ( ), ( ) : ( ) ( ), ( ) ( ) , ( ), ( ), , ( ), ( ) , , ( ), (
)

n n R t R t r t B t B t r t t t Se Sf B t B t t t                         n θθ θ Ψ θ,θ δ ,δ SSe SSf SSe SSf (11) (12) (13)
Hereafter, () n rtbeing the numerical evaluation of the point-valued nominal part 1  , will be referred to as nominal resid- ual and 2 Ψ will be termed as Uncertain Residual Interval Function (URIF).

Residual Based Determination of Degradation Model

In BG framework, the DM of a system parameter under degradation θ d  θ , θ N  θ can be obtained from the time evolution profile of the respective ARR to which it is sensitive, assuming that the rest of the system parameters which are sensitive to the same, do not undergo any kind of progressive fault or degradation [START_REF] Medjaher | Hybrid prognostic method applied to mechatronic systems[END_REF], [START_REF] Borutzky | Failure Prognosis for Hybrid Systems Based on ARR Residuals, Bond Graph Model-based Fault Diagnosis of Hybrid Systems[END_REF]. Here, consider the point valued part of the d th I-ARR ()

d rt , such that with \ θ ( ) d t   θθ , 0, ( ) 0 d n t r t  :   1 () θ ( ), , ( ), ( ), ( ), ( )
d d d n r t t t t t t  
θ SSe SSf Se Sf [START_REF] Chelidze | A dynamical systems approach to failure prognosis[END_REF] Here, the sub-script n denotes nominal value. The computed values of ()

  θ ( ) ( ), , ( ), ( ), ( ), (
)

dd dn t r t t t t t    θ SSe SSf Se Sf ( 15
)
Residual based DM should be obtained prior to prognostics i.e. prior to the phase when system"s health monitoring is of interest.

Fault Detection Using Interval Valued Thresholds

In this work, the fault detection module is constructed for the robust detection of degradation commencement.

Interval Valued Thresholds

Consider point valued parametric deviation(s

)   such that , , θ( ) , , t 
                              θθ , ()
bt is the numerical evaluation of function 2  with point valued arguments (cf.( 16)). Ψ can be expressed as finite sequence of interval arithmetic operations (evaluated as class code during implementation [START_REF] Moore | Introduction to interval analysis[END_REF]) so that it is considered as a Rational Interval Function of 2  and hence, is inclusion isotonic (see Definition A.3 and Definition A.4 in Appendix A). Then, through Fundamental Theorem of Interval Analysis (see Theorem A.1 in Appendix A) the inclusion of ( 17) can be verified.

  2 ( ) ,θ( ), ( ), ( ) b t t SSe t SSf t     (16)
( ) ( ), ( ) b t B t B t    (17)
Now, at all times, due to power conservation at the junction,

( ) ( ) 0 ( ) ( ) n n r t b t b t r t     (18)
Thus, the change in effort/flow brought by deviation at any instant in system is given by negative value of the nominal residual at that time. From ( 16),( 17) and( 18), following is used for fault detection:

Under nominal conditions:

( ) , ( ) n r t B B t    (19) Fault is detected if: ( ) , ( ) n r t B B t    (20)

Robust Fault Detection

The thresholds are generated by the range evaluation of URIF as shown in [START_REF] Daigle | A Model-Based Prognostics Approach Applied to Pneumatic Valves[END_REF] and [START_REF] Daigle | Model-based prognostics with concurrent damage progression processes[END_REF]. In discrete time step k, the algorithm for fault detection is given in Table I. It should be noted that therein, the bounds of URIF 2 Ψ is computed by expressing them as a sequence (computational graph or code list) of real valued functions [START_REF] Moore | Introduction to interval analysis[END_REF].

Fault model

In this work, the system parameter that undergoes degradation is assumed to be known a priori. Let θ ( )

d t
 θ be such a prognostic candidate. The objective in this paper is to estimate the state of θ ( )

d t based upon information (measurement)
provided by the values of nominal residual sensitive to θ ( )

d t , () d n
rt.

State Equation

The parameter under degradation θ ( )

d t is included as a tuple   θ , ,
d t : θ θ ( ) ( ( ), ( )) ;θ ( 0) θ d d d d d d n t g t t t    γ v ( 21
)
where, ()

d N d t  γ γ is DPP vector and θ θ () d d N t  v v
is the respective associated process noise vector. The fault model for

  θ , , d d d g γ
is constructed in state-space form by considering the parameter θ ( ) d t as the state variable augmented with the DPP vector as, ( ) ( ( ), ( ))

d d d d t t t  x x f x v (22)
where, () θ ( ), ( )

T d d d t t t    γ x
is the augmented state vector and d f is state transition function following the Markovian assumption.

Residual Based Observation Equation

Here, the objective is to exploit the nominal residual for the estimation of state variables. This way, the nominal residual used for detection of degradation beginning can be further used furthur for estimation of state of health of the prognostic candidate and associated DPPs. This is possible if the ARR expression is altered to obtain the observation equation in an appropriate way, such that the nominal residual provides the measurements of state variables. For this purpose, a simple algebraic approach is proposed.

Theorem: Under the single degradation hypothesis, assuming the nominal part () [START_REF] Lim | Switching Kalman filter for failure prognostic[END_REF] or,

d
        ( ) , ( ), ( ), , , A (θ ) A (θ ( )) A (θ ) 0 ( ) ( ) A (θ ( )) (θ ) 0 A( θ ( )) (θ ) ( ) d T d T d T d n n n d d T d d nn T d d d nn r t t t Se Sf t r t r t t t r t                θ SSe SSf φ φ φ φφ φφ (25)
Thus, degradation state θ ( ) d t can be linked implicitly to the measurements of ()

d
    1 () θ ( ) d n d rt a t    (27)
The observation equation argument in [START_REF] Daigle | Model-based prognostics of gear health using stochastic dynamical models[END_REF] includes known variables (sensor measurements, system parameters, inputs etc.) and their derivatives.

Degradation Estimation And RUL Prediction

In discrete time step k  , the fault model   θ , ,

d d d d k k k k   x x f x v   d d d d k k k y h w  x (29) (30)
where θ, y . In the following section, the method employed for degradation estimation and consequent prognostics is explained assuming that degradation begins at the start. In reality, information about k d will be given by fault detection module as described in Section 3.

T d d d k k k    γ x , d k f is

Degradation Estimation

In this section, concise details about PF are provided. The related concepts mentioned here can be found detailed in [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF] and [START_REF] Doucet | An introduction to sequential Monte Carlo methods[END_REF]. The state distribution is approximated by set of discrete weighted samples or particles,  

,,

), w

N d i d i i k k k i  γ
, where N is the total number of particles and for i th particle at time k, , θ di k is the estimate of the state (system faulty parameter here) and , di k γ is the estimate of fault progression parameters. The weight associated with each particle is denoted by w i k .

The posterior density at any time step k is approximated as,

0: (θ , ) 1 (θ , | ) w . ( θ ) d d k k N d d d i d d k k k k k k i p y d d     γ γ γ (32)
where (θ ,| )

d d d k k k py γ , assuming that particles   ,, 1 1 1 1 (θ , ), w
N d i d i i k k k i     γ
are available as realizations of posterior probability

1 1 0 : 1 (θ , | ) d d d k k k py    γ at time 1 k  , with the follow-
ing main steps:

 Realizations of prediction 0 : 1 (θ , | ) d d d k k k py  γ , is obtained in form of new set of particles   ,, 1 (θ , ), w
N d i d i i k k k i  γ , with
weights being chosen using the principle of importance sampling. The proposal importance density is chosen as the transitional prior

,,

| θ , ) / ( | θ , ) N i d d i d i d d j d j k k k k k k k j p y p y    γγ , , 1 w (
Note that with the choice of importance density as the prior, the weights were obtained as, ,, 1 w ( | θ , ) . The objective behind resampling is the elimination of particles with small weights and focus on particles with large weights, for estimation. In this work, systematic resampling scheme is preferred as it is easy to implement and takes O(N) time and the algorithm can be referred in [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF]. , dj

i i d d i d i k k k k k w p y   γ (34

 is modeled as a random walk process , , , 11

d j d j d j k k k     where, , 1 dj k
  is sampled from an artificial random zero-mean Gaussian distribution i.e. , makes the task complicated. In this paper, random walk variance is controlled in similar fashion as in [START_REF] Daigle | Model-based prognostics with concurrent damage progression processes[END_REF], however, with the distinguishing feature that variance adaptation is triggered by

, ˆdj k  : , , 0 , 1 mean( ) if ˆ1 mean( ) if l L dj dj k l l k dj k k L L k L                 (36)
with θ d at time k , being the average of the estimation mean , ˆdj  in a running window of previous L estimates. Fig. 2 shows the schematic of the proposed algorithm. The fact that degradation model of θ ( ) The corresponding pseudo-algorithm of the variance adaptation scheme followed in this paper is given in Table III.

d
Fig. 2. Schematic of variance control scheme

Remaining Useful Life Prediction

The critical/failure value of θ ( )

d d d i d i d d d i d i k l k l N d d d i d d kk k l k l k l k l i p y w d d          γ γγ 0: (θ , ) 1 (θ , | ) . ( θ )
where for the i th particle, the corresponding weight during the , di l -step propagation is kept equal to weight w i k . For i th particle, , θ , ,

di d i d i k RUL k l k l     ; the corresponding θ d k RUL is obtained as: θ, θθ 0: () 1 ( | ) w ( ) dd d i k N di k k k k RUL i p RUL y dRUL     (38)
The prediction IV.

Health Monitoring Algorithm

The beginning of degradation is detected by the fault detection module described in Section 2. Subsequently, the joint   decide the bounds of the uniform distribution as shown in [START_REF] Lee | Prognostics and health management design for rotary machinery systems-Reviews, methodology and applications[END_REF]. Such an approximation guarantees to include the true initial state of θ ( ) d t . The complete algorithm is shown in Table V. Fig. 3 shows the schematic description of the proposed methodology.

Fig. 3. Schematic description of the Health Monitoring Methodology

Evaluation Metrics

In this section, various metrics employed to evaluate the performance of estimation, prediction etc. are briefly discussed. [START_REF] Saxena | Metrics for offline evaluation of prognostic performance[END_REF] can be referred for details and [START_REF] Daigle | Model-based prognostics under limited sensing[END_REF][START_REF] Daigle | Model-based prognostics with concurrent damage progression processes[END_REF] for implementation of the same.

Estimation performance

The estimation performance is evaluated using two metrics that quantify the accuracy and spread.

Root mean square error (RMSE): This metric expresses the relative estimation accuracy as:   signifies the desired (pre-fixed) fraction of the RUL prediction probability mass percentage that must fall between the cones of accuracy determined by  , for the respective RUL prediction to be acceptable. In this paper, for all λ (all k), β =0.5 which translates to the requirement of 50% of probability mass distribution of

θ, d p k RUL fall- ing within ** θ , θ , [ (1 ) ,(1 ) ]
dd pp kk RUL RUL   for the prediction at k p , to be acceptable.

Case Study on Mechatronic System

The method presented in this paper is applied on a mechatronic Torsion Bar 1.0 system shown in Fig. 4. [START_REF] Tarasov | Flat control of a torsion bar with unknown input estimation[END_REF][START_REF] Kleijn | Torsion Bar 2.0 Reference Manual[END_REF]; it is integrated with 20SIM, a BG dedicated software [START_REF] Broenink | Modelling, simulation and analysis with 20-sim[END_REF]. Real time implementation is achieved through 20 SIM 4C 2.1, a prototyping environment that enables C-code implementation in real time on ARM-9 processor based torsion bar system [START_REF] Kleijn | SIM 4C 2.1 Reference Manual[END_REF]. The interval computations, estimation, variance control and prediction algorithms are written in Matlab Function Block in Simulink. The embedded code is generated through Simulink Coder in Matlab2013a ® . INTLAB [START_REF] Moore | Introduction to interval analysis[END_REF][START_REF] Rump | INTLAB-interval laboratory[END_REF] is used to implement interval calculations during simulation. For real time C-code generation, relevant/required functionalities are borrowed from INTLAB.

Nominal System

The schematic model of the mechatronic system (detailed in [START_REF] Kleijn | Torsion Bar 2.0 Reference Manual[END_REF]) is shown in The BG model of the nominal system in integral causality [START_REF] Tarasov | Flat control of a torsion bar with unknown input estimation[END_REF] is given in Fig. 6. The control input from PI controller (controlled variable: motor speed m  ) modulates the input voltage MSe:U PI . The measured angular velocities (obtained from angular position measurements) of motor shaft and load disk are represented as :

M Df  and :

Ld Df  respectively.

Belt is considered of high stiffness and the rigidity is not considered in the model. Also, the frictional loss due to the action of belt is lumped with frictional loss at motor bearing; it is modeled as a resistor element : is the ratio between number of teeth on motor disk to motor shaft [START_REF] Kleijn | Torsion Bar 2.0 Reference Manual[END_REF]. The electrical part of the DC motor is not monitorable as there is no sensor installed in it. Only the monitorable part (marked in Fig. 6) is used for analysis. It must be noted that the system is operating in feedback closed loop (Proportional-Integral (PI) control) regime. Analysis or development of the control strategy is not described, as the main interest of the paper does not lie in the same.

Uncertain BG and System Validation

The uncertain BG of the monitorable part in preferred derivative causality is shown in Fig. 7. The parametric uncertainties are modeled and represented in interval form. The global system is considered uncertain with uncertain parameter vector ,

  θθ : , , , , , , , , , , , , , , , , T m
m Md Md Md Md Ld Ld Ld Ld s s s s Jm Jm f f J J b b J J b b k k b b                                        θθ (44)
The monitorable part has input in form of the controlled electrical torque input generated by the DC motor. Both the sensors (Df) are dualized to corresponding source of flows as ( / )

m belt Ld f f f k      .

Interval valued ARRs and Robust Thresholds

From the steps described in Section 2, I-ARR can be generated from the detectable junction 1 1 of Fig. 7 as: ,( )

        ,, , , , , 1 , , , , , , , , . ( /
) , ( / ) ( / ) , .( / ) (1/ ) ( ) , m m m m Md Md Md Md ss PI m n m J J m n m m n m f f m n m Md n m belt J J Md n m belt Md n m belt b b Md n m belt belt m s n Ld k k belt R R J J f f J k J k b k b k k k dt k                                                  , , , ( ) ( )
ss m m m s n Ld s n Ld b b s n Ld belt belt belt k dt b b k k k                                          (
( ), ( ) n R R r t B t B t          )
where,

m m m m n in m n m m n m Md n Md n s n Ld s n Ld belt belt belt belt belt r t J f J b k dt b k k k k k                        (48)     ,, 1, , , , , , , 1 ( ) ( ) ( )
)

Md Md Md Md m m m m s s s s m m J J Md n b b Md n belt belt J J m n m f f m n m m m belt k k s n Ld b b s n Ld belt belt J b k k B t B t J f k k dt b k k                                                                  (49)
Only one I-ARR has been derived here at 1 1 ; it serves the purpose of approach-demonstration. Similarly, another independent I-ARR may be derived from 1 2 junction.

Study by Simulation

Nominal conditions

The nominal parameter values and respective multiplicative interval uncertainty is tabulated in Table VI. Fig. 8 shows

the
Md b t g  , is formulated as, 1, 1 , , 1 , 1 1, 1, 1 1, 1 , , , 1, 2 .
()

k Md t Md k Md k b k k k k mk d k Md k Md n k belt b b e v y b b w k                (53)
where, . The convergence is achieved very quickly but with large initial estimation spread. This is due to the high artificial noise variance set for the desirable quick convergence. As shown in Fig. 11b, the estimation spread is reduced (effective from t=20s) until 1 * v 10%   is achieved at around t=50s and thereafter, * 1  is tracked smoothly with controlled spread and

3.02% RMSE 

. For comparison, Fig. 11c shows performance with no variance adaptation. Therein, the estimation continues with large spread even after the convergence is achieved.

RUL prediction

Using 0.1, 0.5

 

and for all k  , RUL prediction is shown in Fig. 12. The RUL predictions obtained until t=52s

are not good and suffer with large variance spread due to the large corresponding spread in 1  (see Fig. 11a), making them virtually useless. However, after t=52s, the RUL distributions are well within accuracy cone (more than 50% of RUL probability mass lies within accuracy cone). Ignoring the initial period of convergence, the overall prediction performance is very good with 9.8% RUL RMAD  and 97.15% RA  .

A Qualitative Analysis

As seen in the previous section, accuracy and spread of RUL predictions are directly influenced by the estimation quality of DPP. This in turn, depends on choice/initial setting and tuning of the several parameters involved. They are discussed here qualitatively. Note that estimation obtained in Fig. 11a forms the most desirable performance. In this section, only the specified PF parameters are played with; rests are kept same as for Fig. 11a. Although some of the things discussed are intuitive from the perspective of estimation in a state space model, authors have felt the necessity to highlight their concurrence when the residual is used as measurement.

    . It is the most important factor that determines good tuning of parameters in succession. A very high value of the latter, results in bad estimation performance. Fig. 13a shows estimation with  is seen, the estimation continues further with a very wide spread for a long time before it is gradually reduced, owing to variance adaptation scheme. On the contrary, a very low variance will result in very late convergence, if at all.  Proportional Gain P: Proportional gain determines how fast the estimation spread is reduced to the reference * v 

. As observed in Fig. 11a, an appropriate choice of latter was found as P 1 = 0.001. It resulted in smooth tracking after convergence was achieved. A high gain value results in quick reduction of estimation spread; however, it is accompanied with continuous shrink and expansion as shown in Fig. 13c with P 1 =0.005. The latter has also been demonstrated in [START_REF] Daigle | Model-based prognostics with concurrent damage progression processes[END_REF]. Although, a very high gain value may bring down variance spread quickly; however, it may be followed by poor convergence results as shown in Figs. 14c andd, with P 1 =0.01. On the contrary, a very low P 1 renders a non-effective variance adaptation as shown in Fig. 13d with P 1 = 0.0001, adding no significant benefits in RUL prediction.

 Desired RMAD ( * v  ): The pre-fixed 1 * v  for 1  , determines how much freedom is given to 1  after the estima- tion spread is brought under control. An appropriate choice of 1 * v  gives enough freedom for convergence even after actual variance is well under 1 * v  , as shown in Fig. 11a between t=50s and t=80s with 1 * v  =10%.

In extreme cases, where P 1 is chosen of high value (rate of RMAD reduction depends on P 1 ) and 1 * v  is set very low, the estimation may remain stagnant near, but not equal to * 1  . This is shown in Fig. 14c with P 1 =0.01 (read high) and desired RMAD 1 * v  =6% (read very low).

 True DPP interval

** Here, the variance adaption is not effective enough. On the contrary, if the interval width is appropriately set (assuming that initial estimate is outside of it), 1, ˆk  is captured quickly and variance control is triggered early, as shown in Fig. 14b with  Residual noise variance (measurement noise) assumed by PF: Noise that corrupts the residual measurements can be non-Gaussian due to presence of derivative terms. Such noises can be dealt by PF effectively without any restrictions. In this work, the explicit distribution of the residual noise is not found. Instead, it is approximated as normally distributed Gaussian in nature. The related standard deviation and variance is found out from residual measurements. Moreover, generally, the variance of measurement noise (residual noise here) assumed by PF, is greater than the approximated measurement noise. This is done to counter the sample impoverishment problem which happens while very few particles have significant weights and most of the other particles with nonsignificant weights are abandoned during the resampling process [START_REF] Li | Fight sample degeneracy and impoverishment in particle filters: A review of intelligent approaches[END_REF]. Higher residual noise variance assumed by PF allows higher particles being sampled for estimation, thus, reducing the problem of sample degeneracy and consequent impoverishment. As followed in other works [START_REF] Daigle | Model-based prognostics with concurrent damage progression processes[END_REF], in this work too, the residual noise assumed by PF is greater than actual residual noise.

  ** 1, 1, , 0.01

Computational complexity

The time taken per step for estimation and RUL prediction depends on the number of particles used. With N=500, on an average, 0.03s was consumed per step. Fig. 15 shows the RUL prediction computation time per step for the RUL prediction performance of Fig. 12. In addition to the number of particles N, computational time for RUL prediction varies:

 Inversely with the time at which prediction is made: The farther is the time from EOL at which RUL prediction is made, the longer it takes to simulate to EOL. This makes the computational time large.

 Inversely with estimated DPP  : At a certain time of prediction, higher is the rate of damage progression, smaller is time taken to simulate to EOL. As seen in Fig. 11a, before t=50s, the estimation value of 1  is lower than true value accompanied with large variance. Therefore, for a specific N, the computation time per step before t=50s is higher and with large variations. After t=50s, with a nearly uniform 1  estimation and lesser spread (see Fig. 11a), the computation time follows an almost uniform monotonic decreasing trend (see Fig. 15).

Simulations were run on a 2.49-GHz dual core processor with 8GB RAM. With N=500, and sample time of 0.1s (which translates to 10 computational steps per second); it took on an average 32 minutes to simulate system dynamics, estimation and RUL prediction till 100s. With N=50, the same took 110 seconds. This indicates that through employment of lesser number of particles, the RUL predictions could be achieved in real time, for experimental purposes. Moreover, for real experiments run on complied C, the run time reduces drastically by an order of magnitude.

with Md r as the radius of the motor disk. In the BG model, it is incorporated as non-linear resistance element at motor disk as shown in [START_REF] Benmoussa | Bond Graph Approach for Plant Fault Detection and Isolation: Application to Intelligent Autonomous Vehicle[END_REF]; the corresponding characteristic equation is obtained as shown in [START_REF] Chatti | Functional and Behavior Models for the Supervision of an Intelligent and Autonomous System[END_REF]

        (56)
The corresponding I-ARR

                       (57)     , , , , , , 2 , , , , , sgn( / 1
m m J J Md n b b Md n n n Md m belt belt belt J J m n m f f m n m m m belt k k s n Ld R R s n belt be J b M gr k k k B t B t J f k k dt R k k                                                                     ) Ld lt               (58)
The nominal value of  , n  is found out from 2, () n rt and ( 57) by suspending a known load mass. Fig. 16 shows the residual profile under nominal conditions. Fig. 17 shows the effect of adding load (or frictional toque) in a discrete way on the system.

Md

 is controlled at 30 rad/s. Each time load is added, there is PI controller enabled compensation due to which Md  settles to the reference velocity. However, The experiments involve only non-destructive procedures so that there is no degradation (wear) of the surfaces. In other words,  is assumed constant. Experiments involve variation of suspended load mass M in a uniform way till the limit fail M , is reached. () Mtis treated as system parameter under degradation. The experiments were conducted in two distinct phases:

 Offline: In this phase, multiple tests were done with the load being added uniformly. As explained in Section 2.2, variations of () Mt were obtained from the evolution of 2, () n rt found in [START_REF] Harabi | Pseudo bond graph for fault detection and isolation of an industrial chemical reactor part I: bond graph modeling[END_REF]. This provided the time dependent DM of the system parameter () Mt .  Online health monitoring: The maximum limit of additive load mass fail M was pre-decided keeping in mind the safety of the system. Load was varied until fail M ; this was performed in the similar environment as of the offline phase. In real time, estimation of () Mt and associated DPPs, and subsequent RUL predictions were obtained.

Case I :Linear variation of mass

Linear degradation models are frequently employed where incipient degradation does not accelerate subsequent degradation. Here, such a scenario was created through experiments and tested in real time.

Degradation Test and Degradation model

Load is varied linearly. Ten experiments are carried out wherein; sand is poured with same environmental conditions to maintain the uniformity. Fig. 18 shows the experimental data and the data mean found at each instant. A linear fit over data mean is obtained using linear regression. The DM can be expressed as,

1 2, 1 2 2, 2, 2, k k k k M k k k M M t v             (60)
Observation equation is obtained from the nominal part of I-ARR

) ( ) ( ) n Md Md k k n k k k n k belt gr y r w t M M w k            (61) , 2, 2, , 2, 2, sgn( )
For the experiment, load is varied until ( )

1.5 fail M t M Kg 
. Fig. 19 shows the nominal residual profile under degradation.

Estimation

The prognostics module is triggered at t=22 s; estimation and predictions are performed with N=50 particles, t  =0.1s,  is not guaranteed to remain constant; the DM provides only an approximate idea of its magnitude order. Fig. 20b shows the estimation with large initial variance. The estimation spread is reduced effectively from t=40s.

Thereafter, the estimation mean remains around * 2  with RMAD of 6%.

RUL prediction

Prediction of RUL is shown in Fig. 21

Fault Model

The tuple ( () Mt, 3  , 3 g )is formulated in state space as,

) ( ) ( ) k t k k M k k k k n Md Md k k n k k k n k belt M M e v gr y r w t M M w k                      3, 1 1 3, 1 3, 3, 1 3, 1 , 3, 2, , 3, 3, . sgn( )

Estimation and RUL Prediction

The prognostic module is triggered at t=22 s. It is performed with N=50, t  =0.1s, . In fact, the EOL is achieved slightly before than that predicted by DM.

It should be noted that RMSE  in real time experiments is higher than that obtained in simulations as *  does not remaining perfectly constant in real cases. Also, usage of lesser number of particles leads to worse estimation performance.

However, overall prediction and estimation performances are very good and satisfactory. model of the uncertain system can be used for detection of system parameter"s degradation. Subsequent estimation of the state of health and associated degradation progression parameter(s), and prediction of the remaining useful life of the prognostic candidate can be obtained using particle filtering algorithms. This leads to an efficient integration of the benefits of BG modeling framework and Monte Carlo framework. The uncertain part of the I-ARRs is used for robust threshold generation over the nominal part. This enables efficient detection of the degradation commencement, robust to parametric uncertainty. Further, the same nominal residual can be used for obtaining the measurements of state variables in the fault model while the observation equation is developed from the nominal part of the I-ARR. For the latter, a novel algebraic approach is proposed so that the robust detection of degradation and further estimation of state variables of the fault model can be achieved using the same nominal residual in a unified framework. Moreover, this methodology can be extended in presence of multiple degradations which forms a potential future work. In future, the work will be effectively explored for large systems with multiple prognostic candidates. Being sensitive to the control inputs, nominal residual is able to capture the parametric degradation profile even while the system outputs remain in feedback closed loop regime.

This makes the approach effective for system level health management. Approximation of the distribution of noise present in residuals can be difficult or impossible, due to presence of derivative or integral terms in the arguments. As such, employed Particle filter algorithms form the best choice in this regard, supporting non-Gaussian noises. The novel variance adaptation scheme leads to very good estimation results and involves less complexity in terms of tuning of the involved factors. In future, the latter will be developed further and exploited for similar purposes. Through simulations, this approach has the capability of generating long term and very long term predictions.

Through experiments, capability of obtaining RUL predictions in real time has been shown, although, in very short time window. The associated computational complexity prevents the long and very long-term RUL predictions in real time.

In future, additional ways to obtain the same in sliding time windows will be explored. The method will be extended to achieve very long term predictions in multiple stages, comprising of small time windows, in real time. Although, robustness of the methodology has not been analyzed quantitatively, a qualitative analysis has been presented which helps in an efficient tuning of the PF parameters. As this work forms an effective initial step towards prognostics in BG framework, the same methodology will be applied over complex non-linear thermochemical-hydraulic systems such as fuel cells and vapor generator systems.

:

  MSf w  ) are added to represent the introduction of an additional uncertain effort (or uncertain flow) generated by the interval uncertainty on the system. The virtual detectors * De (or * Df ) are used to represent the information exchange/transfer. The star "*" is added as super-script for distinguishing the fictitious detectors (signals) from the real ones. In general, symbol ":" is used alongside a generic BG element to indicate the value in its respective characteristic equation. For instance, (see Fig. 1b) : n RR indicates that the system component modeled as resistor R has the resistance value of n R in its characteristic equation

.

  As it is clear, the associated effort (or flow) information z.

  the virtual detector De * .For better illustration of the effort/flow transfer, Fig.1cshows the equivalent block diagram representation of the uncertain BG in Fig.1b.

Fig. 1 .

 1 Fig.1. (a). Nominal R element (resistance causality), (b) Uncertain R element (resistance causality) in Interval form, (c) Equivalent Block Diagram Representation of Uncertain R element (Illustration of Signal Transfer)

d

  rtat time sample points give an implicit relation of the degradation profile of θ ( ) d t in time. Assuming that implicit function theorem is satisfied [79], (14) gives a real valued function d  such that:

2 Ψ 2 

 22 can be considered as the Natural Interval Extension Function of point valued function (see Definition A.2 in Ap- pendix A), with the point valued arguments and operators replaced by the corresponding interval arguments (timeinvariant here) and interval operators in the syntactic expression of the function 2  [80, 81]. Then, 2

  the linear/non-linear degradation progression function (DPF) obtained from the corresponding DM. The latter models the way degradation progresses in θ ( )

γ.

  denotes the Dirac delta function located at (θ , ) dd kk γ and sum of the weights In this paper, sampling importance resampling (SIR) PF is employed for estimation of 0:

  that particles are generated by sampling from probability distribution of system noise 1 updated. The weight w i k is associated to each of the particles based on the likelihood of observation d k y made at time k as, , ,

)

  To avoid the degeneracy problem, a new set of particles is resampled (with replacement) from the approximation of

X

  artificial random walk noise permits the estimation of θ ( ) d t to converge to its true value during the estimation process. Selection of the variance of the random walk noise is essentially a tradeoff between values that are big enough to allow the convergence in reasonable amount of time, yet small enough to let the parameter values be tracked smoothly once convergence is reached[START_REF] Daigle | Model-based prognostics under limited sensing[END_REF]. One of the efficient ways of ensuring good estimation of θ ( ) d t is to reduce the random walk noise variance , a suitable convergence is reached. In this regard, performance enhancement has been achieved by the usage of proportional control law type variance adaptation method; it is proposed, demonstrated and implemented in[START_REF] Daigle | Model-based prognostics with concurrent damage progression processes[END_REF]. Therein,Variance (spread) is quantified by the statistically robust metric Relative Median Absolute Deviation (RMAD) is an element for a data set X.The variance is adapted in a proportional control law way where the normalized error between the current RMAD , g. 10%) is multiplied by a proportional gain constant , is then increased or decreased by that amount. Thereafter, current random walk noise ,dj k is sampled from a zero mean Gaussian distribution with the modified variance , -in, the adaptation that progresses in arbitrarily decided multiple stages, requires a proper tuning of refereach stage. Such a procedure can be a tedious task especially in presence of multiple DPP. Although the objective that rests in achieving proper convergence and subsequent smooth tracking is clearly achievable, availability of no guidelines for a proper selection of number of stages and ,

t

  is known, leads to an approxi- mate knowledge of the true value of ,* dj can be decided based upon the approximate knowledge of ,* dj , obtained from the DM. The main objective rests in letting the variance be regulated in an automatic way.

dt

  is θ d fail ; it is specified beforehand. The corresponding RUL prediction at time k is θ is obtained with d l =1,… d Tk  , where d T is the time horizon of interest i.e. time until θθ ahead state distribution is computed by propagating each of the particles

y

   which are not available. Pseudo algorithm for RUL prediction is given in Table

Fig. 5 ..

 5 It consists of the Maxon® servo motor that provides the controlled actuation (rotation) to the disks; it is equipped with voltage amplifier A m , inductance La, resistance Ra, rotor inertia J m . The associated motor friction coefficient is m f and torque constant is k m . The high stiffness transmission belt provides the torque transmission with the transmission ratio of belt k , to the motor disk with rotational inertia Md J .The motor disk is connected to load disk with rotational inertia Ld J , through a flexible shaft that constitutes the drive train. The shaft is modeled as spring-damper element having damping coefficient s b and spring constant as s k The friction in the bearings of the motor disk and load disk is modeled as viscous friction with damping parameters as Md b and Ld b , respectively. Friction arising due to belt action is lumped with the viscous friction coefficient at motor disk in Md b . The setup is equipped with motor encoder and load encoder that measure, respectively, the angular position of motor shaft and load disk (2000 pulses per revolution). Angular position of the motor disk is obtained by dividing the motor encoder counts by belt ratio.

MdRb.

  GY element models the conversion of electrical current to electrical torque in the DC motor with corresponding coefficient of gyration being GY m mk  . TF element models the transmission of velocity through the belt from motor shaft to the motor disk. The corre-

Fig. 4 .Fig. 5 .Fig. 7 .

 457 Fig. 4. (a). Mechatronic Torsion Bar 1.0 system Fig. 4. (b). Fabricated Mechanical Lever type arrangement for Load (Mass) Suspension

.

  C element remains in integral causality with the initial condition given by the flow at respective 0-junction, provided by encoder readings as 10 9 13

2 (

 2 the interval threshold bounds determined in[START_REF] Mukherjee | Bond graph in modeling, simulation and fault identification[END_REF].

Fig. 8 Fig. 9 . 6 . 3 . 2 .

 89632 Fig. 8 Nominal Conditions (a) Motor disk speed (b) Input voltage (c) Nominal residual and Interval valued thresholds Fig. 9. Simulation of Degradation (a) Injected Degradation (b) Motor Disk Speed (c) Input Voltage to the System (d) Nominal Residual 1, () n rt

1 

 1 additive random walk noise. The estimation of state of parameter () Md bt is triggered at t d =10s. Initial estimate , 10 ~(0.045, 0.055) .s/rad. The true value of DPP is kept as * =0.05 Nm/rad so that , Md stop  is reached at 100s. Here, 0.1 ts  and N=500. Simulation is run until t f =100s when ˆMd b reaches the failure value , fact, estimation spread decreases as the estimation progresses, indicating the desirable performance. Estimation of Md b largely depends upon quality of estimation achieved with 1  . Fig. 11 shows the estimation of 1 ./rad, proportional gain P 1 = 0.001, assumes measurement noise variance equal to 4 times that of residual noise variance

Fig. 12 .

 12 Fig. 12. RUL Prediction performance with respect to estimation in Fig. 10 and Fig. 11. (a).

Fig. 10 .Fig. 11 .

 1011 Fig. 10. State estimation of the prognostic candidate system param-

(

  read high) wherein, although, quick convergence of mean 1

Fig

  late convergence. For tuning of other related parameters in this paper, an initial high value of variance

  : 1 

 1 The main objective of the latter remains in triggering the variance adaptation. As such, if width of band. This may lead to a very insignificant effect of variance adaptation on the estimation performance. Fig.14ashows the estimation with , which can considered "too tight" around * =0.05 Nm/rad.

  . This leads to early reduction in variance. However, a bad choice of P 1 (read high) and early variance adaptation, may lead to a rapid reduction in spread, followed by stagnation of estimation around ,* dj , before converging slowly to the same, as shown in Fig.14dwith P 1 = 0

Fig. 13 . 0001 Fig. 14 .Fig. 15 Fig 7 .

 13000114157 Fig. 13. Estimation of DPP 1  for qualitative analysis (a) Large initial random walk variance 1 ,0 0.04

2 Fig. 18 Degradation 2 2

 2182 Fig. 18 Degradation Test Data (linear variation) Fig.19 Nominal residual 2, () n rt  while system is under degradation. (linear variation of mass)

1

 1

 2 w

 2 =5x10 -3 V. For estimation, particle filter assumes measurement noise variance 9 times that of measurement variance 2 to counter sample impoverishment problem. Estimation of M is shown in Fig. 20a. The true M * is the residual based measurement of M(t) (as described in section 2.2, cf. (15)). State is estimated very accurately with M RMSE = 3.98%.

Fig. 20 .Fig 21 : 2 

 20212 Fig. 20. (a). Estimation of M (b) Estimation of 2 

7. 2 . 3 Fig. 22 Fig

 2322 Fig. 22a shows the experimental data. Fig. 22b shows the exponential fit over the experimental data mean from which the approximate value of DPP * 3  =0.05 Kg/s, is obtained. Regression residuals provide

 and the approximation of 3 w.

 3 is determined from 2, () n rt values during degradation tests. The structure of the observation equation remains same as in (62). For the experiment, load mass is varied until Fig 23 shows the profile of nominal residual under exponential degradation.

3 

 3   =4x10 -6 and 3 w  =5x10 -3 V.For estimation, particle filter assumes measurement noise variance 9 times that of measurement variance 2 3 w  to counter sample impoverishment problem during the experimentation. As shown in Fig.24a, state of parameter is estimated accu-=[1x10 -2 , 9x10 -2 ] Kg/s. Estimation is achieved with 3 RMSE  =7.6%. It must be noted that in reality, * cannot be claimed to be the accurate true value of 3  . Fig.24cshows the RUL initial predictions until t=32 (due to large spread),

Fig. 24 .

 24 Fig. 24. (a) State estimation for Trail 1 (b) Estimation of DPP 3  (c) RUL prediction for case II

Algorithm 5 :

 5 Health monitoring of 0 θ d with respect to d n r while system is running do Detect the beginning of degradation using Algorithm 1 if fault detection =true then /

, I, R TY,GY with

											 θ CN   ,   θ	. Sub-script n denotes the nom-
	inal value of the parameters. The sensor vector is formed by ( ) [ ( ), ( )] T t t t  Y De Df	with	De	() t 	N	De	being effort sensor
	vector and () t  Df	N	Df	being the flow sensor vector. The control/input vector is formed by ( ) [ ( ), ( )] T t t t  U Se Sf	with
	Se	() t 	N	Se	and	Sf	() t 	N	Sf	being respectively the source of effort and source of flow vectors. There is no uncertainty
	considered on system input (actuator/load/control input). The global system is considered uncertain with system parame-
	ters in interval form		,    θθ	N	m	where m NN   . The system parameter θ is modeled in interval form as θ θ,θ    ,
	where θ θ,θ    	θ θ θ   	and θ,θ     ,      θθ . Then,

2.1.1. Uncertainty Modeling and Representation on BG

	Uncertain system parameters can be represented in interval form on an uncertain BG by decoupling the nominal pa-
	rameter value	θ	n	  C , I , R TY ,GY , from the uncertain interval part    , n n n nn   where for the notational simplicity, n
	        ,   	. The additional uncertain effort (or flow) is brought-in at the junction 1 (or 0) by interval uncertainty	i     .
	It is represented on uncertain BG by a combination of: virtual effort (or flow) detectors	θ De z (or *:	θ Df z ) and fictitious *:
	source of effort input

   power conservation at the BG junction where the corresponding I-ARR is derived, gives,

	where | 1, 2... i i m 	,	A	1 mT 12 [ ... ] m a a a  		is a vector of known (measured system variables) with
	( ,   θ SSe 	( ),	SSf	( ), ( ), ( )) Se Sf	and	φ	1 	12 (θ ( )) [ (θ ( )), (θ ( )),.... (θ ( ))]    	is the vector of non-linear
	n rt , the measurement of the state θ ( ) d t \ θ ( ) d t   linear combination of non-linear functions of θ ( ) d t value of () d n rt. Proof: Let θ ( ) d t be the prognostic candidate and θθ . Assuming () d n rtcan be expressed as, of an I-ARR can be expressed as a can be obtained from the negative   ( ) ,SSe( ), ( ), ( ), ( ) A (θ) d T d n n n r t t t t t     θφ SSf Se Sf (23) i i n a t t t t m d d d d T m t t t t functions of θ ( ) 0 t d T d n r t t t t t t d t . Then,     ( ) ,SSe( ), ( ), ( ), ( ) A θ ( ) 0      θφ SSf Se Sf

  state transition function (possibly non-linear) and is described by first order Markov model.

	problem in particle filter (PF) framework, where the estimation at time k is obtained as probability density function (pdf)
	(θ , | d d k k py γ	: k k d d	)
	( )	( )  x	( )	(28)
	Measurements   | θ , γ d d d k k k p y d d d y t h	d k  1 2 y are assumed conditionally independent, given the state process       2 2 exp 2     d k d k d d d k k w w y h x d t w t	x .The likelihood function becomes as, (31) d k

With the beginning of degradation being detected by the FDI module as a fault at time step d k , the prediction of EOL/RUL at prediction time k , requires the estimate of θ d k , d k γ . This problem is cast as joint state-parameter estimation   , based upon history of measurements from the time of beginning of degradation k d up to k , : d d kk

Random Walk Noise Variance Adaptation

	Consider the DPP vector	γ	N d  γ	d	such that	{1,.. } d jN  γ	, , d   γ , , j d ˆdj  is the estimated value, and
		d	, * j	*  γ * , N d d γγ	d	being the respective true value. Also consider the interval vector	   d* d* lu γ , γ	N	γ	d	, consisting of
	intervals	, * j lu , d d   	, * j	    	d* lu d* γ , γ	 	, that contain the true value , * d j 	, * j lu , d d   , * j    . Moreover, for every	, d j   γ , con-d
	sider an associated constant (proportional gain)	P	, dj	such that	, P  P , d j d

The prediction, update and resample procedures form a single iteration step and are applied at each time step k. The algorithm for SIR filter is given in Table

II

. Details about other variants of sequential importance sampling PFs can be referred in

[START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF]

. 4.1.1. d N d  γ P .

  The latter is the time point at which system is considered to have obtained the failure state. The residual is sensitive to the input torque and hence, the input voltage captures the degradation evolution throughout the system"s lifetime.6.3.3. Fault modelState measurement is obtained from the observation equation which is developed using the Nominal Part of I-ARR

	t=100s. 1 [ , ] RR , 1, () n rt (cf. (25) and (27)) as,
	0		1, r t ( ) ( ( ) n Md b t b Md n ,  	).	    1, () n Md rt b  		(51)
	y	1, dd ( ) n r t w t ( ) ( ( ) Md b t b Md n ,     	)	2 m belt () t  k 		1 w t ( )	(52)
	where	1 2 1 ( ) ~(0, ) w wt  approximates the noise which corrupts 1, () n rt 	. For estimation, the fault model denoted as tuple
	11 ( ( ), , )
									, Md stop	=3rad/s at

The fault is detected at t d =10s when residual crosses the interval thresholds. Note that Md  is controlled at 30 rad/s until t=44.2s while the PI controller is effective. Thereafter, as the saturation value of actuator (motor) input voltage (12V) is reached, the speed Md  starts to decrease and reaches

  .

	R	b 	( ). t r  / | | M g						(55)
		Md			Md
	88 ( ) e R f	Md b	Md	M	Md r ( ). t	g	(	Md	/ |	Md	|)

  For safety reasons, the disk is kept at stop condition for few seconds after which the load is removed; this brings back the controller action into play. It is clearly visible that residual captures the variation of friction (variation of mass) while controller remains effective or otherwise.

		n 2, () rt	being sensitive to increase in current (and thus, voltage) decreases
	and settles to a different value. Saturation value for input voltage is reached around t=65s as the total load suspended is 1.6
	Kg. Thereafter (t>65s), controller is unable to compensate the change in

Md

 . Addition of more sand leads to reduction in motor disk speed; it stops at around t=70s.

  During the last 3 seconds of experimentation, the sand inflow is stopped gradually (and not abruptly) bringing in certain non-uniformity. As such, RUL predictions at t=58s, 59s and 60s, do not fall under the

							with	 	0.2	and	 	0.5	.The initial predictions have a very large spread due to
	the large corresponding spread in 2  . However, after t=35s, the RUL is within the	(1			* ) RUL	bounds with
	RA 	98.64%	,	9.4%
	(1			* ) RUL	bounds that are based upon the ideal linear degradation model.

RUL RMAD

 .

Table I Algorithm 1: Fault detection with d th I-ARR

 I

	Input:	    	 FC_ 1 i 	 ( ), ( ) kk ( ), , kk ( ) SSf   θθ Se Sf, SSe  2 , , i n θ       Ψ θ,θ , δ ,δ SSe SSf		    
	Output: fault detection
	n r	( ) k		1		( ) k Se, Sf, , SSe , SSf n θ	( ) k	
	( ), ( ) k B k  B 		Ψ	2	k     k       θθ θ,θ , δ ,δ ,SSe( ),   SSf( )
	if	( ) k  r	B	( ) k	and	( ) k  r	B	( ) k

fault detection  false else fault detection  true end if Table II Algorithm 2: Estimation using SIR filter

	Inputs:  (θ , ,, 1 1 d i d k k   γ	i	), w	i k	1 	 i N  1	k , d y
	Output:  (θ , ,, d i d k k γ	i	), w	i k	 i  1 N
	for i=1 to N do
	γ	d k	,	i	~( p	γ	d k	,	i	|	γ	d k	,  1 i	)
	, θ ~(θ | θ , , , 11 , d i d i d i d i k k k k p  γ	)
	,, w ~( | θ , ) i d d i d i py γ
			k					k					k	k
	end for
							N
	W		1    i	w	i k
	for i=1 to N do
	w	w / ii kk W 
	end for
	 (θ , , d i k	 ii    , , 11 RESAMPLE (θ , ), w NN d i d i i k k k  γγ , ), w d i i k k

Table III Algorithm 3 :

 III3 Adaptation

Acknowledgements

The first author is thankful to Prof. Genevieve Dauphin-Tanguy and Ecole Centrale de Lille, for funding the research work.

Definition A.4:

A rational interval function is an interval-valued function whose values are defined by a specific finite sequence of interval arithmetic operations.

Lemma A.3.1: All rational interval functions are inclusion isotonic. Theorem A.1 (Fundamental Theorem of Interval Analysis): If F is an inclusion isotonic interval extension of f , then

Corollary A.1.1: If F is a rational interval function and an interval extension of f, then ( , ... ) ( , ... ) nn f X X X F X X X  .