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.
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 [9].
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 by w(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 x 1 and x 2 are equal if and only if x 1 = x 2 and x 1 = x 2 . Real arithmetic operations have been extended to intervals [8]:
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 IR n and the set of n × m real interval matrices is denoted by IR n×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].
Given x a box of IR n and a function f from IR n to IR m , 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:
where ν is a real vector of the same dimension as x. Then, the interval function
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 .
Throughout the paper and unless explicitly mentioned, variables are assumed to take values in IR d , where d is the dimension of the variable. Exception is made for overlined and underlined variables that are assumed to take values in R d , where d is the dimension of the variable. Bold symbols are used to denote multi-dimensional variables (vector or matrices).
3 Principle of the Set-Membership Health Management Method
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: where x(t) is the state vector of the system of dimension n x , u(t) is the input vector of dimension n u , y(t) is the output vector of dimension n y , θ(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.:
where w is a wearing parameter vector of dimension n θ .
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 of unit 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, [2] defines a unit cycle as the opening of the valve, the filling of the tank and the valve closing when the tank is full. In the case of an aircraft, an unit cycle may be chosen to be a flight: it starts with the plane take-off, a cruising stage and landing.
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.
Diagnosis is achieved through SM parameter estimation. This problem assumes that measured outputs y m (t i ) generated by the real system on a time horizon t i = t 0 , . . . , t H of length H × δ, where δ is the sampling period, are corrupted by bounded-error terms that may originate from the system parameters varying within specified bounds, bounded noise, or sensor precision. The y m (t i )'s are hence interval vectors of IR ny . The SM parameter estimation problem for the system Σ is formulated as finding the set Θ ⊆ R n θ of real parameter vectors such that the arising outputs y(t i , θ) ∈ R ny hit all the output data sets, i.e.:
Θ is called the feasible parameter set (FPS). SM parameter estimation problems are generally solved with a branchand-bound algorithm like SIVIA [10] that enumerates candidate box solutions thanks to a rooted tree and assumes the full parameter space as the root. At every node, the set of ŷ(t i ), t i = t 0 , . . . , t H , arising from the considered box parameter vector [θ] * , i.e. solution of Σ for any real θ ∈ [θ] * , is checked for consistency against the measurements and labelled feasible, unfeasible or undetermined. Unfeasible candidates are rejected while undetermined candidates are split and checked in turn until the set precision of the candidate solutions is below a given threshold ε provided by the user.
Such algorithms return an overestimation of the FPS given by the convex union of the candidates that have been labelled feasible and undetermined [11]. Interestingly, the convex union may consist of one set or more, which means that the systems does not need to be identifiable in the classical sense [12].
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(t i ) at instants t i = t 0 , . . . , t H , arising from a given candidate parameter vector [θ] * , from which the [ŷ(t i )], t i = t 0 , . . . , t H 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.
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(t i ) j , t i = t 0 , . . . , t H . The output vector can now be estimated as:
We then keep track of the parameters vectors for which the ŷ(t i ) j 's are consistent with the measurements, for all t i = t 0 , . . . , t H , 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:
then there is no consistency between the estimation and the measured input u m (t) and output y m (t) with the tested parameter vector [θ] j . The parameters box [θ] j is unfeasible and hence rejected.
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:
If none of the two previous conditions is true, i.e.:
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:
where ./ denotes the division of two vectors term by term. Given two partitions P i and P j , one can evaluate the precision gain as: G(P j /P i ) = ω(P j )./ω(P i ).
(9)
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:
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 times of the system's life T 0 , . . . , T F . The parameter vector search space depends on the time and is hence denoted by S(θ) Ti = [S(θ) Ti , S(θ) Ti ], where T i ∈ {T 0 , . . . , T F }. Initially, for the first inspection time T 0 , we set S(θ) T0 = Ω(θ). Diagnosis then returns the estimated parameter value θ(T 0 ). 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 sup( θk (T 0 ), θk (T 0 )), 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 T i , S(θ) Ti is hence obtained with θ(T i−1 ) as follows:
) is replaced by sup( θk (T i ), θk (T i )).
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 t i = t 0 , . . . , t H , where t H = T 0 , as many times as the number of intervals in the partition P 1 . The number of intervals is defined by the partition factor (P 1 ), which equals 1/15 in our example (cf. Fig. 2). We start with [θ] 1 = [θ bol , θ bol + pw], where pw = (P 1 )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 t i = t 0 , . . . , t H , denoted as x(t 0 . . . t N ) j , and obtain ŷ(t 0 . . . t N ) j thanks to the observation equation (1). This latter is tested for consistency against the measurements y m (t 0 . . . t N ).
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). The convex union of feasible and undetermined intervals provides a guaranteed estimation θ = θ, θ of the admissible values for θ. We iterate the process by creating a new partition P 2 of θ, θ with a precision (P 2 ) = 1/10 (cf. Fig. 3). We proceed as above for each interval of P 2 in order to refine the bounds of θ, θ and find a more precise enclosure of feasible parameter solution (see Fig. 3). We iterate the process until the precision gain G(P i+1 /P i ) 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. n P i=1 1/ (P i ), where n P 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.
The prognosis phase consists is calculating the number of cycles remaining before anomaly, which is also called the Remaining Useful Life or RUL. To optimally adapt this calculation to the system's life requires the knowledge of the health status of the system at the current time, which was the topic of Section 4.
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).
Yes
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 ∈ R n θ , D(p) is approximated as follows:
where a=
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 ).
The set-membership test implemented with the order relation may take several forms. For instance, if the test θ i θ eol is interpreted as:
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:
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 Case study
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.
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:
where k is the spring stiffness constant (N/m), m is the mobile mass (kg), and c is the damping coefficient (Ns/m). ( 18 The impulse response of such system depends on the value of ζ:
, then the answer is a sinusoid;
, then the answer is a damped sinusoid;
, then the answer is a decreasing exponential. The state model is given by the equation:
with X(t) = [x(t), ẋ(t)] T , and the transfer function is:
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.
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 t 0 +5s, where t 0 is the cycle starting time. The force lasts 20s and cancels at t 0 + 25s as shown by the red curve of Fig. 9b. The cycle ends at t 0 + 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.
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]:
The more the spring is used, the weaker it becomes, characterized by the change in the damping coefficient.
The FRP parameter estimation method presented in Section 4 has been used with the measures shown in Fig. 9c. These measures were obtained for
The goal is to estimate the damping coefficient c 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:
The partition P 1 is achieved with a precision (P 1 ) = 1/10 for the two parameters to be estimated c and k. 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. The results for partition P 1 are presented in Fig. 11a and we obtain a first estimation for θ: .
The estimation precision for partition P 1 is given by: ω(P 1 ) = mid( θ) ./( mid( θ) + w( θ)/2) = 0.85 0.94 1
The first estimation for θ is used as the search space for partition P 2 , whose precision is increased by a factor of 10, i.e. (P 2 ) = 0, 1. The obtained estimation results are shown in Fig. 11b.
The estimation is refined as: The precision is now ω(P 2 ) = [0.9, 0.96, 1] T , and the precision gain is G(P 2 /P 1 ) = [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 P 3 , whose precision is increased by a factor of 5, i.e. = 0.02 (cf. Fig. 11c). The new estimation for θ is θ = [4.548, 5.526], [3.872, 4.132], 2
T , and the precision gain is G(P 3 /P 2 ) = [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 coefficient c. [0, 0.9755] (25) After proceeding to the linear interpolation given by ( 14), the graphical representation of c i+1 as a function of c i is given by Fig. 12. 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
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.
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:
∆ 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 u C (t), t ∈ [τ, τ + d C ], 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) = u C (τ ) to compute the state x(τ ) and the output y(τ ); The above algorithm defines the function:
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 N k intervals. N k 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=1 N k 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:
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: