=Paper=
{{Paper
|id=Vol-1507/dx15paper11
|storemode=property
|title=Condition-based Monitoring and Prognosis in an Error-Bounded Framework
|pdfUrl=https://ceur-ws.org/Vol-1507/dx15paper11.pdf
|volume=Vol-1507
|dblpUrl=https://dblp.org/rec/conf/safeprocess/Trave-MassuyesP15
}}
==Condition-based Monitoring and Prognosis in an Error-Bounded Framework==
https://ceur-ws.org/Vol-1507/dx15paper11.pdf
Proceedings of the 26th International Workshop on Principles of Diagnosis
Condition-based Monitoring and Prognosis in an Error-Bounded Framework
L. Travé-Massuyès1,2 and R. Pons1,2 and P. Ribot1,3 and Y. Pencole1,2 and C. Jauberthie1,3
1
CNRS, LAAS, 7 avenue du colonel Roche, F-31400 Toulouse, France
2
Univ de Toulouse, LAAS, F-31400 Toulouse, France
3
Univ de Toulouse, Université Paul Sabatier, LAAS, F-31062 Toulouse, France
e-mail: louise@laas.fr, renaud.pons@gmail.com, pribot,ypencole,cjaubert@laas.fr
Abstract stress conditions. These uncertainties are commonly taken
into account through appropriate assumptions about noise
Condition-based maintenance is recognized as a and model error distributions, which are difficult to acquire.
better health management strategy than regularly An alternative approach is to frame the problem in a set-
planned inspections as used nowadays by most membership framework and make use of recent advances in
companies. In practice, it is however difficult to the field of interval analysis and interval constraint propaga-
implement because it means being able to predict tion.
the time to go before a failure occurs. This predic- This paper demonstrates the applicability of interval-
tion relies on knowing the current health status of based tools — briefly introduced in Section 2 — in inte-
the system’s components and on predicting how grated health management architectures, providing an in-
components age. This paper demonstrates the teresting alternative to the standard statistical approach [5;
applicability of interval-based tools in integrated 6]. It proposes a two stages set-membership (SM) condition-
health management architectures, hence propos- based monitoring method whose principle is presented in
ing an alternative to the standard statistical ap- Section 3. The first stage is diagnosis that provides an esti-
proach1 . mation of the system’s health status. It takes the form of SM
Keywords: Interval analysis, diagnosis, prognosis. parameter estimation using Focused Recursive Partitioning
(FRP) and is the subject of Section 4. The second stage con-
cerns prognosis in the form of the estimation of the remain-
1 Introduction ing system’s lifespan. It is based on the use of a damaging
Nowadays system’s availability is a key ingredient of eco- table and is detailed in Section 5. The case study of a shock
nomical competitiveness. A typical example is the civil absorber is used to illustrate the method and is presented in
aircraft industry for which the unavailability of passenger Section 6 before the concluding Section 7.
carriers generate great costs and considerable economical
losses. Technical inspections are generally planned on reg- 2 Interval analysis
ular bases. If a component fails and needs to be replaced Interval analysis was originally introduced to obtain guar-
between two successive inspections, the plane is taken to a anteed results from floating point algorithms [7] and it was
standstill and the company has to re-schedule the aircraft then extended to validated numerics [8]. A guaranteed re-
fleet, implying money loss during this unplanned immo- sult first means that the solution set encloses the actual so-
bilization. This is why condition-based maintenance is a lution. It also means that the algorithm is able to conclude
preferable strategy that means predicting at inspection time about the existence or not of a solution in limited time or
the time to go before a new failure occurs. In this case the number of iterations [9].
aircraft company can replace the part whose failure is es-
timated during the current inspection and then prevent an 2.1 Interval
extra immobilization of the plane. This strategy not only A real interval x = [x, x] is a closed and connected subset
saves a lot of money but also increases reliability and safety. of R where x and x represent the lower and upper bound of
Integrated systems health management architectures per- x, respectively. x and x are real numbers. The width of an
forming condition-based monitoring naturally couple fault interval x is defined by w(x) = x − x, and its midpoint by
diagnosis and prognosis mechanisms [1; 2; 3; 4]. Diagnosis mid(x) = (x + x)/2. If w(x) = 0, then x is degenerated
is used to assess the current state of the system and is used and reduced to a real number. x is defined as positive (resp.
to initialize a prediction mechanism based on ageing mod- negative), i.e. x ≥ 0 (resp. x ≤ 0), if x ≥ 0 (resp. x ≤ 0).
els that aims to estimate the remaining useful life (RUL). The set of all real intervals of R is denoted IR. Two inter-
In the prognosis process several sources of uncertainty can vals x1 and x2 are equal if and only if x1 = x2 and x1 = x2 .
be identified, in particular the ageing models and the future Real arithmetic operations have been extended to intervals
[8]:
1
This work was supported by the CORALIE Project IA 2012–
◦ ∈ {+, −, ∗, /}, x1 ◦ x2 = {x ◦ y | x ∈ x1 , y ∈ x2 }.
01–06 of the Council for Research in French Civil Aeronautics
(CORAC), WP1 "Contrôle Santé", and by the French National Re- An interval vector or box [x] is a vector with interval com-
search Agency (ANR) in the framework of the project ANR-11- ponents. An interval matrix is a matrix with interval com-
INSE-006 (MAGIC-SPS). ponents. The set of n−dimensional real interval vectors is
83
Proceedings of the 26th International Workshop on Principles of Diagnosis
denoted by IRn and the set of n × m real interval matrices Ageing model
is denoted by IRn×m . The width w(.) of an interval vector Behavioural Damaging
model Σ model Δ
(or of an interval matrix) is the maximum of the widths of
its interval components. The midpoint mid(.) of an interval
vector (resp. an interval matrix) is a vector (resp. a matrix)
Diagnosis
composed of the midpoints of its interval components. Prognosis RUL
Classical operations for interval vectors (resp. interval
matrices) are direct extensions of the same operations for
real vectors (resp. real matrices) [8].
System
2.2 Inclusion function
Given x a box of IRn and a function f from IRn to IRm , an Figure 1 – Health management architecture.
inclusion function of f aims at getting a box containing the
image of x by f . The range of f over x is given by:
where
f (x) = {f (ν) | ν ∈ x}, x(t) is the state vector of the system of dimension nx ,
u(t) is the input vector of dimension nu ,
where ν is a real vector of the same dimension as x. Then,
y(t) is the output vector of dimension ny ,
the interval function [f ] from IRn to IRm is an inclusion
θ(t) is the parameters vector of dimension nθ .
function for f if:
Σ represents the system’s nominal behavior as it is sup-
∀x ∈ IRn , f (x) ⊂ [f ](x). posed to act when its parameters have not yet suffered any
ageing.
An inclusion function of f can be obtained by replacing The damaging model ∆ represents the dynamics of the
each occurrence of a real variable by its corresponding in- behavioral model parameters. It is described by the dy-
terval and by replacing each standard function by its interval namic state equation (2) where the equation states are the
evaluation. Such a function is called the natural inclusion system parameters. The equation models how the parame-
function. A function f generally has several inclusion func- ters evolve over time because of the wearing, leakage, etc.:
tions, which depend on the syntax of f .
∆ : θ̇(t) = g(t, θ(t), w, x(t)) (2)
2.3 Notations
where w is a wearing parameter vector of dimension nθ .
Throughout the paper and unless explicitly mentioned, vari-
ables are assumed to take values in IRd , where d is the di- 3.2 Unit Cycles
mension of the variable. Exception is made for overlined Predicting the evolution of the system’s behavior requires
and underlined variables that are assumed to take values in to know a priori how the system will be solicited either by
Rd , where d is the dimension of the variable. Bold sym- the control system or by external causes (e.g. environmen-
bols are used to denote multi-dimensional variables (vector tal conditions, temperature, humidity, etc.) This knowledge
or matrices). is generally difficult to obtain. In our approach, we make
assumptions about the future solicitations of the system by
3 Principle of the Set-Membership Health determining the most usual way the system is intended to be
Management Method used and we define the notion of unit cycles. A unit cycle C
is defined as a solicitation that repeats in time and that leads
3.1 Method Architecture to a behavioral sequence that is known to impact system’s
The architecture of the preventive maintenance method is ageing.
shown in Fig. 1. The method relies on two modules: For example, in the case of a pneumatic valve from the
• A diagnosis module that uses the system measured in- Space Shuttle cryogenic refueling system, [2] defines a unit
puts and outputs to compute an estimation of the sys- cycle as the opening of the valve, the filling of the tank and
tem’s health status; this is performed by estimating the the valve closing when the tank is full. In the case of an
value of the system parameter vector θ by the means of aircraft, an unit cycle may be chosen to be a flight: it starts
a behavioral Model Σ of the system. with the plane take-off, a cruising stage and landing.
One may simultaneously use unit cycles at different time
• A prognosis module that predicts the parameter evo- scales, depending on the dynamics of the system and its sub-
lution over time by using a Damaging Model ∆ and systems. As an example, one may define a “global” unit cy-
computes the Remaining Useful Life or RUL of the un- cle for a bus as being the journey from the starting station to
derlying subsystems. the terminus, and another unit cycle for the subsystem “bus
The global model representing the progressive evolution doors” as being the opening and the closing of the doors at
of the system over time, e.g. the Ageing Model, is obtained each station.
by putting together the behavioral model Σ and the dam-
aging model ∆. The behavioral model takes the form of a 4 SM Diagnosis
state space model, i.e. a state equation modeling the dynam-
ics of the system state vector, and an observation equation Diagnosis is achieved through SM parameter estimation.
that links the state variables to the observed variables: This problem assumes that measured outputs ym (ti ) gener-
ated by the real system on a time horizon ti = t0 , . . . , tH of
ẋ(t) = f (t, x(t), θ(t), u(t)) length H × δ, where δ is the sampling period, are corrupted
Σ : (1) by bounded-error terms that may originate from the system
y(t) = h(t, x(t), θ(t), u(t))
84
Proceedings of the 26th International Workshop on Principles of Diagnosis
parameters varying within specified bounds, bounded noise, and hence rejected.
or sensor precision. The ym (ti )’s are hence interval vectors
of IRny . The SM parameter estimation problem for the sys- If ŷ(ti )j ⊆ ym (ti ), ∀ti ∈ {t0 , . . . , tH }, (5)
tem Σ is formulated as finding the set Θ ⊆ Rnθ of real pa- then [θ]j is a parameter vector for which the estimation is
rameter vectors such that the arising outputs y(ti , θ) ∈ Rny consistent with the measurements. The box is added to the
hit all the output data sets, i.e.: list of the solution parameter boxes:
θ ∈ Θ ⇔ y(ti , θ) ∈ ym (ti ), ∀ti ∈ {t0 , . . . , tH }.
P = P ∪ [θ]j . (6)
Θ is called the feasible parameter set (FPS). SM parame-
ter estimation problems are generally solved with a branch- If none of the two previous conditions is true, i.e.:
and-bound algorithm like SIVIA [10] that enumerates can- [ŷ(ti )]j ∩ [ym (ti )] 6= ∅, ∀ti ∈ {t0 , . . . , tH }, (7)
didate box solutions thanks to a rooted tree and assumes the
full parameter space as the root. At every node, the set of it means that the parameter box [θ]j is undetermined and
ŷ(ti ), ti = t0 , . . . , tH , arising from the considered box pa- that it partially contains solutions. The box is also added to
rameter vector [θ]∗ , i.e. solution of Σ for any real θ ∈ [θ]∗ , P. Two different labels allow us to keep track of the boxes
is checked for consistency against the measurements and la- that are feasible or undetermined. The convex union of these
belled feasible, unfeasible or undetermined. Unfeasible can- boxes provides the estimation θ̂ that encloses the feasible
didates are rejected while undetermined candidates are split parameter set Θ, i.e. θ̂ ⊇ Θ.
and checked in turn until the set precision of the candidate The quality of the enclosure depends on the size of the
solutions is below a given threshold ε provided by the user. boxes of the partition, in other words on the partition preci-
Such algorithms return an overestimation of the FPS sion. A way to improve the enclosure is to proceed with a
given by the convex union of the candidates that have been partition of the obtained solution θ̂ and run another round of
labelled feasible and undetermined [11]. Interestingly, the consistency tests over the new boxes, and so on recursively.
convex union may consist of one set or more, which means The process of iterating the partition ends when the gain in
that the systems does not need to be identifiable in the clas- precision is low with respect to the SM integration and con-
sical sense [12]. sistency tests computational cost. The method is detailed
When considering a SIVIA-based algorithm for dynam- for a one dimension parameter vector in the next section.
ical systems like Σ, a critical step is the determination of The estimation precision ω(P ) obtained for a given par-
the inclusion function for the state vector x̂(ti ) at instants tition P can be evaluated by the following percentage:
ti = t0 , . . . , tH , arising from a given candidate parameter
vector [θ]∗ , from which the [ŷ(ti )], ti = t0 , . . . , tH can ω(P ) = mid(θ̂) ./ mid(θ̂) + w(θ̂)/2 (8)
be computed using the observation equation of Σ. This
step relies on set-membership integration for which we have where ./ denotes the division of two vectors term by term.
chosen the interval Taylor series integration scheme imple-
Given two partitions Pi and Pj , one can evaluate the pre-
mented in the VNODE-LP solver [13]. Although quite well
cision gain as:
optimized [14], it is well-known that this method is compu-
tationally stable only for [θ]∗ of very small size. SIVIA-like G(Pj /Pi ) = ω(Pj )./ω(Pi ). (9)
parameter estimation algorithms are hence particularly in-
efficient as they enumerate candidate parameter subspaces 4.2 Parameter search space
starting with the full parameter space. This is why we pro- The domain value of the parameter vector θ is given by
pose the FRP schema presented in the next section. Ω(θ) = [.inf(θ bol , θ eol ), .sup(θ bol , θ eol )], where .inf and
4.1 Principle of FRP-based SM Parameter .sup denote the operators inf and sup applied term by term
and θ bol and θ eol are the real vectors whose components
Estimation are given by θk,bol and θk,eol for each parameter θk , k =
The principle of the FRP method is based on partitioning 1, . . . , nθ :
the parameter search space S(θ). Each part of the partition
represents a candidate parameter vector [θ]j for which SM • θk,bol (bol: “beginning-of-life”) is the factory setting
integration of the state equation of Σ provides a conservative defined by the design specification,
numerical enclosure x̂(ti )j , ti = t0 , . . . , tH . The output • θk,eol (eol: “end-of-life”) is the maximal/minimal ad-
vector can now be estimated as: missible value, i.e. the value above/below which the
ŷ(ti )j = h (t, x̂j , [θ]j , um ) , ∀ti ∈ {t0 , . . . , tH } (3) component is considered to have failed and the func-
tion is no longer guaranteed.
We then keep track of the parameters vectors for which
During the system’s life, the impact of ageing results in the
the ŷ(ti )j ’s are consistent with the measurements, for all
parameter vector value evolving in Ω(θ). Depending on the
ti = t0 , . . . , tH , i.e. the unfeasible ones are discarded. Com-
impact of ageing, its value may decrease or increase with
puting the convex hull then provides us with a minimal and
time:
maximal value for the admissible parameter vectors.
θ(ti ) = αθ(tj ), ti ≥ tj (10)
The consistency test is defined as testing the intersection
of the estimated output vector with the measurements: where α is an nθ dimensional real vector whose compo-
nents αk are greater or lower than 1 depending on the impact
If ∃ ti ∈ {t0 , . . . , tH } s.t. ŷ(ti )j ∩ ym (ti ) = ∅, (4)
of ageing on the change direction of the parameter.
then there is no consistency between the estimation and the We assume a health management strategy, which means
measured input um (t) and output ym (t) with the tested pa- that diagnosis (and prognosis) is performed according to a
rameter vector [θ]j . The parameters box [θ]j is unfeasible given inspection planning at some chronologically ordered
85
Proceedings of the 26th International Workshop on Principles of Diagnosis
times of the system’s life T0 , . . . , TF . The parameter vec- We proceed as above
h for
i each interval of P2 in order to re-
tor search space depends on the time and is hence denoted fine the bounds of θ̂, θ̂ and find a more precise enclosure
by S(θ)Ti = [S(θ)Ti , S(θ)Ti ], where Ti ∈ {T0 , . . . , TF }.
Initially, for the first inspection time T0 , we set S(θ)T0 = of feasible parameter solution (see Fig. 3).
Ω(θ). Diagnosis then returns the estimated parameter value
θ̂(T0 ). For the next inspection time, S(θ) is updated by tak- P1
ing the parameter value estimation into account as follows:
• if αk > 1, θk,bol is replaced by inf(θˆk (T0 ), θˆk (T0 )), P2
• if αk < 1, θk,bol is replaced by sup(θˆk (T0 ), θˆk (T0 )), Figure 3 – Partition P2 and test results for this partition.
and one of the bounds of the components of S(θ)T1 remains
equal to θk,eol . We iterate the process until the precision gain
In the general case, when considering the inspection time G(Pi+1 /Pi ) is greater then a given threshold, as it is
Ti , S(θ)Ti is hence obtained with θ̂(Ti−1 ) as follows: shown in Fig. 4.
• if αk > 1, inf(θˆk (Ti−1 ), θˆk (Ti−1 )) is replaced by
P1
inf(θ̂(Ti ), θ̂(Ti )),
• if αk < 1, sup(θˆk (Ti−1 ), θˆk (Ti−1 )) is replaced by
P2
sup(θˆk (Ti ), θˆk (Ti )).
4.3 FRP Parameter Estimation for a Single P3
Parameter
Figure 4 – Test results for partition P3 .
In this section, we consider one single parameter θ whose
evolution is monotonically increasing. As an example, let’s
state that θ is a bearing friction coefficient that grows with Remarks
the bearing wearing and the clogging of the environment. In The method can be easily generalized to a system whose
the general case, this kind of knowledge must be brought by parameter vector has dimension nθ > 1. The computing
an expert of the system and/or the manufacturer. costPis proportional to the number of boxes that are tested,
Let us consider the first inspection time and the initial nP
i.e. i=1 1/(Pi ), where nP is the number of partitions.
search space S(θ)T0 given by the domain value of the pa- Let’s notice that the partition may be non-regular. For ex-
rameter Ω(θ) = [θbol , θeol ]. The search space is partitioned ample, for a slowly ageing parameter, one may choose small
into boxes, in our case intervals (cf. Fig. 2). boxes for the values of θ that are close to θbol and larger
The dynamic equation of Σ is integrated on the time win- ones for the values close to θeol . The result is guaranteed
dow ti = t0 , . . . , tH , where tH = T0 , as many times as even if the partition has not been properly chosen or if the
the number of intervals in the partition P1 . The number parameter has evolved in a non expected way, although the
of intervals is defined by the partition factor (P1 ), which computation cost may be higher.
equals 1/15 in our example (cf. Fig. 2). We start with The convex union provides a poor result if the set of ad-
[θ]1 = [θbol , θbol + pw], where pw = (P1 )w(S(θ)T0 ) is missible values is made of several mutually disjoint con-
the width of the partition intervals, then proceed with the nected sets, as shown in Fig. 5. The algorithm may test
subsequent intervals [θ]j . For each interval, we get an esti- some boxes that have already been rejected by the tests of
mation of the state vector at times ti = t0 , . . . , tH , denoted the previous partition. This drawback could be addressed
as x̂(t0 . . . tN )j , and obtain ŷ(t0 . . . tN )j thanks to the ob- by defining the solution as a list of boxes whose labels (un-
servation equation (1). This latter is tested for consistency feasible, feasible, or undetermined) are inherited by the next
against the measurements ym (t0 . . . tN ). partition boxes.
Depending on the output of the tests (4), (5), and (7), the
parameter interval [θ] is rejected or added to the solution as
feasible or undetermined (red-colored, green-colored, and P
yellow-colored parts, respectively, in Fig. 2). Solution
Figure 5 – The returned solution is the convex hull of mutu-
P1 ally disjoint connected intervals.
Figure 2 – Partition P1 and test results for this partition.
5 SM prognosis
The convex union of feasible and undetermined
h i intervals The prognosis phase consists is calculating the number of
provides a guaranteed estimation θ̂ = θ̂, θ̂ of the admis- cycles remaining before anomaly, which is also called the
sible values for hθ. Wei iterate the process by creating a new Remaining Useful Life or RUL. To optimally adapt this cal-
culation to the system’s life requires the knowledge of the
partition P2 of θ̂, θ̂ with a precision (P2 ) = 1/10 (cf. health status of the system at the current time, which was
Fig. 3). the topic of Section 4.
86
Proceedings of the 26th International Workshop on Principles of Diagnosis
5.1 Component degradation partitioned into NΠ = Πnk=1θ
Nk possible boxes that must be
The global model (Σ + ∆) assumes that the parameters of fed as input to D. Let us for instance consider a two param-
the behavior model Σ given by (1) evolve in time, and that eters vector and its beginning-of-life and end-of-life values
their evolution is represented by the degradation model ∆ as follows:
given by the dynamic equation (2) that is recalled below: θ 1 4 3
θ = 1 , θ bol = , θ eol = and N = , (12)
∆ : θ̇(t) = g(t, θ(t), w, x(t)). θ2 1 9 2
∆ provides the dynamics of the parameter vector as a func- then, if we select the partition landmarks as {5} for θ1 and
tion of the state of the system x(t) and of a degradation {2, 3} for θ2 2 , D must be run for the following 6 box values:
parameter vector w that allows one to tune the degradation
for each of the considered parameters. [1, 2] [2, 3] [3, 4]
[θ]1 = , [θ]2 = , [θ]3 = ,
The global model (Σ + ∆), in the form of a dynamic [1, 5] [1, 5] [1, 5]
model with varying parameters, cannot be directly inte-
[1, 2] [2, 3] [3, 4]
grated by VNODE-LP. An original method, coupling the [θ]4 = , [θ]5 = , [θ]6 = . (13)
[5, 9] [5, 9] [5, 9]
two models Σ and ∆ iteratively is proposed in the following.
The method is illustrated by Fig. 6 and used to determine For each of these box values taken as input for cycle i, i.e.
the degradation suffered by each parameter during one unit θ i = [θ]l , l = 1, . . . , 6, D returns the (box) value θ i+1 after
cycle as defined in Section 3.2. one unit cycle. This computation is then projected on each
Let us denote uC (t), t ∈ [τ, τ + dC ], the system input dimension to obtain a set of nθ tables, Dθk , k = 1, . . . , nθ ,
stress during one unit cycle C . As shown in Fig. 6, the that provide the degradation of each individual parameter θk
following steps are iteratively executed, every iteration cor- after one unit cycle.
responding to a computation step given by the sampling pe-
riod δ: 5.2 RUL determination
1. The normal behavior model Σ is used first with input The RUL, understood as a RUL for the whole system, can
u(t) = uC (τ ) to compute the state x(τ ) and the output now be determined by computing the number of cycles that
y(τ ); are necessary for the parameters to reach the threshold defin-
2. The parameters are updated with the degradation ing the end-of-life (cf. Fig. 7).
model ∆ using the value of the state determined pre-
viously, i.e. θ(τ ) is computed;
Diagnosis at
3. The parameters of the behavior model Σ are updated inspection
No
with θ(τ ); time Tk
4. The next stress input value uC (τ + δ) is considered, Yes
and so on until the end of the cycle, i.e. until the last RUL = i
value of the cycle uC (τ + dC ) is reached.
Figure 7 – RUL computation
Σ
For the cycle i = 0, θ 0 is initialized with θ̂, which is
the result of the parameter estimation computed by the di-
agnosis engine. θ̂ is given as input to D, which returns
D(θ̂) = θ 1 . i is incremented by 1 and θ 1 is given as input
Δ
to D and so on until the set-membership test θ i θ eol is
achieved, which provides the stopping condition. This test
Figure 6 – Computation of the degradation parameters dur- may take several forms as explained in Section 5.3. If the
ing one unit cycle. test is true, then the index i is the number of cycles required
to reach the degradation threshold, so RUL = i.
The above algorithm defines the function: For a given cycle i, the box value θ i that must be given
as input to D is not necessarily among the values [θ]l , l =
D : IRnθ → IRnθ (11)
1, . . . , NΠ , of the partition. We propose to compute θ i+1
where nθ is the number of parameters of the system. Let’s by assuming that the mapping between θ i and θ i+1 is linear
assume the cycle i, then D maps θ i into D(θ i ) = θ i+1 , in every domain l of the partition. Considering p ∈ Rnθ ,
which is the value of θ after one unit cycle. D(p) is approximated as follows:
D is nonlinear. Thus the value of the parameter vector
∀θ ∈ [θ]l , D(θ) ≈ a θ + b, l=1, . . . , NΠ (14)
after one cycle θ i+1 depends on the initial value θ i . Indeed,
we know that a system generally degrades in a nonlinear where a= w(D([θ]l ))./w([θ]l ), b= D([θ]l ) − a[θ]l , and
fashion. We must hence compute θ i+1 for all possible val- is the product of two vectors term by term.
ues of the parameter vector θ i . i
Equation (14) is applied to θ and θ i to obtain an approx-
For this purpose, the domain value Ω(θk ) of each param- i
imation of D(θ ).
eter θk is partitioned into Nk intervals. Nk is chosen suffi-
ciently large to reduce non conservatism of the interval func- 2
Notice that the intervals issued from the partitioning are not
tion D. The domain value of the parameter vector θ is hence required to be of equal length.
87
Proceedings of the 26th International Workshop on Principles of Diagnosis
5.3 Set-membership test for the RUL • if ζ = 0, then the answer is a sinusoid;
The set-membership test implemented with the order rela- • if 0 < ζ < 1, then the answer is a damped sinusoid;
tion may take several forms. For instance, if the test
• if ζ ≥ 1, then the answer is a decreasing exponential.
θ i θ eol is interpreted as:
The state model is given by the equation:
∃k ∈ {1, . . . , nθ } |
i 0 1 0
θk ≥ θk,eol if αk > 1 or θik ≤ θk,eol if αk < 1, (15) Ẋ(t) = −k −c X(t) +
U (t)
1
then it means that the bound of the interval value of at least m m (19)
one parameter θk is above or below its end-of-life threshold Y (t) = 1 0 X(t)
value θk,eol . The RUL is then qualified as the “worst case 0 1
RUL”, which means that the RUL indicates the earliest cycle with X(t) = [x(t), ẋ(t)]T , and the transfer function is:
at which the system may fail.
One can also test whether the value higher bound of one X(p) 1
of the parameters is higher than its end-of-life threshold, that = 2 c k
. (20)
U (p) p +m p+ m
is to say:
An example of bounded error step response obtained with
∃k ∈ {1, . . . , nθ } |
VNODE-LP with a sampling parameter δ = 0.1 s, c = 1,
i
θik ≥ θk,eol if αk > 1 or θk ≤ θk,eol if αk < 1. (16) m = 2 and k = [3, 9 ; 4, 1] is shown in Fig. 9a. There,
The RUL then represents the cycle at which it is certain that ζ ' 0.177 and the step response is a damped sinusoid. Be-
the system will fail. cause k is assumed to have an uncertain value bounded by
It is obviously possible to combine these different tests an interval, the outputs are in the form of envelops.
applied to the different individual parameters depending on 6.2 Unit cycle
their criticality.
In the case study, a unit cycle is defined by the application
6 Case study of a power unit for a determined time. The force is applied
at time t0 +5s, where t0 is the cycle starting time. The force
6.1 Presentation lasts 20s and cancels at t0 + 25s as shown by the red curve
The case study is a shock absorber that consists of a moving of Fig. 9b. The cycle ends at t0 + 50s.
mass connected to a fixed point via a spring and a damper as Fig. 9b presents the system’s response for a spring con-
illustrated by Fig. 8. The movement of the mass takes place stant k = [3.9, 4.1] N/m, a mass m = 2 kg, a damping co-
in the horizontal plane in order to eliminate the forces due efficient c = 10 Ns/m, and initial speed and position equal
to gravity. Aerodynamic friction forces are neglected. to zero. The response is a decreasing exponential.
x 6.3 Degradation model
k The degradation model chosen is the ageing of the damper
m cylinder. It is represented by a reduction of the damping
c
coefficient proportional to the velocity of the mass [15]:
ċ = β ẋ, β < 0. (21)
Figure 8 – Spring and damper system The more the spring is used, the weaker it becomes, charac-
terized by the change in the damping coefficient.
The Newton’s second law is written as:
6.4 Diagnosis
m~a = ΣF~ = F~r + F~c + ~u (17)
The FRP parameter estimation method presented in Section
where m is the mass, ~a is the acceleration, F~k is the spring 4 has been used with the measures shown in Fig. 9c. These
biasing force, F~c is the friction force exerted by the damper measures were obtained for
and ~u is the force applied on the mass. Expressing the forces " # " #
c 5
and the acceleration as a function of the position of the mass
θ= k = 4 (22)
x(t), we get: m 2
c k
ẍ(t) + ẋ(t) + x(t) = u(t) (18) The goal is to estimate the damping coefficient c and the
m m
where k is the spring stiffness constant (N/m), m is the mo- stiffness constant k. The search space is defined by the inter-
bile mass (kg), and c is the damping coefficient (Ns/m). (18) val [4 9] for c and [3.5, 9] for k. The value of m is assumed
is a second order ODE. Let us rewrite to be known m = 2. Using the notation introduced above,
c k we have: " # " #
= 2ζω 0 and = ω 20 4 9
m m θ bol = 3.5 , θ eol = 9 (23)
and we get 2 2
r
k c The partition P1 is achieved with a precision (P1 ) = 1/10
ω0 = and ζ = √ .
m 2 km for the two parameters to be estimated c and k. Fig. 10
The impulse response of such system depends on the value presents two examples of prediction results with two param-
T
of ζ: eter boxes of P1 : [θ]i = [4, 1, 4, 2], [4, 7, 4, 8], 2 on the
88
Proceedings of the 26th International Workshop on Principles of Diagnosis
(a) Step response for ζ = 0.177. (b) Unit cycle for the case study. (c) Measured input and output.
Figure 9 – Cases study simulation and data plots.
T
left and [θ]j = [5, 5, 1], [4, 4, 1], 2 on the right. On the 6.5 RUL computation
left figure, one can see that there is no intersection between In this section we apply the set-membership method de-
the estimate and the measurement for the position, hence the scribed in Section 5.2 to compute the RUL for the damping
box used for the simulation is rejected. On the right, there coefficient c.
is an intersection between the measurement and the estima- The damper is assumed to fail when c ≤ ceol = 2. The
tion for all time points, but the estimate is not included in degradation model (21) with β = −0, 13 allows us to deter-
the measure envelop, hence the parameter box is considered mine the degradation table Dc for the parameter c for a unit
undetermined. cycle:
ci D(ci ) = ci+1 ci D(ci ) = ci+1
[9, 10] [8.917, 9.977] [4, 5] [3.874, 4.977]
Dc = [8, 9] [7.911, 8.978] [3, 4] [2.863, 3.973]
[7, 8] [6.898, 7.979] [2, 3] [1.721, 2.97]
[6, 7] [5.814, 6.982] [1, 2] [0, 1.98]
[5, 6] [4.859, 5.979] [0, 1] [0, 0.9755]
(25)
Figure 10 – Estimation results with a rejected parameter box After proceeding to the linear interpolation given by (14),
(left) and an indetermined box (right) the graphical representation of ci+1 as a function of ci is
given by Fig. 12.
The results for partition P1 are presented in Fig. 11a and
we obtain a first estimation for θ:
" #
[4, 1, 5.8]
θ̂ = [3, 7, 4, 2] .
2
The estimation precision for partition P1 is given by:
" #
0.85
ω(P1 ) = mid(θ̂) ./( mid(θ̂) + w(θ̂)/2) = 0.94 (24)
1
The first estimation for θ is used as the search space for
partition P2 , whose precision is increased by a factor of 10, Figure 12 – Approximated degradation of the damping co-
i.e. (P2 ) = 0, 1. The obtained estimation results are shown efficient c
in Fig. 11b.
The estimation is refined as: The number of elements of the partition has been chosen
"
[4, 51, 5, 57]
# relatively small to better illustrate the method. In a real sit-
θ̂ = [3, 85, 4, 14] . uation, this number should be high in order to obtain less
2 conservative predictions.
The value of c has been previously estimated and is
The precision is now ω(P2 ) = [0.9, 0.96, 1]T , and the ĉ = [4.548, 5.526].
precision gain is G(P2 /P1 ) = [0.056, 0.025, 0]T . The val-
ues for the gain indicate that partitioning a third time might The graph of Fig. 12 allows us to approximate the predicted
be quite inefficient. To confirm this fact, let us perform a value after one unit cycle:
third partition P3 , whose precision is increased by a factor D(ĉ) = c1 = [4.4787, 5.4481].
of 5, i.e. = 0.02 (cf. Fig. 11c). The new estimation for θ is
T The next iteration of the algorithm allows us to compute c2 ,
θ̂ = [4.548, 5.526], [3.872, 4.132], 2 , and the precision
etc. After 30 iterations, we obtain c30 = [1.7665, 3.4235].
gain is G(P3 /P2 ) = [0.0073, 0.0036, 0]T . As expected, the
gain is quite negligible with respect to the computation time 3
The coefficient β has been chosen arbitrarily to illustrate the
increase. approach; it does not represent the real ageing of a damper.
89
Proceedings of the 26th International Workshop on Principles of Diagnosis
(a) Partition P1 (b) Partition P2 (c) Partition P3
Figure 11 – Partitions and estimation results (red, yellow and green boxes are resp. rejected, undetermined, accepted param-
eters values).
Since c30 < ceol , we get RU L = 30 cycles. After the 44th [5] D. Gucik-Derigny, R. Outbib, and M. Ouladsine. Es-
iteration, we get c44 = [0.037591, 1.985928]. We then have timation of damage behaviour for model-based prog-
c44 < ceol and hence RU L = 44 cycles. The RUL of the nostic. In Fault Detection, Supervision and Safety of
damper is hence given by: Technical Processes, pages 1444–1449, 2009.
[6] X. Guan, Y. Liu, R. Jha, A. Saxena, J. Celaya, and
RU L = [30, 44] cycles.
K. Geobel. Comparison of two probabilistic fatigue
damage assessment approaches using prognostic per-
7 Conclusion formance metrics. International Journal of Prognos-
This paper addresses the condition-based monitoring and tics and Health Management, 1(005), 2011.
prognostic problems with a new focus that trades the tra- [7] R.E. Moore. Automatic error analysis in digital com-
ditional statistical approach by an error-bounded approach. putation. Technical report LMSD-48421, Lockheed
It proposes a two stages method whose principle is to first Missiles and Space Co, Palo Alto, CA, 1959.
determine the health status of the system and then use this [8] R.E. Moore. Interval Analysis. Prentice-Hall, Engle-
result to compute the RUL of the system. This study uses
wood Cliffs, 1966.
advanced interval analysis tools to obtain guaranteed results
in the form of interval bounds for the RUL. [9] L. Jaulin, M. Kieffer, O. Didrit, and E. Walter. Ap-
The results for the case study demonstrate the feasibility plied Interval Analysis, with examples in parameter
of the approach. The next step is to adapt the FRP-based and state estimation, Robust control and robotics.
SM parameter estimation algorithm in order to output a list Springer, Londres, 2001.
of boxes instead of a single box given by the convex hull [10] L. Jaulin and E. Walter. Set inversion via interval anal-
of the boxes. The convex hull is indeed a very conservative ysis for nonlinear bounded-error estimation. Automat-
approximation when the solution set is not convex. ica, 29:1053–1064, 1993.
The second stream of work is to consider contextual con-
[11] L. Jaulin, M. Kieffer, O. Didrit, and E. Walter. Ap-
ditions and their associated uncertainties. Environmental
conditions, like weather, different usage, etc. may indeed plied Interval Analysis, with examples in parameter
significantly affect the stress input and prognostics results. and state estimation, Robust control and robotics.
Springer, Londres, 2001.
[12] C. Jauberthie, N. Verdière, and L. Travé-Massuyès.
References
Fault detection and identification relying on set-
[1] Indranil Roychoudhury and Matthew Daigle. An inte- membership identifiability. Annual Reviews in Con-
grated model-based diagnostic and prognostic frame- trol, 37:129–136, 2013.
work. In Proceedings of the 22nd International Work- [13] N. Nedialkov. VNODE-LP, a validated solver for ini-
shop on Principle of Diagnosis (DX’11). Murnau,
tial value problems for ordinary differential equations.
Germany, 2011.
[14] R. J. Lohner. Enclosing the solutions of ordinary initial
[2] Matthew J Daigle and Kai Goebel. A model-based
and boundary value problems. In E. W. Kaucher, U. W.
prognostics approach applied to pneumatic valves. In-
Kulisch, and C. Ullrich, editors, Computer Arithmetic:
ternational Journal of Prognostics and Health Man-
Scientific Computation and Programming Languages,
agement, 2:84, 2011.
pages 255–286. Wiley-Teubner, Stuttgart, 1987.
[3] J. Luo, K. R Pattipati, L. Qiao, and S. Chigusa. Model- [15] Ian M Hutchings. Tribology: friction and wear of
based prognostic techniques applied to a suspension engineering materials. Butterworth-Heinemann Ltd,
system. Systems, Man and Cybernetics, Part A: Sys- 1992.
tems and Humans, IEEE Transactions on, 38(5):1156–
1168, 2008.
[4] Q. Gaudel, E. Chanthery, and P. Ribot. Hybrid par-
ticle petri nets for systems health monitoring under
uncertainty. International Journal of Prognostics and
Health Management, 6(022), 2015.
90