Condition-based maintenance is recognized as a better health management strategy than regularly planned inspections as used nowadays by most companies. In practice, it is however difficult to implement because it means being able to predict the time to go before a failure occurs. This prediction relies on knowing the current health status of the system's components and on predicting how components age. This paper demonstrates the applicability of interval-based tools in integrated health management architectures, hence proposing an alternative to the standard statistical approach1.
Nowadays system’s availability is a key ingredient of economical competitiveness. A typical example is the civil aircraft industry for which the unavailability of passenger carriers generate great costs and considerable economical losses. Technical inspections are generally planned on regular bases. If a component fails and needs to be replaced between two successive inspections, the plane is taken to a standstill and the company has to re-schedule the aircraft fleet, implying money loss during this unplanned immobilization. This is why condition-based maintenance is a preferable strategy that means predicting at inspection time the time to go before a new failure occurs. In this case the aircraft company can replace the part whose failure is estimated during the current inspection and then prevent an extra immobilization of the plane. This strategy not only saves a lot of money but also increases reliability and safety.
Integrated systems health management architectures performing condition-based monitoring naturally couple fault diagnosis and prognosis mechanisms [1; 2; 3; 4]. Diagnosis is used to assess the current state of the system and is used to initialize a prediction mechanism based on ageing models that aims to estimate the remaining useful life (RUL). In the prognosis process several sources of uncertainty can be identified, in particular the ageing models and the future stress conditions. These uncertainties are commonly taken into account through appropriate assumptions about noise and model error distributions, which are difficult to acquire. An alternative approach is to frame the problem in a setmembership framework and make use of recent advances in the field of interval analysis and interval constraint propagation.
This paper demonstrates the applicability of intervalbased tools — briefly introduced in Section 2 — in integrated health management architectures, providing an interesting alternative to the standard statistical approach [5; 6]. It proposes a two stages set-membership (SM) conditionbased monitoring method whose principle is presented in Section 3. The first stage is diagnosis that provides an estimation of the system’s health status. It takes the form of SM parameter estimation using Focused Recursive Partitioning (FRP) and is the subject of Section 4. The second stage concerns prognosis in the form of the estimation of the remaining system’s lifespan. It is based on the use of a damaging table and is detailed in Section 5. The case study of a shock absorber is used to illustrate the method and is presented in Section 6 before the concluding Section 7. 2
Interval analysis was originally introduced to obtain
guaranteed results from floating point algorithms [7] and it was
then extended to validated numerics [8]. A guaranteed
result first means that the solution set encloses the actual
solution. It also means that the algorithm is able to conclude
about the existence or not of a solution in limited time or
number of iterations [
A real interval x = [x, x] is a closed and connected subset of R where x and x represent the lower and upper bound of x, respectively. x and x are real numbers. The width of an interval x is defined byw(x) = x − x, and its midpoint by mid(x) = (x + x)/2. If w(x) = 0, then x is degenerated and reduced to a real number. x is defined as positive (resp. negative), i.e. x ≥ 0 (resp. x ≤ 0), if x ≥ 0 (resp. x ≤ 0).
The set of all real intervals of R is denoted IR. Two intervals x1 and x2 are equal if and only if x1 = x2 and x1 = x2. Real arithmetic operations have been extended to intervals [8]:
◦ ∈ {+, −, ∗, /}, x1 ◦ x2 = {x ◦ y | x ∈ x1, y ∈ x2}. An interval vector or box [x] is a vector with interval components. An interval matrix is a matrix with interval components. The set of n−dimensional real interval vectors is denoted by IRn and the set of n × m real interval matrices is denoted by IRn×m. The width w(.) of an interval vector (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) composed of the midpoints of its interval components.
Classical operations for interval vectors (resp. interval matrices) are direct extensions of the same operations for real vectors (resp. real matrices) [8]. 2.2
Given x a box of IRn and a function f from IRn to IRm, an 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: f (x) = {f (ν) | ν ∈ x}, where ν is a real vector of the same dimension as x. Then, the interval function [f ] from IRn to IRm is an inclusion function for f if:
∀x ∈ IRn, f (x) ⊂ [f ](x).
An inclusion function of f can be obtained by replacing each occurrence of a real variable by its corresponding interval and by replacing each standard function by its interval evaluation. Such a function is called the natural inclusion function. A function f generally has several inclusion functions, which depend on the syntax of f . 2.3
Throughout the paper and unless explicitly mentioned, variables are assumed to take values in IRd, where d is the dimension of the variable. Exception is made for overlined and underlined variables that are assumed to take values in Rd, where d is the dimension of the variable. Bold symbols are used to denote multi-dimensional variables (vector or matrices). 3 3.1
The architecture of the preventive maintenance method is shown in Fig. 1. The method relies on two modules: • A diagnosis module that uses the system measured inputs and outputs to compute an estimation of the system’s health status; this is performed by estimating the value of the system parameter vector θ by the means of a behavioral Model Σ of the system. • A prognosis module that predicts the parameter evolution over time by using a Damaging Model Δ and computes the Remaining Useful Life or RUL of the underlying subsystems.
The global model representing the progressive evolution of the system over time, e.g. the Ageing Model, is obtained by putting together the behavioral model Σ and the damaging model Δ. The behavioral model takes the form of a state space model, i.e. a state equation modeling the dynamics of the system state vector, and an observation equation that links the state variables to the observed variables: Σ : x˙ (t) = f (t, x(t), θ(t), u(t)) y(t) = h(t, x(t), θ(t), u(t)) (1)
Ageing model Behavioural model Σ where x(t) is the state vector of the system of dimension nx, u(t) is the input vector of dimension nu, y(t) is the output vector of dimension ny, θ(t) is the parameters vector of dimension nθ. Σ represents the system’s nominal behavior as it is supposed to act when its parameters have not yet suffered any ageing.
The damaging model Δ represents the dynamics of the behavioral model parameters. It is described by the dynamic state equation (2) where the equation states are the system parameters. The equation models how the parameters evolve over time because of the wearing, leakage, etc.: Δ : θ˙ (t) = g(t, θ(t), w, x(t)) (2) where w is a wearing parameter vector of dimension nθ. 3.2
Predicting the evolution of the system’s behavior requires to know a priori how the system will be solicited either by the control system or by external causes (e.g. environmental conditions, temperature, humidity, etc.) This knowledge is generally difficult to obtain. In our approach, we make assumptions about the future solicitations of the system by determining the most usual way the system is intended to be used and we define the notion ofunit cycles. A unit cycle C is defined as a solicitation that repeats in time and that leads to a behavioral sequence that is known to impact system’s ageing.
For example, in the case of a pneumatic valve from the
Space Shuttle cryogenic refueling system, [
One may simultaneously use unit cycles at different time scales, depending on the dynamics of the system and its subsystems. As an example, one may define a “global” unit cycle for a bus as being the journey from the starting station to the terminus, and another unit cycle for the subsystem “bus doors” as being the opening and the closing of the doors at each station. 4
Diagnosis is achieved through SM parameter estimation. This problem assumes that measured outputs ym(ti) generated by the real system on a time horizon ti = t0, . . . , tH of length H × δ, where δ is the sampling period, are corrupted by bounded-error terms that may originate from the system and hence rejected.
If yˆ(ti)j ⊆ ym(ti), ∀ti ∈ {t0, . . . , tH }, then [θ]j is a parameter vector for which the estimation is consistent with the measurements. The box is added to the list of the solution parameter boxes:
P = P ∪ [θ]j .
If none of the two previous conditions is true, i.e.: [yˆ(ti)]j ∩ [ym(ti)] 6= ∅, ∀ti ∈ {t0, . . . , tH }, it means that the parameter box [θ]j is undetermined and that it partially contains solutions. The box is also added to P. Two different labels allow us to keep track of the boxes that are feasible or undetermined. The convex union of these boxes provides the estimation θˆ that encloses the feasible parameter set Θ, i.e. θˆ ⊇ Θ.
The quality of the enclosure depends on the size of the boxes of the partition, in other words on the partition precision. A way to improve the enclosure is to proceed with a partition of the obtained solution θˆ and run another round of consistency tests over the new boxes, and so on recursively. The process of iterating the partition ends when the gain in precision is low with respect to the SM integration and consistency tests computational cost. The method is detailed for a one dimension parameter vector in the next section.
The estimation precision ω(P ) obtained for a given partition P can be evaluated by the following percentage: ω(P ) = mid(θˆ) ./ mid(θˆ) + w(θˆ)/2 parameters varying within specified bounds, bounded noise, or sensor precision. The ym(ti)’s are hence interval vectors of IRny . The SM parameter estimation problem for the system Σ is formulated as finding the setΘ ⊆ Rnθ of real parameter vectors such that the arising outputs y(ti, θ) ∈ Rny hit all the output data sets, i.e.:
θ ∈ Θ ⇔ y(ti, θ) ∈ ym(ti), ∀ti ∈ {t0, . . . , tH }.
Θ is called the feasible parameter set (FPS). SM
parameter estimation problems are generally solved with a
branchand-bound algorithm like SIVIA [
Such algorithms return an overestimation of the FPS
given by the convex union of the candidates that have been
labelled feasible and undetermined [
When considering a SIVIA-based algorithm for dynamical systems like Σ, a critical step is the determination of the inclusion function for the state vector xˆ(ti) at instants ti = t0, . . . , tH , arising from a given candidate parameter vector [θ]∗, from which the [yˆ(ti)], ti = t0, . . . , tH can be computed using the observation equation of Σ. This step relies on set-membership integration for which we have chosen the interval Taylor series integration scheme implemented in the VNODE-LP solver [13]. Although quite well optimized [14], it is well-known that this method is computationally stable only for [θ]∗ of very small size. SIVIA-like parameter estimation algorithms are hence particularly inefficient as they enumerate candidate parameter subspaces starting with the full parameter space. This is why we propose the FRP schema presented in the next section. 4.1
The principle of the FRP method is based on partitioning the parameter search space S(θ). Each part of the partition represents a candidate parameter vector [θ]j for which SM integration of the state equation of Σ provides a conservative numerical enclosure xˆ(ti)j , ti = t0, . . . , tH . The output vector can now be estimated as: yˆ(ti)j = h (t, xˆj , [θ]j , um) , ∀ti ∈ {t0, . . . , tH } (3) We then keep track of the parameters vectors for which the yˆ(ti)j ’s are consistent with the measurements, for all ti = t0, . . . , tH , i.e. the unfeasible ones are discarded. Computing the convex hull then provides us with a minimal and maximal value for the admissible parameter vectors.
The consistency test is defined as testing the intersection of the estimated output vector with the measurements: If ∃ ti ∈ {t0, . . . , tH } s.t. yˆ(ti)j ∩ ym(ti) = ∅, (4) then there is no consistency between the estimation and the measured input um(t) and output ym(t) with the tested parameter vector [θ]j . The parameters box [θ]j is unfeasible (5) (6) (7) (8) (9) where ./ denotes the division of two vectors term by term.
Given two partitions Pi and Pj , one can evaluate the precision gain as:
G(Pj /Pi) = ω(Pj )./ω(Pi). 4.2
The domain value of the parameter vector θ is given by Ω(θ) = [.inf(θbol, θeol), .sup(θbol, θeol)], where .inf and .sup denote the operators inf and sup applied term by term and θbol and θeol are the real vectors whose components are given by θk,bol and θk,eol for each parameter θk, k = 1, . . . , nθ: • θk,bol (bol: “beginning-of-life”) is the factory setting defined by the design specification, • θk,eol (eol: “end-of-life”) is the maximal/minimal admissible value, i.e. the value above/below which the component is considered to have failed and the function is no longer guaranteed.
During the system’s life, the impact of ageing results in the parameter vector value evolving in Ω(θ). Depending on the impact of ageing, its value may decrease or increase with time: θ(ti) = αθ(tj ), ti ≥ tj (10) where α is an nθ dimensional real vector whose components αk are greater or lower than 1 depending on the impact of ageing on the change direction of the parameter.
We assume a health management strategy, which means that diagnosis (and prognosis) is performed according to a given inspection planning at some chronologically ordered P2 P2
P3 times of the system’s life T0, . . . , TF . The parameter vector search space depends on the time and is hence denoted by S(θ)Ti = [S(θ)Ti , S(θ)Ti ], where Ti ∈ {T0, . . . , TF }. Initially, for the first inspection time T0, we set S(θ)T0 = Ω(θ). Diagnosis then returns the estimated parameter value θˆ(T0). For the next inspection time, S(θ) is updated by taking the parameter value estimation into account as follows: • if αk > 1, θk,bol is replaced by inf(θˆk(T0), θˆk(T0)), • if αk < 1, θk,bol is replaced by sup(θˆk(T0), θˆk(T0)), and one of the bounds of the components of S(θ)T1 remains equal to θk,eol.
In the general case, when considering the inspection time Ti, S(θ)Ti is hence obtained with θˆ(Ti−1) as follows: • if αk > 1, inf(θˆk(Ti−1), θˆk(Ti−1)) is replaced by • if αk < 1, sup(θˆk(Ti−1), θˆk(Ti−1)) is replaced by inf(θˆ(Ti), θˆ(Ti)), sup(θˆk(Ti), θˆk(Ti)). 4.3
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 the bearing wearing and the clogging of the environment. In the general case, this kind of knowledge must be brought by an expert of the system and/or the manufacturer.
Let us consider the first inspection time and the initial search space S(θ)T0 given by the domain value of the parameter Ω(θ) = [θbol, θeol]. The search space is partitioned into boxes, in our case intervals (cf. Fig. 2).
The dynamic equation of Σ is integrated on the time window ti = t0, . . . , tH , where tH = T0, as many times as the number of intervals in the partition P1. The number of intervals is defined by the partition factor (P1), which equals 1/15 in our example (cf. Fig. 2). We start with [θ]1 = [θbol, θbol + pw], where pw = (P1)w(S(θ)T0 ) is the width of the partition intervals, then proceed with the subsequent intervals [θ]j . For each interval, we get an estimation of the state vector at times ti = t0, . . . , tH , denoted as xˆ(t0 . . . tN )j , and obtain yˆ(t0 . . . tN )j thanks to the observation equation (1). This latter is tested for consistency against the measurements ym(t0 . . . tN ).
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 yellow-colored parts, respectively, in Fig. 2).
P1
The convex union of feasible and undetermined intervals provides a guaranteed estimation θˆ = hθˆ, θˆi of the admissible values for θ. We iterate the process by creating a new partition P2 of hθˆ, θˆi with a precision (P2) = 1/10 (cf. Fig. 3).
P1 P1 P
We proceed as above for each interval of P2 in order to refine the bounds ofhθˆ, θˆi and find a more precise enclosure of feasible parameter solution (see Fig. 3).
We iterate the process until the precision gain G(Pi+1/Pi) is greater then a given threshold, as it is shown in Fig. 4.
The method can be easily generalized to a system whose parameter vector has dimension nθ > 1. The computing cost is proportional to the number of boxes that are tested, i.e. Pin=P1 1/ (Pi), where nP is the number of partitions.
Let’s notice that the partition may be non-regular. For example, for a slowly ageing parameter, one may choose small boxes for the values of θ that are close to θbol and larger ones for the values close to θeol. The result is guaranteed even if the partition has not been properly chosen or if the parameter has evolved in a non expected way, although the computation cost may be higher.
The convex union provides a poor result if the set of admissible values is made of several mutually disjoint connected sets, as shown in Fig. 5. The algorithm may test some boxes that have already been rejected by the tests of the previous partition. This drawback could be addressed by defining the solution as a list of boxes whose labels (unfeasible, feasible, or undetermined) are inherited by the next partition boxes.
Solution The global model (Σ + Δ) assumes that the parameters of the behavior model Σ given by (1) evolve in time, and that their evolution is represented by the degradation model Δ given by the dynamic equation (2) that is recalled below: Δ : θ˙ (t) = g(t, θ(t), w, x(t)). Δ provides the dynamics of the parameter vector as a function of the state of the system x(t) and of a degradation parameter vector w that allows one to tune the degradation for each of the considered parameters.
The global model (Σ + Δ), in the form of a dynamic model with varying parameters, cannot be directly integrated by VNODE-LP. An original method, coupling the two models Σ and Δ iteratively is proposed in the following. The method is illustrated by Fig. 6 and used to determine the degradation suffered by each parameter during one unit cycle as defined in Section 3.2.
Let us denote uC (t), t ∈ [τ, τ + dC ], the system input stress during one unit cycle C . As shown in Fig. 6, the following steps are iteratively executed, every iteration corresponding to a computation step given by the sampling period δ: 1. The normal behavior model Σ is used first with input u(t) = uC (τ ) to compute the state x(τ ) and the output y(τ ); 2. The parameters are updated with the degradation model Δ using the value of the state determined previously, i.e. θ(τ ) is computed; 3. The parameters of the behavior model Σ are updated with θ(τ ); 4. The next stress input value uC (τ + δ) is considered, and so on until the end of the cycle, i.e. until the last value of the cycle uC (τ + dC ) is reached.
Σ Δ where nθ is the number of parameters of the system. Let’s assume the cycle i, then D maps θi into D(θi) = θi+1, which is the value of θ after one unit cycle.
D is nonlinear. Thus the value of the parameter vector after one cycle θi+1 depends on the initial value θi. Indeed, we know that a system generally degrades in a nonlinear fashion. We must hence compute θi+1 for all possible values of the parameter vector θi.
For this purpose, the domain value Ω(θk) of each parameter θk is partitioned into Nk intervals. Nk is chosen sufficiently large to reduce non conservatism of the interval function D. The domain value of the parameter vector θ is hence partitioned into NΠ = Πnθ
k=1Nk possible boxes that must be fed as input to D. Let us for instance consider a two parameters vector and its beginning-of-life and end-of-life values as follows: θ = θ1 , θbol = θ2 1 1 , θeol = 94 and N = 32 , (12) then, if we select the partition landmarks as {5} for θ1 and {2, 3} for θ2 2, D must be run for the following 6 box values:
For each of these box values taken as input for cycle i, i.e. θi = [θ]l, l = 1, . . . , 6, D returns the (box) value θi+1 after one unit cycle. This computation is then projected on each dimension to obtain a set of nθ tables, Dθk , k = 1, . . . , nθ, that provide the degradation of each individual parameter θk after one unit cycle. The RUL, understood as a RUL for the whole system, can now be determined by computing the number of cycles that are necessary for the parameters to reach the threshold defining the end-of-life (cf. Fig. 7).
Diagnosis at inspection time Tk
Yes RUL = i
No
For the cycle i = 0, θ0 is initialized with θˆ, which is the result of the parameter estimation computed by the diagnosis 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 may take several forms as explained in Section 5.3. If the test is true, then the index i is the number of cycles required to reach the degradation threshold, so RUL = i.
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 = 1, . . . , NΠ, of the partition. We propose to compute θi+1 by assuming that the mapping between θi and θi+1 is linear in every domain l of the partition. Considering p ∈ Rnθ , D(p) is approximated as follows: ∀θ ∈ [θ]l, D(θ) ≈ a θ + b, l=1, . . . , NΠ (14) where a= w(D([θ]l))./w([θ]l), b= D([θ]l) − a[θ]l, and is the product of two vectors term by term.
Equation (14) is applied to θi and θi to obtain an approximation of D(θi).
2Notice that the intervals issued from the partitioning are not required to be of equal length.
The set-membership test implemented with the order relation may take several forms. For instance, if the test θi θeol is interpreted as: ∃k ∈ {1, . . . , nθ} | i θk ≥ θk,eol if αk > 1 or θik ≤ θk,eol if αk < 1, (15) then it means that the bound of the interval value of at least one parameter θk is above or below its end-of-life threshold value θk,eol. The RUL is then qualified as the “worst case RUL”, which means that the RUL indicates the earliest cycle at which the system may fail.
One can also test whether the value higher bound of one of the parameters is higher than its end-of-life threshold, that is to say: ∃k ∈ {1, . . . , nθ} |
θik ≥ θk,eol if αk > 1 or θik ≤ θk,eol if αk < 1. (16) The RUL then represents the cycle at which it is certain that the system will fail.
It is obviously possible to combine these different tests applied to the different individual parameters depending on their criticality. 6
The case study is a shock absorber that consists of a moving mass connected to a fixed point via a spring and a damper as illustrated by Fig. 8. The movement of the mass takes place in the horizontal plane in order to eliminate the forces due to gravity. Aerodynamic friction forces are neglected. k c x m where m is the mass, ~a is the acceleration, F~k is the spring biasing force, F~c is the friction force exerted by the damper and ~u is the force applied on the mass. Expressing the forces and the acceleration as a function of the position of the mass x(t), we get: c
k x¨(t) + x˙ (t) +
x(t) = u(t) m m where k is the spring stiffness constant (N/m), m is the mobile mass (kg), and c is the damping coefficient (Ns/m). (18) is a second order ODE. Let us rewrite c k (18) m = 2ζω0 and m = ω02 and we get r k
c ω0 = m and ζ = 2√km .
The impulse response of such system depends on the value of ζ: • if ζ = 0, then the answer is a sinusoid; • if 0 < ζ < 1, then the answer is a damped sinusoid; • if ζ ≥ 1, then the answer is a decreasing exponential. The state model is given by the equation: X˙ (t) = 0 −k m 1 −c X(t) + m 0 1 U (t) 1 0 Y (t) = 0 1 X(t) with X(t) = [x(t), x˙ (t)]T , and the transfer function is: X(p) U (p) =
1 p2 + mc p + mk .
An example of bounded error step response obtained with VNODE-LP with a sampling parameter δ = 0.1 s, c = 1, m = 2 and k = [3, 9 ; 4, 1] is shown in Fig. 9a. There, ζ ' 0.177 and the step response is a damped sinusoid. Because k is assumed to have an uncertain value bounded by an interval, the outputs are in the form of envelops. 6.2
In the case study, a unit cycle is defined by the application 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 lasts 20s and cancels at t0+ 25s as shown by the red curve of Fig. 9b. The cycle ends at t0 + 50s.
Fig. 9b presents the system’s response for a spring constant k = [3.9, 4.1] N/m, a mass m = 2 kg, a damping coefficient c = 10 Ns/m, and initial speed and position equal to zero. The response is a decreasing exponential. 6.3
The degradation model chosen is the ageing of the damper cylinder. It is represented by a reduction of the damping coefficient proportional to the velocity of the mass[15]: c˙ = βx˙ , β < 0.
The more the spring is used, the weaker it becomes, characterized by the change in the damping coefficient. 6.4
The FRP parameter estimation method presented in Section 4 has been used with the measures shown in Fig. 9c. These measures were obtained for θ = " c # "5# k = 4 m 2
The goal is to estimate the damping coefficientc and the stiffness constant k. The search space is defined by the interval [4 9] for c and [3.5, 9] for k. The value of m is assumed to be known m = 2. Using the notation introduced above, we have: θbol = " 4 # 3.5 , θeol = 2 (a) Step response for ζ = 0.177. (b) Unit cycle for the case study. (c) Measured input and output. left and [θ]j = [5, 5, 1], [4, 4, 1], 2 T on the right. On the left figure, one can see that there is no intersection between the estimate and the measurement for the position, hence the box used for the simulation is rejected. On the right, there is an intersection between the measurement and the estimation for all time points, but the estimate is not included in the measure envelop, hence the parameter box is considered undetermined.
Figure 10 – Estimation results with a rejected parameter box (left) and an indetermined box (right)
The results for partition P1 are presented in Fig. 11a and
we obtain a first estimation forθ:
θˆ =
" [4, 1, 5.8] #
[
2 The estimation precision for partition P1 is given by: ω(P1) = mid(θˆ) ./( mid(θˆ) + w(θˆ)/2) =
The precision is now ω(P2) = [0.9, 0.96, 1]T , and the precision gain is G(P2/P1) = [0.056, 0.025, 0]T . The values for the gain indicate that partitioning a third time might be quite inefficient. To confirm this fact, let us perform a third partition P3, whose precision is increased by a factor of 5, i.e. = 0.02 (cf. Fig. 11c). The new estimation for θ is T θˆ = [4.548, 5.526], [3.872, 4.132], 2 , and the precision gain is G(P3/P2) = [0.0073, 0.0036, 0]T . As expected, the gain is quite negligible with respect to the computation time increase.
In this section we apply the set-membership method described in Section 5.2 to compute the RUL for the damping coefficientc.
The damper is assumed to fail when c ≤ ceol = 2. The degradation model (21) with β = −0, 13 allows us to determine the degradation table Dc for the parameter c for a unit cycle:
ci D(ci) = ci+1 ci D(ci) = ci+1
[
The number of elements of the partition has been chosen relatively small to better illustrate the method. In a real situation, this number should be high in order to obtain less conservative predictions.
The value of c has been previously estimated and is
The graph of Fig. 12 allows us to approximate the predicted value after one unit cycle:
D(cˆ) = c1 = [4.4787, 5.4481].
The next iteration of the algorithm allows us to compute c2, etc. After 30 iterations, we obtain c30 = [1.7665, 3.4235].
3The coefficient β has been chosen arbitrarily to illustrate the approach; it does not represent the real ageing of a damper. (a) Partition P1 (b) Partition P2 (c) Partition P3
Since c30 < ceol, we get RU L = 30 cycles. After the 44th iteration, we get c44 = [0.037591, 1.985928]. We then have c44 < ceol and hence RU L = 44 cycles. The RUL of the damper is hence given by:
RU L = [30, 44] cycles. 7
This paper addresses the condition-based monitoring and prognostic problems with a new focus that trades the traditional statistical approach by an error-bounded approach. It proposes a two stages method whose principle is to first determine the health status of the system and then use this result to compute the RUL of the system. This study uses advanced interval analysis tools to obtain guaranteed results in the form of interval bounds for the RUL.
The results for the case study demonstrate the feasibility of the approach. The next step is to adapt the FRP-based SM parameter estimation algorithm in order to output a list of boxes instead of a single box given by the convex hull of the boxes. The convex hull is indeed a very conservative approximation when the solution set is not convex.
The second stream of work is to consider contextual conditions and their associated uncertainties. Environmental conditions, like weather, different usage, etc. may indeed significantly affect the stress input and prognostics results. [13] N. Nedialkov. VNODE-LP, a validated solver for initial value problems for ordinary differential equations. [15] Ian M Hutchings. Tribology: friction and wear of engineering materials. Butterworth-Heinemann Ltd, 1992.