Pascal Paulmier

Jean-François Maire

Myriam Kaminski

Jean-Michel Roche

Keywords: Pronostic à base de modèles, filtrage, propagation d'incertitudes, filtres stochastiques, méthodes ensemblistes, analyse de fiabilité, propagation de fissure Model-based prognosis, filtering, uncertainty propagation, stochastic filter, set-membership framework, reliability analysis, fatigue crack growth . . . Model-based prognosis, filtering, uncertainty propagation, stochastic filter, set-membership framework, reliability analysis, fatigue crack growth

Les travaux présentés dans ce mémoire concernent le développement de méthodes de pronostic à base de modèles. Le pronostic à base de modèles a pour but d'estimer le temps qu'il reste avant qu'un système ne soit défaillant, à partir d'un modèle physique de la dégradation du système. Ce temps de vie restant est appelé durée de résiduelle (RUL) du système.

Le pronostic à base de modèle est composé de deux étapes principales : (i) estimation de l'état actuel de la dégradation et (ii) prédiction de l'état futur de la dégradation. La première étape, qui est une étape de filtrage, est réalisée à partir du modèle et des mesures disponibles. La seconde étape consiste à faire de la propagation d'incertitudes. Le principal enjeu du pronostic concerne la prise en compte des différentes sources d'incertitude pour obtenir une mesure de l'incertitude associée à la RUL prédite. Les principales sources d'incertitude sont les incertitudes de modèle, les incertitudes de mesures et les incertitudes liées aux futures conditions d'opération du système. Afin de gérer ces incertitudes et les intégrer au pronostic, des méthodes probabilistes ainsi que des méthodes ensemblistes ont été développées dans cette thèse.

Dans un premier temps, un filtre de Kalman étendu ainsi qu'un filtre particulaire sont appliqués au pronostic de propagation de fissure, en utilisant la loi de Paris et des données synthétiques. Puis, une méthode combinant un filtre particulaire et un algorithme de détection (algorithme des sommes cumulatives) a été développée puis appliquée au pronostic de propagation de fissure dans un matériau composite soumis à un chargement variable. Cette fois, en plus des incertitudes de modèle et de mesures, les incertitudes liées aux futures conditions d'opération du système ont aussi été considérées. De plus, des données réelles ont été utilisées. Ensuite, deux méthodes de pronostic sont développées dans un cadre ensembliste où les erreurs sont considérées comme étant bornées. Elles utilisent notamment des méthodes d'inversion ensembliste et un observateur par intervalles pour des systèmes linéaires à temps discret. Enfin, l'application d'une méthode issue du domaine de l'analyse de fiabilité des systèmes au pronostic à base de modèles est présentée. Il s'agit de la méthode Inverse First-Order Reliability Method (Inverse FORM).

Pour chaque méthode développée, des métriques d'évaluation de performance sont calculées dans le but de comparer leur efficacité. Il s'agit de l'exactitude, la précision et l'opportunité.

À mon père, ma mère, mes frères, et ma moitié.

«Le hasard gouverne un peu plus de la moitié de

Introduction Context

With the increasing competition in aircraft industry, manufacturers engaged in a fierce commercial battle which drove them to find an optimum balance between aircraft performance and cost-effectiveness. However, users safety and mission success remain of prime importance and should be taken into account as well. To this end, the implementation of optimized maintenance strategies under cost, availability and safety constraints has become a priority for many industrials.

There exist two major types of maintenance strategies: reactive maintenance and proactive maintenance. Reactive maintenance is performed after the occurrence of failures, while proactive maintenance is performed before the onset of failures.

Reactive maintenance, also known as corrective maintenance, consists in repairing or replacing the failed equipment. Therefore, costs and investment required to implement a maintenance program are minimized. However, unexpected failures can lead to a significant loss of revenue because of downtime, equipment, and labor costs. Moreover, it is an inefficient use of time as maintenance staff is forced to work overtime or to stop their current activities in order to fix the problem. For all these reasons, this kind of maintenance strategy is restricted to inexpensive components that are easy to replace and whose failure does not involve safety issues.

Proactive maintenance is particularly used on the critical components of the system and is performed while the latter is still working. The risk of failures are anticipated and eliminated to avoid unexpected system breakdown. With such a strategy, the cost and safety issues encountered by applying a reactive maintenance strategy could be minimized [START_REF] Mobley | An introduction to predictive maintenance[END_REF]]. Proactive maintenance policies can be classified into three categories: preventive maintenance, condition-based maintenance and predictive maintenance. Each of these policies is governed by a decision variable that is used to trigger the maintenance procedure.

Preventive maintenance is regularly performed following a predefined schedule, time is therefore the decision variable. Periodic inspections, repairs or replacement are realized on equipments that were previously selected, even if no fault has been detected. The planning is established from equipment history, reliability analysis and mean time between failures (MTBF). It may be noted that even if time is the most frequently used unit of measurement, depending on the situation this may be the number of cycles, the distance traveled, the number of flights, etc. However, this maintenance approach is not optimal because in some cases the equipments ar replaced either too early or too late, therefore increasing the costs.

Condition-based maintenance relies on condition-monitoring which consists in monitoring indicators that reflect the current health state of in-service equipment. It allows to determine if the performance of one or more components is deteriorating, which is a sign of future failure. The system is maintained in-service until one or more indicators reach a predefined threshold, which triggers maintenance actions. Each of these thresholds correspond to a decision variable. The indicators depend on the system and can be observed through different processes: vibration analysis, oil analysis (measure the number and size of particles), acoustic measurements (detection of gas or liquid leaks), electrical measurements (motor current), etc. Therefore, condition-based maintenance enables one to repair faulty equipment and not prematurely as in preventive maintenance. However, the moment at which a failure occurs is uncertain, if not unknown. It would be more efficient to know this time in order to keep equipments working as long as possible.

Predictive maintenance also involves condition-monitoring and solves the aforementioned problem related to condition-based maintenance. Indeed, it takes a step further by forecasting the evolution of the health indicators and estimating when they will cross a critical threshold causing the system failure. The time available until the monitored component no longer behaves properly is its remaining useful life (RUL). Maintenance is triggered depending of the estimated RUL, which is thus the decision variable. In this way, maintenance staff could schedule when exactly an operation shall be initiated, which means neither too early nor too late. This allows to reduce unnecessary preventive actions and to increase system reliability. Therefore, predictive maintenance is the perfect candidate for an optimized maintenance strategy under cost, availability and safety constraints.

Predicting the evolution of the system degradation and estimating its RUL is the core of predictive maintenance, and is referred to as prognosis [START_REF] Vachtsevanos | Intelligent fault diagnosis and prognosis for engineering systems[END_REF]]. The work presented in this thesis is a contribution to this research field.

Problem statement

Prognosis methods are usually classified into three categories [START_REF] Liu | A review on fault prognostics in integrated health management[END_REF]]: knowledge-based, data-driven and model-based. In the knowledge-based approaches, expert judgment, historical and empirical failure data are used to deduce the RUL. Data-driven approaches extract features from operating data such as current, temperature, or vibration signals, and then employ statistical and machine learning techniques to estimate the RUL. Finally, model-based prognosis approaches compute an estimation of the RUL by using a physics-based mathematical model of the system degradation.

Based on the available knowledge about the system, one can choose among these three approaches or combine them. Moreover, various prognosis approaches can even be applied within the same system, depending on the considered subsystem. Datadriven approaches are the most used ones in the literature [START_REF] Lee | Introduction to data-driven methodologies for prognostics and health management[END_REF]. Indeed, the degradation models are built on the basis of available data and machine learning methods (neural networks, Gaussian processes, support vector regressions, fuzzy logic, etc.). Therefore, whatever the system or subsystem under study, data-driven techniques can always be applied as long as enough sensors and data are accessible. However, when a physics-based model of the system is available, a model-based approach is recommended. Indeed, the ability of model-based methods to adapt the model to the evolution of the system degradation ensures an accurate prognosis when more information about the degradation becomes available. Moreover, it avoids the recording and storage of a massive amount of data. In this context, this thesis is focused on model-based prognosis approaches.

One major challenge in model-based prognosis is the measure of the uncertainty affecting RUL prediction [START_REF] Sankararaman | Significance , interpretation , and quantification of uncertainty in prognostics and remaining useful life prediction[END_REF]]. Indeed, because of cost and safety issues, the predicted RUL must be as reliable as possible to allow risk-based decisions. For this purpose, prognosis methods should account for the different uncertainty sources that affect the estimation of the current and future health state of the system. The main types of uncertainty that inevitably influence RUL prediction are measurement uncertainty, modeling uncertainty and future loading uncertainty. Measurement uncertainty is due to sensor inaccuracy, modeling uncertainty is characterized by the difference between the degradation model and its real behavior, and future loading uncertainty is caused by various environmental factors that could possibly affect the evolution of the degradation. In order to deal with these uncertainties, appropriate uncertainty representation and uncertainty propagation techniques should be chosen. Depending on the desired modeling and simulation framework, one can choose among various approaches [START_REF] Ayyub | Uncertainty modeling and analysis in engineering and the sciences[END_REF] which are for example based on interval analysis, fuzzy set theory, probability theory, etc.

In the literature, modeling and measurement uncertainty are mainly handled with probabilistic methods. The random variables are represented using probability distributions. Particle filtering-based techniques are commonly used as they can deal with nonlinearities. However, Monte-Carlo simulation-based approaches require huge computational efforts and are not therefore practicable for online RUL prediction.

Introduction

Objective and outline of the thesis

Based on the aforementioned insights, the objective of this thesis is to adapt different uncertainty propagation methods to model-based prognosis. These methods are suitable for online RUL prediction and the problem of future loading uncertainty is handled without any computational expensive technique. To demonstrate the efficiency of the proposed methodologies, various degradation mechanisms from linear to nonlinear ones are considered.

In order to highlight the contributions of this work, the following outline was adopted:

• Chapter 1: A detailed presentation of prognosis principles is given. The objective of this chapter is to provide a comprehensive overview of the main issues and challenges to consider when developing model-based prognosis methods.

• Chapter 2: Two classical filtering-based methods, namely an extended Kalman filter (EKF) and a particle filter (PF), are investigated in the framework of nonlinear fatigue crack growth prognosis based on the Paris' law and synthetic data. This example is widely used in the literature as it is an interesting test case for aerospace industry. Although the chosen case study and methods are classical, the aim of this chapter is to:

1. clearly establish each step of the prognosis process, 2. understand how to implement and to interpret the proposed prognosis performance evaluation metrics (accuracy, precision and timeliness),

3. identify the challenges involved in RUL prediction, both when estimating the current degradation state and during its future state prediction.

Generally speaking, this chapter has allowed the establishment of an improved model-based prognosis method based on particle filtering that takes loading uncertainty into account.

• Chapter 3: Based on the knowledge gained in Chapter 1, a more realistic case study that takes loading uncertainty in addition to modeling and measurement uncertainty is presented. Two scenarios where the system is subject to both constant and variable unknown and unmeasured input loading is considered.

The proposed methodology allows to jointly estimate the degradation state and the loading parameters. For this purpose, a particle filter and a cumulative sum (CUSUM) algorithm are combined to deal with potential abrupt changes in the estimated loading parameters values. To illustrate the efficiency of the method, fatigue crack growth prognosis in fiber-reinforced composite materials is performed and real data from fatigue tests are used. The contributions in this chapter are interesting not only for the prognosis research field but also for the aerospace industry. Indeed, composite materials are widely used for aircraft structures, and the estimation and prediction of unmeasured loading conditions is an essential task that considerably improves the applicability of a prognosis algorithm. Moreover, neither computationally expensive technique such as finite element method nor surrogate models that leads to additional uncertainty were introduced. Although the results obtained with the particle filtering-based method are satisfactory, another kind of uncertainty propagation method is investigated in the next chapter.

• Chapter 4: In this chapter, two set-membership state estimation techniques based on interval analysis are applied to the problem of RUL prognosis. In this context, uncertainty propagation is performed without any stochastic assumption and only the bounds of the noises and disturbances are required. This allows the estimation of the degradation state in a guaranteed way.

In the first approach, a set-membership methodology using constraints satisfaction is adopted to estimate and then forecast the degradation state. For this purpose, the feasible domain of the model parameters is first estimated based on the available measurements and the Waltz filtering algorithm. Then, using the natural inclusion function of the degradation model and the projection of the outer approximation of the feasible domain of the parameters, the future degradation behavior is predicted and the interval that contains the RUL of the system is finally computed. The proposed methodology is illustrated using the same nonlinear crack growth example as in Chapter 2.

The second method belongs to the category of set-membership observers. In this case, the degradation cannot be directly measured. It is considered as a disturbance that affects the system dynamics and is referred to as an unknown input. Thus, an unknown input interval observer for linear time-invariant discrete systems is proposed for the first time to tackle joint state and unknown input estimation. Once the bounds of the degradation state are reconstructed, they are propagated through time to obtain the interval containing the RUL of the system. The efficiency of the prognosis method is demonstrated on a suspension system subject to crack propagation that affects its stiffness constant.

• Chapter 5: In this chapter, another kind of uncertainty propagation method based on probability theory is presented. Indeed, Monte Carlo simulation-based approaches require huge computational efforts and are time consuming. An alternative to sampling-based methods, namely the inverse first-order reliability method (Inverse FORM), is then proposed. This is an analytical method originally developed in the structural reliability analysis field, which aims at evaluating the failure probability of a structure. Here, a method to compute the RUL probability density function (pdf) based on the Inverse FORM algorithm is proposed. This methodology can be applied if an expression of the RUL Introduction derived from the explicit dynamical degradation model is available. In a first step, a filter is employed for the estimation of the current degradation state with the available measurements. Then, the Inverse FORM algorithm is used to find for which parameter values the system fails, with reference to a specified failure probability level. Finally, these parameter values and the explicit expression of the RUL allow to calculate the desired percentiles of the RUL pdf in terms of the cumulative distribution function according to the user's requirements. The proposed methodology is less time consuming as the computation of the entire pdf of the RUL is not required, and there is no need to propagate the model equation step by step until the threshold is reached. For the sake of illustration, the crack growth model presented in Chapter 2 is used again.

• Chapter 6: A general conclusion which can be drawn from the work presented in this thesis and discusses potential future work arising from this research is proposed in this last chapter.

Chapter 1

Prognosis and uncertainty propagation

Critical systems such as an aircraft or a spacecraft are made of complex components whose malfunction and failure could have unacceptable impacts on the users safety, the mission success and the costs related to maintenance operations. To address these safety and cost issues, a prognosis module should be integrated to these systems in order to continuously assess their state of health and estimate their remaining useful life (RUL). Thus, prognosis is gaining more and more importance in structural and health monitoring (SHM) systems, in condition-based maintenance (CBM) systems and more recently in prognosis and health management (PHM) architectures [START_REF] Tinga | Aligning phm, shm and cbm by understanding the physical system failure behaviour[END_REF].

In this chapter, a review on existing prognosis approaches is provided first. Then, a focus on model-based prognosis methods is made. Finally, performance evaluation metrics that are used all along this thesis to evaluate the proposed prognosis methods are given.

Prognosis approaches

Prognosis aims at estimating the remaining useful life (RUL) of a system, which is the time available until this latter reaches a failure state. Because of the safety and cost benefits involved in adopting a predictive maintenance strategy, prognosis has become an increasingly popular research field. Thus, prognosis methods arising from various research areas have emerged these last few decades. Although all of these methods aim at computing the RUL of a system, they can be distinguished by the type of tools and models that are used.

Many prognosis methods classifications have been proposed, but the one that is both the oldest and the most common identifies three main prognosis approaches Chapter 1. Prognosis and uncertainty propagation [START_REF] Lebold | Open standards for condition-based maintenance and prognostic systems[END_REF][START_REF] Liu | A review on fault prognostics in integrated health management[END_REF]: knowledge-based approaches, data-driven approaches and model-based approaches.

In knowledge-based approaches [START_REF] Biagetti | Automatic diagnostics and prognostics of energy conversion processes via knowledge-based systems[END_REF]], the actual situation is compared to historical failure database or design recommendations of the system considered under similar operating conditions in order to predict the RUL of the system. For this purpose, expert and fuzzy systems are commonly used. Data-driven approaches [START_REF] Si | Remaining useful life estimation -A review on the statistical data driven approaches[END_REF] explore data collected from sensors (current, temperature, vibration signals, etc.) to extract features which will then be used to model the degradation behavior. Statistical and artificial intelligence techniques are used to track, approximate and forecast the evolution of the degradation state. Finally, model-based approaches [START_REF] Luo | Modelbased prognostic techniques [maintenance applications[END_REF] rely on a physics-based model of the degradation to assess the RUL of the system.

Each of these approaches has its advantages and drawbacks. Knowledge-based approaches are easy to implement, but frequent updates are required as new forms of faults that are not yet listed in the database can occur. Moreover, it could be difficult in some cases to convert expert domain knowledge into rules. Data-driven approaches have the ability to transform high-dimensional noisy data into lower-dimensional information for prognosis decisions. However, they are highly-dependent on the quantity and quality of operational data and therefore require a significant storage space. Model-based prognosis approaches can be difficult to set up as a physics-based degradation model is seldom available. However, when it is the case, they can be particularly effective. Indeed, their ability to incorporate physical knowledge of the system (operational conditions such as temperature, etc.) allows continuous observation of their effect on the system behavior.

The choice of the approach to use will depend on the knowledge about the system behavior, its degradation mechanisms, the available data, or also on the knowledge domains of the user. Depending on the available information, it may be interesting to combine different prognosis approaches in order to benefit from the advantages of each. Such strategy has led to the introduction of a new class of prognosis approaches, namely hybrid approaches. As a matter of fact, a new classification of prognosis approaches has been increasingly considered by the PHM community [START_REF] Lee | Prognostics and health management design for rotary machinery systems-reviews, methodology and applications[END_REF][START_REF] Medjaher | Contribution au pronostic de défaillances guidé par des données[END_REF]]: data-driven approaches, model-based approaches and hybrid approaches.

A review of the prognosis methods developed in the context of data-driven and model-based approaches is given in the following. Hybrid approaches are not discussed as they are not yet widespread in the literature, but the interested reader may for example consult [START_REF] Zio | Fatigue crack growth estimation by relevance vector machine[END_REF] or [START_REF] Daroogheh | A hybrid prognosis and health monitoring strategy by integrating particle filters and neural networks for gas turbine engines[END_REF] to have an overview of such approaches.

Data-driven approaches

Data-driven methods use condition monitoring data and historical run-to-failure data in order to track, estimate and predict the degradation state. It is performed by first learning the degradation model with a so-called training dataset. Then, the learned model is used to estimate the current degradation state and predict the RUL of the system. For this purpose, statistical and artificial intelligence techniques can be used. Statistical techniques mainly include Hidden Markov Models (HMM) and regression models based on relevance/support vector machine (RVM/SVM) or Gaussian Process (GP) regression, while for artificial intelligence techniques, Neural Networks (NN) is the most commonly used one for data-driven prognosis [START_REF] An | Practical options for selecting data-driven or physics-based prognostics algorithms with reviews[END_REF].

In the following, some data-driven prognosis methods are presented in order to have an overview of how data-driven prognosis is performed. As data-driven prognosis is not the main topic of this thesis, the interested reader may refer to [START_REF] Si | Remaining useful life estimation -A review on the statistical data driven approaches[END_REF] or [An, Choi, and Kim 2013a] to have a detailed review of the existing data-driven prognosis approaches.

SVM and RVM.

In [START_REF] Caesarendra | Application of relevance vector machine and logistic regression for machine degradation assessment[END_REF], a methodology combining logistic regression (LR) and RVM is developed for the prognosis of bearing degradation, using vibration data. The first step consists in extracting a feature that characterizes the degradation evolution. Kurtosis of the vibration data was chosen. After feature calculation, LR is applied to estimate the failure probability of run-to-failure bearing data. Then, RVM is trained by input data obtained from kurtosis and target vector of failure probability. For this purpose, Gaussian kernel and appropriate kernel-width are used to achieve optimized training process and thus obtain the best model of RVM. Finally, RVM is employed to predict the failure time of the bearing. Simulated and experimental data downloaded from Prognostics Center of Excellence (PCoE) are used to validate the performance of the proposed method, which gives satisfactory results. [START_REF] Patil | A novel multistage support vector machine based approach for li ion battery remaining useful life estimation[END_REF] proposed a twostage approach (classification and regression) for Li-ion battery state pf health (SOH) prognosis. The first stage consists in giving a rough estimate of the RUL using a SVM classification model. Then, based on the results obtained during this first stage, an accurate RUL estimate is computed using a regression model based on support vector regression (SVR). In order to build the classification and regression models, a reduced set of features extracted from voltage and temperature data is used as input variable. Moreover, a multi-layer perceptron based kernel is chosen. The performance of the proposed approach is demonstrated with the PCoE battery cycling data. The results are satisfactory because the the addition of the classification stage enabled the authors to obtain a more accurate RUL estimate and to reduce computational time as the regression is performed on a reduced number of data. Moreover, the algorithm allows simultaneous RUL estimation of multiple batteries. However, the authors have highlighted some drawbacks of the method: the accuracy highly depends on the number of data available for training, and the model is trained off-line. Indeed, developing an on-line learning strategy could improve the model when new data is collected.

In the literature, other data-driven prognosis approaches using SVM or RVM have been applied to various systems. For example, [START_REF] Benkedjouh | Health assessment and life prediction of cutting tools based on support vector regression[END_REF] have used nonlinear feature reduction and SVR to predict the RUL of cutting tools in order to maintain a high quality of the manufactured pieces. In [START_REF] Chang | Prognostics-based led qualification using similarity-based statistical measure with rvm regression model[END_REF], a RVM regression model was used to predict the RUL of light-emitting diodes (LEDs).

Gaussian Process (GP).

In [START_REF] Liu | Data-driven prognostics for lithium-ion battery based on gaussian process regression[END_REF], a methodology for lithium-ion battery prognosis based on GP regression is presented. The authors proposed first to tackle the self-recharge phenomena during the battery lifetime. For this purpose, they selected a covariance function which is a combination of the periodic and diagonal squared exponential covariance functions. Then, they addressed the problem of improving the precision of the prediction when only few data are available. Thus, they adopted a dynamic model approach which consists in updating the GP model when new data are recorded. The methodology was tested using datasets from the PCoE, and its effectiveness in battery health prognosis was demonstrated. The authors somehow suggested to use methods other than the conjugate gradient to estimate the hyper-parameters. [START_REF] Richardson | Gaussian process regression for forecasting battery state of health[END_REF] investigates the ability of GP to perform battery capacity forecasting using datasets from PCoE and data extracted from previous works. First, a classical GP is studied in order to highlight the importance of choosing a suitable kernel function. For this purpose, different combinations of three stationary kernels are considered (squared exponential, Matérn and periodic), and the negative log marginal likelihood of the GP is used to evaluate them. Thus, the additive kernel composed of Matérn 3/2 and Matérn 5/2 was chosen. The simulation results of RUL evaluation showed that GP with the considered additive kernel gives the best results. However, it is suitable for short-term predictions but not for long-term predictions, i.e when there are limited training data available. In order to improve the long-term prediction capacity, the authors proposed in a second study to use a GP with an explicit mean function based on a battery degradation model. The results showed that this second approach could improve the prediction performance of the basic GP. Finally, it was illustrated in a third case study that using a multi-output GP which considers data from multiple cells could improve the forecasting performance. However, the computational cost is increased. Even if the authors studied various advantages of GP for RUL estimation, they proposed other possible improvements such as considering regime changes in a real-time procedure. [START_REF] Liu | Prognostics for state of health estimation of lithium-ion batteries based on combination gaussian process functional regression[END_REF] applied a linear Gaussian process functional regression (LGPFR) model to perform multi-step-ahead SOH prognosis of lithium-ion battery while taking the regeneration phenomena into account. [START_REF] Li | A new prognostics method for state of health estimation of lithiumion batteries based on a mixture of gaussian process models and particle filter[END_REF] proposed a method based on mixture of GP and particle filtering for lithium-ion battery SOH prognosis.

Hidden Markov Models (HMM).

In [START_REF] Zhang | An integrated approach to bearing fault diagnostics and prognostics[END_REF], a methodology for bearing health monitoring is developed. Three main components are combined to achieve the RUL computation: principal component analysis (PCA), HMM and an adaptive prognosis model. PCA is a multivariate statistical approach that analyzes a set of data and catches patterns in order to convert the initial set into a reduced one containing only the most important information, i.e the principal components [START_REF] Smith | A tutorial on principal components analysis[END_REF]]. Thus, PCA has allowed the authors not only to obtain a reduced number of data, but also to generate residual models that can be used for anomaly detection. Then, the features extracted by PCA are used by HMM to generate the state sequence estimate which provides on-line information about the current degradation state. Moreover, the HMM also generates the degradation index which is used by the adaptive prognosis model to perform RUL prediction. For this purpose, a stochastic model of the fault progression that depends on the degradation index was established, then the adaptive algorithm updates the model parameters based on the received degradation index. Real bearing vibration data were used to demonstrate the efficiency of the method, however the authors suggested that more data intensive validation and performance evaluation should be performed.

[ [START_REF] Tobon-Mejia | A data-driven failure prognostics method based on mixture of gaussians hidden markov models[END_REF]] also proposed a data-driven method for bearing RUL prediction. This time, wavelet packet decomposition (WPD) is used for feature extraction and Mixture of Gaussians HMM (MoG-HMM) allows the continuous assessment of the degradation state of the component and the estimation of its RUL. MoG-HMM was used rather than HMM because of its ability to deal with condition monitoring data that are continuous signals. First, during an off-line phase, the WPD technique is used to extract the energy in the different spectrum bands [START_REF] Zarei | Bearing fault detection using wavelet packet transform of induction motor stator current[END_REF]] from vibration measurements. The resulting features are then used to learn several MoG-HMM, each one corresponding to different initial states and operating conditions. Finally, the parameters and the temporal parameters (stay duration in each state) of the MoG-HMM are estimated and the learned MoG-HMM are stored in a database as references of degradation patterns. Then, during the on-line phase, the processed data and the extracted features are continuously fed to the learned models, then model selection is performed to estimate the current degradation state. Using the current state and the corresponding stay duration learned off-line, the RUL and its confidence value are estimated. Data provided by NASA are used to prove the efficacy of the proposed methodology.

[ [START_REF] Dong | A segmental hidden semi-markov model (hsmm)-based diagnostics and prognostics framework and methodology[END_REF], proposed a prognosis method based on Hidden Semi-Markov Model (HSMM) where the stay duration in each state is variable, and is estimated

Chapter 1. Prognosis and uncertainty propagation

during the learning phase. [START_REF] Zaidi | Prognosis of gear failures in dc starter motors using hidden markov models[END_REF] proposed a method based on HMM for gear fault prognosis in DC machines with sparse motor current data.

Neural Networks (NN).

In [START_REF] Wang | Prognosis of machine health condition using neuro-fuzzy systems[END_REF], the ability of two data-driven prognosis methods based on NN to perform prognosis of machine health condition is evaluated. These methods are referred to as recurrent neural networks (RNN) and neuro-fuzzy (NF) systems. First, two popular datasets in the field of time-series prediction are used to compare two methods. One of the datasets is non-linear, non-Gaussian and non-stationary while the other is chaotic and non-periodic. The simulation results highlighted that if the NF system was trained with sufficient data, it can give better forecasting performance than the RNN. Then, the two algorithms were compared through an on-line fault prognosis in a rotary machinery example. The results showed that the NF still gives higher forecasting accuracy than the RNN. [START_REF] Vachtsevanos | Intelligent fault diagnosis and prognosis for engineering systems[END_REF] proposed a data-driven prognosis method that uses dynamic wavelet NN (DWNN). The method is based on two main components that both use a DWNN model as a mapping tool: a virtual sensor that performs fault evaluation (fault dimension, location, etc.) from measurements, and a predictor which forecasts the fault growth using the output of the virtual sensor as inputs. In order to illustrate the feasibility of the proposed method, RUL prognosis of a rolling-element bearing fault was performed. The results were satisfactory, however the authors suggested that more failure data should be used for validation.

More recently, [START_REF] Vatani | Health monitoring and degradation prognostics in gas turbine engines using dynamic neural networks[END_REF] proposed two methods based on NN for degradation prognosis of gas turbine engines. The first one uses a RNN architecture and the second one uses a nonlinear autoregressive neural network (NARNN). The main challenge highlighted by the authors concerns the choice of the parameters required to implement both algorithms (number of input and output delays, percentage of data for training and testing, etc.). The NARNN needs more parameters but less training data than the RNN. The simulation results showed that the NARNN gives better prediction performance because the predicted output of the neural network is closer to the real data with lower error levels.

Many NN-based data-driven prognosis methods have been developed in the literature. In [START_REF] Parthasarathy | Neural network models for usage based remaining life computation[END_REF], model reduction technique based on dynamic NN was applied for RUL estimation of engines while taking the actual operating conditions (temperatures and stresses) into account. [START_REF] Zhang | Fault diagnosis and prognosis using wavelet packet decomposition, fourier transform and artificial neural network[END_REF]] used a combination of WPD, Fourier transform and NN for fault prognosis in critical components of machines.

Prognosis approaches

Summary

To conclude, data-driven approaches might be very effective to perform on-line RUL prognosis. Indeed, no physical model is required, only data related to the studied system are used. However, these techniques highly rely on the availability of a significant amount of run-to failure data set containing several degradation modes to be learned. Thus, when a physical model of the system behavior and its degradation state are available, model-based prognosis methods can be used.

Model-based approaches

When a physics-based degradation model of the system or component is available, model-based prognosis approaches may be used to estimate the current degradation state with measurements and then predict the RUL of the system.

The degradation model contains parameters (e.g. material properties) that might be unknown. Therefore, model-based prognosis methods involve state and parameter estimation. Moreover, real-time estimation and update of the degradation state is required for an accurate RUL prognosis. Thus, taking these requirements into account, most of the methods developed for model-based prognosis are based on stochastic filtering techniques. Indeed, stochastic filters allow joint state-parameter estimation, and they are based on Bayesian inference [Box and Tiao 2011] which uses prior knowledge and measurement data in order to estimate the degradation state pdf. Model-based prognosis methods mainly include Kalman filters [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF]], particle filters [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF], and their derivatives such as unscented Kalman filter [START_REF] Wan | The unscented kalman filter for nonlinear estimation[END_REF]] and so on. The Kalman filter is suitable for linear systems with Gaussian noise terms. If nonlinear systems are concerned, an extended Kalman filter can be used. However, when nonlinear systems with non-Gaussian noise terms are studied, the particle filter is best suited.

Sometimes, the degradation may not be directly measured, but its influence on the system's behavior can be tracked. In this case, models that simultaneously describe the system and the degradation states are adopted. [START_REF] Chelidze | Procedure for tracking damage evolution and predicting remaining useful life with application to an electromechanical experiment system[END_REF] proposed to consider the degradation as a hidden state which evolves slower than the system dynamics. Therefore, the goal of the methods developed in this framework is to track the hidden degradation state in order to estimate its current state and then predict the RUL of the system. They are mainly based on interactive multiple models (IMM) [START_REF] Phelps | Useful lifetime tracking via the imm[END_REF] and unknown input observers (UIO) [START_REF] Darouach | Full-order observers for linear systems with unknown inputs[END_REF].

Finally, another kind of model-based prognosis methods based on reliability analysis drew the attention of some researchers. Indeed, as RUL prognosis is naturally associated to system reliability, the methods developed in this research Chapter 1. Prognosis and uncertainty propagation field were therefore exploited [START_REF] Lemaire | Structural reliability[END_REF]].

In the following, an overview of the above mentioned model-based prognosis methods that can be found in the literature is given. As particle filter-based methods are the most prevalent ones, they will be reviewed in more details.

Stochastic filtering

Particle filter. Because of its ability to deal with nonlinear systems and non-Gaussian noises, the particle filter (PF) is the most used tool in model-based prognosis.

A complete study of the use of PF for on-line failure prognosis was first introduced in the PhD thesis [START_REF] Orchard | A particle filtering-based framework for on-line fault diagnosis and failure prognosis[END_REF]]. The same author proposed in [START_REF] Orchard | A particle filtering-based framework for real-time fault diagnosis and failure prognosis in a turbine engine[END_REF]] a methodology using a PF for real-time failure prognosis in two stages. The first one consists in projecting the current particle population in time without any new measurements. The state equation was simply propagated by considering particles with unchanged weights. Then, during the second stage, the RUL probability density function (pdf) was computed by combining the weights of the predicted trajectories and the specifications for the hazard zone containing the failure threshold. Fatigue crack growth prognosis on a turbine engine blade using real data was used to demonstrate the efficiency of the proposed method and to show the effects of model uncertainty on the algorithm performance.

Later, the authors extended their method in [START_REF] Orchard | A particle-filtering approach for on-line fault diagnosis and failure prognosis[END_REF] by proposing a new approach to predict the future degradation state without new measurements. In order to take uncertainty in future state transitions into account, they suggested to use Epanechnikov kernels to approximate the state pdf and then to resample this predicted pdf at each future time instants. Thus, information about the future state pdf is given by the position of the particles and not their weights anymore. This procedure allowed them to reduce the impact of model errors. [START_REF] Zio | Particle filtering prognostic estimation of the remaining useful life of nonlinear components[END_REF] proposed a PF-based method to estimate accurately the RUL of a system. Biased estimators using credibility weights (interval weights, likelihood weights and RMSE weights) to measure the credibility of the RUL estimates were introduced. The addition of these estimators to the PF algorithm allowed the authors to obtain more precise RUL predictions. A fatigue crack growth prognosis case study based on the Paris'law was used to demonstrate the efficiency of the method.

[ Daigle and Goebel 2011a] presented a method for pneumatic valves RUL prediction with a PF. First, a detailed physics-based nominal model of a pneumatic valve was established. Then, the relevant damage mechanisms were modeled, and the model parameters that evolve with the damage were selected. Finally, a PF was used for joint state-parameter estimation. A sequential importance resampling (SIR) PF with systematic resampling was chosen, and the parameter evolution was described by a Gaussian random walk. In order to compute the RUL pdf, the state 1.1. Prognosis approaches equation was propagated with unchanged particle weights until the failure threshold was reached. The proposed methodology was demonstrated with historical pneumatic valve data from the refueling system, and gave satisfactory results. However, the authors have assumed a single known future input trajectory, which is not the case in general. Moreover, the case where only a single damage mode is progressing was considered.

This latter issue was treated by the authors in [Daigle and Goebel 2011b]. Indeed, they developed a PF-based method for multiple damage progression paths prognosis. In addition, they developed a variance control technique to maintain an uncertainty bound around the parameters that are jointly estimated with the state in order to reduce prediction uncertainty. The methodology was demonstrated through simulation-based experiments aiming at performing multiple damage prognosis in a centrifugal pump.

In [START_REF] Olivares | Particle-filtering-based prognosis framework for energy storage devices with a statistical characterization of state-of-health regeneration phenomena[END_REF], a PF-based approach for simultaneous SOH and RUL estimation in energy storage devices was proposed. The authors implemented a PF-based detection and prognosis algorithm [START_REF] Orchard | Risk-sensitive particle-filtering-based prognosis framework for estimation of remaining useful life in energy storage devices[END_REF]] capable of detecting and including the effect of self-regeneration phenomena in long-term predictions. Experimental data from accelerated degradation tests were used to demonstrate the efficiency of the proposed approach. [START_REF] Miao | Remaining useful life prediction of lithium-ion battery with unscented particle filter technique[END_REF] proposed a method to tackle particle impoverishment and to obtain a more accurate state pdf estimate by using an unscented particle filter (UPF). The choice of the proposal distribution which allows the computation of the posterior probability when a new measurement is available can have significant impact on the accuracy and computational efficiency of the PF. Usually, the prior distribution is chosen as proposal distribution. In this paper, the authors considered the output of an unscented Kalman filter (UKF) as proposal distribution. Thus, the UPF algorithm is based on two steps: (i) application of UKF algorithm to compute the proposal distribution and (ii) application of standard PF algorithm to obtain the results. The proposed methodology was tested with battery capacity data from the Center for Advanced Life Cycle Engineering (CALCE) of University of Maryland. The results showed that the prediction error with the UPF-based method is lower than the one obtained with a classical PF.

As the use of the PF algorithm for on-line failure prognosis is gaining more and more popularity, a complete Matlab tutorial for model-based prognosis using a PF was given in [An, Choi, and Kim 2013b]. The authors detailed the PF-based prognosis framework and proposed two case studies: battery SOH prognosis and fatigue crack growth prognosis. The authors aimed at helping the beginners in prognosis by providing a standard of PF-based prognosis.

In [START_REF] Daroogheh | A novel particle filter parameter prediction scheme for failure prognosis[END_REF], a method to propagate the degradation equation in future time instants without any new measurements while updating Chapter 1. Prognosis and uncertainty propagation the particle weights and perform particle resampling is proposed. For this purpose, future observation profiles are generated using a fixed-lag dynamic linear model (DLM) structure [START_REF] West | Mixture models, monte carlo, bayesian updating, and dynamic models[END_REF]] with adaptive length moving window. More precisely, the fixed-lag DLM model aims at fitting the historical measurements in each time window to a linear autoregressive moving average (ARMA) model which is then used to forecast the future observations. The DLM is updated every time new measurements are collected. In order to highlight the benefits of such strategy, it was compared to PF-based RUL prediction approach with invariant particle weights through a gas turbine failure prognosis example. The resulting percentage RMSE showed the effectiveness of the proposed methodology. [START_REF] Hu | A particle filtering and kernel smoothingbased approach for new design component prognostics[END_REF] proposed a solution to deal with particle parameter impoverishment when joint-state parameter estimation is performed, by introducing a kernel smoothing (KS) technique. Usually, instead of considering unchanged parameter distributions, a Gaussian random walk is added to the parameters evolution equation in order to avoid particle impoverishment and obtain more accurate parameter estimates. However, such a strategy tends to increase the parameters covariance and thus decrease the precision of the estimated values. Moreover, the variance of this artificial noise must be set up, which is in most cases done with a trial and error procedure. When complete degradation trajectories are not available, as in the presented work, this strategy can not be adopted. Therefore, to bypass these two issues; the authors used a KS algorithm which consists in two steps: (i) shrinkage to force the particle parameters to move toward their estimated values, and (ii) perturbation to increase their variance and thus maintain a desired particle parameters covariance. For RUL prediction, the degradation equation was propagated in time with invariant particle weights. The proposed algorithm was applied to the prognosis of lithium-ion battery degradation and its performance was evaluated in terms of precision, accuracy and steadiness indexes. The results showed that the use of the KS algorithm lead to better results than the Gaussian random walk strategy or the classical PF algorithm with unchanged parameter distributions. [START_REF] Wang | Remaining useful life prediction of lithium-ion batteries based on spherical cubature particle filter[END_REF] developed a strategy other than the one presented in [START_REF] Miao | Remaining useful life prediction of lithium-ion battery with unscented particle filter technique[END_REF] to improve the choice of the proposal distribution. This time, a spherical cubature PF (SCPF) was used for lithium-ion battery degradation prognosis, i.e a PF which takes the results of a SC Kalman filter (SCKF) as proposal distribution. The method based on spherical cubature [START_REF] Särkkä | Bayesian filtering and smoothing[END_REF]] is able to give more weights to the recent measured data when updating the proposal distribution. Capacity degradation from [START_REF] Park | Dual features functional support vector machines for fault detection of rechargeable batteries[END_REF] were used to compare the proposed method to classical PF. The results allowed the authors to conclude that RUL predictions are improved with this method.

More recently, [START_REF] Mo | A remaining useful life prediction approach for lithium-ion batteries using kalman filter and an improved particle filter[END_REF] proposed a PF-based method that handles both high observation noise and particle impoverishment during the estimation step. The authors integrated a Kalman filter (KF) to the basic PF 1.1. Prognosis approaches algorithm to obtain an optimal state pdf estimate in the presence of high observation noise. In order to deal with particle impoverishment, they used the particle swarm optimization (PSO) algorithm. PSO forces the particles to move toward optimal state space regions before the update of the particle weights when a new measurement arrives, then resampling is performed. This PSO-PF based approach was introduced in [START_REF] Zhang | Particle filter based on pso[END_REF]. Battery degradation datasets from NASA AMES Center [START_REF] Saha | Battery data set[END_REF] were used. Two batteries were first used as training data to estimate the degradation model parameters. Then, a third dataset was used to perform RUL prediction by propagating the degradation equation in time with unchanged particle weights. Results with the proposed KF-PSO-PF, the PSO-PF and the classical PF were compared. The calculated RMSE showed the KF-PSO-PF gives more precise and accurate RUL predictions than the two other methods.

The same authors developed later in [START_REF] Yu | Remaining useful life prediction for lithium-ion batteries using a quantum particle swarm optimization-based particle filter[END_REF]] a quantum PSO-based PF (QPSO-PF) to predict the RUL of lithium-ion batteries. A QPSO algorithm [START_REF] Sun | Particle swarm optimization with particles having quantum behavior[END_REF] was used instead of a PSO algorithm, and no Kalman filtering algorithm was included. The QPSO algorithm aims at obtaining a global optimum at a faster convergence. Moreover, fewer particles and tuning parameters are required. The proposed QPSO-PF method was tested with the same NASA AMES Center battery dataset, and joint state-parameter estimation with a Gaussian random walk to describe the parameter evolutions was performed. The results showed that higher precision and accuracy than the PSO-PF and classical PF were obtained. Moreover, the authors concluded that the real-time computation could be accelerated with such strategy. However, further investigations should help to learn much about the control and determination of the parameter involved in the QPSO algorithm.

There exist several papers that present PF-based model-based prognosis algorithms, and a more complete review can be found in [START_REF] Jouin | Particle filter-based prognostics: Review, discussion and perspectives[END_REF] as it is not the purpose of this chapter. However, regarding the PF-based prognosis methods that have been reviewed above, some conclusions can be drawn:

• Particle filter is the most widely used tool in model-based prognosis;

• Most of the case studies concern battery degradation and fatigue crack growth prognosis. Indeed, physics-based degradation models were widely developed in these two areas;

• The more accurate the state and parameter estimations are at the prediction time, the more accurate the RUL prediction will be. Various solutions have been provided to obtain the most accurate pdf estimates as possible:

-Optimize the proposal distribution (UPF, SCPF, etc.); -Choose a suitable particle parameters evolution equation (Gaussian random walk, kernel smoothing, etc.) when joint state-parameter estimation is involved; -Move the particles into an optimal state space region before the update and resampling steps (PSO, QPSO, etc.);

• Most of the time, the particle weights are considered unchanged when propagating the degradation equation in future time instants. Even if such a strategy can help to obtain correct results, RUL predictions can be improved by adopting a particular method when predicting the future degradation state without any new measurements:

-Use Epanechnikov kernels to approximate the state pdf and then resample the predicted pdf in the future time instants; -Generate future observations profiles with available historical measurements and time series forecasting methods.

Although the PF was widely used in model-based prognosis mainly because of its ability to deal with nonlinear degradation models, other approaches based on Bayesian recursive estimation have, however, been developed. Kalman filter. The use of KF-based methods for RUL prognosis is limited due to the fact that they cannot deal with high nonliearities and non-Gaussian noises. [START_REF] Celaya | Uncertainty representation and interpretation in model-based prognostics algorithms based on kalman filter estimation[END_REF] discussed the ability of the Kalman filter (KF) to represent uncertainty in model-based prognosis of electronic components. [START_REF] Daigle | A comparison of filter-based approaches for model-based prognostics[END_REF] compared three nonlinear filtering-based techniques, namely Daum filter, unscented KF and PF for damage prognosis in a centrifugal pump. [START_REF] Lim | Switching kalman filter for failure prognostic[END_REF]] used a switching KF for degradation prognosis of a helicopter gearbox bearing. It was assumed that the degradation could evolve through time and each degradation process was therefore modeled with a KF. The SKF was used to infer the underlying process and apply the most probable filter for RUL prediction.

In [START_REF] Ompusunggu | Kalman-filtering-based prognostics for automatic transmission clutches[END_REF], an extended KF (EKF) was used to predict the RUL of wet friction clutches used in vehicles equipped with automatic transmissions. [Bressel, Hilairet, Hissel, and Bouamama 2016b] performed RUL prognosis of a Proton Exchange Membrane Fuel Cell with an EKF.

In [START_REF] Wang | Predictive airframe maintenance strategies using model-based prognostics[END_REF], the authors proposed a fatigue crack growth prognosis method that combines an EKF for current degradation state estimation, and first-order perturbation (FOP) [START_REF] Yiwei | A cost driven predictive maintenance policy for structural airframe maintenance[END_REF] to predict the future damage distribution.

In the previous studies, the authors concluded that the KF-based methods could give satisfactory prognosis results. However, this is true as long as the models remain close to linear conditions, and only a few number of unknown model parameters are involved. Moreover, the Gaussian noise assumption may not always be accurate.

Hidden damage

Model-based prognosis methods were developed for the case where a hidden damage process is evolving in the system. The tracking of the degradation state relies on the assumption that the damage process is growing slowly in a directly observable dynamic system with a fast behavior.

In [START_REF] Chelidze | A dynamical systems approach to failure prognosis[END_REF], a procedure for hidden damage tracking and RUL prediction is presented. First, based on phase space reconstruction data from the system, a damage tracking algorithm is used to compute a damage tracking metric. This latter is then used as input to a nonlinear filtering algorithm (unscented filter) together with the damage evolution model in order to forecast the future degradation state and finally predict the RUL. The proposed algorithm was applied to an electromechanical oscillator, and the battery discharge was considered to be a hidden damage process.

In [START_REF] Luo | Model-based prognostic techniques applied to a suspension system[END_REF], another model-based prognosis method for hidden damage tracking is proposed. First, a degradation model of the system combining the fast-time system behavior and the slow-time damage evolution is build. Then, using this model and Monte Carlo simulations under different load conditions, a prognosis model combining multiple models is built. Finally, an IMM filter is used on-line to track the degradation measure which is a hidden variable, and the RUL of the system is estimated. They illustrated their method on an automotive suspension system subject to a hidden damage evolution which is a crack growing in the suspension spring.

Another kind of model-based prognosis methods that are used in the context where the degradation of the system is not directly measured are based on unknown input observers (UIOs). In this case, the degradation is considered as a perturbation that affects the system behavior, and is referred to as an unknown input. Therefore model-based prognosis methods based on UIOs aim at reconstructing the degradation signal and predict the RUL of the system. In [Gucik-Derigny, Outbib, and Ouladsine 2011a], a methodology to perform RUL prognosis with UIOs is introduced. A nonlinear unknown input high-gain observer for state and unknown input estimation was synthesized, and then used to predict the RUL of an electromechanical oscillator. The same authors compared the ability of two types of UIO to perform RUL prognosis in [Gucik-Derigny, Outbib, and Ouladsine 2016a], namely unknown input high-gain observer and unknown input sliding-mode observer. Both methods were tested on an electromechanical system with one battery then with two batteries. The results showed that the high-gain UIO gives more accurate RUL predictions than the sliding-mode UIO.

Reliability analysis methods

In the context of reliability analysis, dedicated methods are used to estimate a failure probability which characterizes the safe/unsafe behavior of a system regarding the uncertainties affecting the input variables. Recently, some researchers have extended these methods to model-based prognosis. [START_REF] Sankararaman | Analytical algorithms to quantify the uncertainty in remaining useful life prediction[END_REF] investigated the use of three reliability analysis methods within a state-space framework and then extended them to the problem of RUL prognosis. These methods are the first-order second moment method (FOSM) [START_REF] Madsen | Methods of Structural Safety[END_REF], the first-order reliability method (FORM) [START_REF] Lemaire | Structural reliability[END_REF]] and the Inverse FORM [START_REF] Der Kiureghian | Inverse reliability problem[END_REF]. They are referred to as analytical methods, as opposed to sampling-based methods such as adaptive sampling, Latin hypercube sampling and so on. Indeed, they are able to quantify uncertainty in RUL predictions with a reduced number of calculations. FOSM can be used to calculate the first two moments (mean and variance) of the RUL distribution. However, it cannot be used to calculate tail probabilities as it does not estimate the type of probability distribution of the RUL. FORM aims at computing any value of the cumulative distribution function (CDF) at a particular RUL value. In other words, it calculates the probability that the RUL is smaller than a given value. However, it is difficult to guess this RUL value as the RUL distribution is not known in advance. Rather, the authors suggested to use the Inverse FORM algorithm which allows the calculation of the RUL value corresponding to a given probability level. Inverse FORM can be used to compute only the desired RUL CDF values (median, probability bounds) or the entire RUL CDF by applying it as many times as required. The efficiency of the proposed Inverse FORM algorithm was demonstrated by calculating the median and the 95% probability bounds of the RUL of a lithium-ion battery. The results were compared to those obtained with 1000 Monte Carlo simulations (MCS) and showed that maximum difference between the two CDFs is less than 0.5%. Moreover, the computational cost of the Inverse FORM is about 10% of the cost of MCS.

The same authors have discussed in more details the practical challenges involved in RUL prediction with inverse FORM in [Sankararaman, Daigle, and Goebel 2014b]. They have especially included a time-invariant process noise uncertainty which was not considered in their previous work. The same lithium-ion battery model was used to illustrate their methodology which gave accurate RUL estimated in comparison with MCS.

In [Bressel, Hilairet, Hissel, and Bouamama 2016a], the Inverse FORM was applied for RUL prognosis of Proton Exchange Membrane Fuel Cell.

More recently, [START_REF] Chiachío | A new algorithm for prognostics using subset simulation[END_REF] proposed an adaptation of a sampling-based reliability analysis method to RUL prognosis. Their strategy combines a particle filter for state estimation until the prediction time and then, the use of Subset Simulation (SubSim) algorithm [START_REF] Au | Estimation of small failure probabilities in high dimensions by subset simulation[END_REF] to propagate uncertainties from the current state to the failure threshold. The authors' starting point was the fact the uncertainty affecting multi-step ahead predictions with the particle filter due to particle approximation is propagated in time and thus increases the uncertainty in RUL prediction. The longer this multi-step ahead prediction is, the larger the RUL uncertainty gets. Thus, the authors proposed to consider intermediate failure domains between the current state and the failure threshold. The SubSim algorithm was then applied to adaptively draw samples over the nested sequence of intermediate failure domains until the failure threshold was reached. The basic idea of SubSim is, in the reliability context, to efficiently estimate a small failure probability (thus, difficult to estimate) by considering nested failure events associated to a product of stronger probabilities (thus, easier to estimate). The authors applied their methodology to fatigue damage prognosis of carbon-fibre composite coupons. They obtained satisfactory results, but highlighted that this method is especially efficient is the case of rare-event prediction.

Conclusion

Examples of data-driven and model-based prognosis approaches that can be found in the literature were given in this section. The main criteria that determines the use of one or other of these approaches is the availability or not of data and a degradation model. The main drawback of data-driven approaches is the need of a significant amount of data to train the prognosis model that will be used for on-line RUL prognosis. As for model-based prognosis methods, a degradation model is not always available depending on the case study.

In this thesis, the case where a physics-based model of the degradation is available is considered. Thus, in the following, the principle of model-based prognosis is presented, as well as the main challenges involved in the development of such methods.

Model-based prognosis

Model-based prognosis methods rely on the availability of a mathematical model of the degradation. The main advantage of this approach is its precision since the predictions are achieved based on a mathematical model of the degradation. Indeed, the models include the operational conditions of the system and therefore the effect of their variability on the degradation evolution.

In this section, the general principle of model-based prognosis methods and the challenges associated to their development, mainly concerning uncertainty propagation, are presented. Finally, the performance evaluation metrics that will be considered all along this manuscript are given.

Principle of model-based prognosis

The central idea of model-based prognosis is to use a dynamic mathematical model that describes the evolution of a degradation within a system or a component. Usually, a discrete state-space representation of the model is adopted for simulation purpose.

The RUL calculation procedure can be divided into two main steps: (i) current degradation state estimation and (ii) future degradation state prediction for RUL calculation. The model-bases process, which is described in what follows, is illustrated in Figure 1.1.

Current degradation state estimation: a filtering problem

The first step consists in estimating the current degradation state using the model and the measurements from sensors. This procedure is referred to as state estimation and is addressed with filtering, which is based on the use of two equations, namely the state equation and the measurement equation.

The state equation is the discretized mathematical model of the degradation and therefore describes its evolution according to the following expression:

x k = f (x k-1 , θ k-1 , u k , w k ) (1.1)
where x ∈ R n denotes the degradation state, u ∈ R q is the input vector, θ contains the model parameter vector that are possibly unknown and k ∈ N is a discrete time step. The function f describes the evolution of the state which may be nonlinear, and w is the process noise which represents the model uncertainties.

As for the measurement equation, it relates the measurements to the state according to:

y k = h(x k , v k ) (1.2)
where y ∈ R p is the measurement vector and h describes the evolution of the measurements over time and may be nonlinear. The variable v is the measurement noise which represents the measurements uncertainty.

There exist various filtering techniques for state estimation. Based on the literature review in the previous section, the most used one is the particle filter because of its ability to handle nonlinear systems and non-Gaussian noises. If linear systems and Gaussian noises are involved, the Kalman filter gives optimal solutions, and the

Model-based prognosis

extended Kalman filter can be used for nonlinear systems with Gaussian noises. Thus, using these two equations and appropriate filtering techniques, the current degradation state is estimated from the initial time step k 0 until the so called prediction time step k p . The time interval between k 0 and k p is referred to as the observation interval.

Future degradation state prediction and RUL computation: an uncertainty propagation problem

From the prediction time k p , the future degradation state must be predicted until the critical degradation state, i.e the failure threshold, is reached. When the failure threshold is finally crossed, the failure time k f is obtained. The time interval between the prediction time k p and the failure time k f is called the prediction interval. Finally, when the failure time k f is obtained, the RUL at time step k p can be computed as:

RU L(k p ) = k f -k p (1.3)
The difficulty in this step is that the degradation state prediction is performed without any new measurements. Moreover, the forecasting process is naturally affected by uncertainty. There is therefore a need to use uncertainty propagation techniques to take input uncertainty into account in order to compute a measure of the uncertainty in the RUL prediction, in the form of a probability distribution function (pdf) or an interval. As new measurements are collected, predictions are improved via model parameters and degradation state estimates updating. Thus, the uncertainties are reduced over time.

The main uncertainty sources and their impact on the prognosis process will be discussed in the next section. Moreover, uncertainty propagation techniques to tackle RUL prediction problem will be presented.

Uncertainty in model-based prognosis

Because of economic and human safety challenges, the prediction of the RUL must be as reliable as possible to allow risk-based decisions. For this purpose, prognosis methods should account for the different uncertainty sources that affect both the current and future degradation state, and thus provide a measure of the uncertainty affecting the predicted RUL [START_REF] Tang | Methodologies for uncertainty management in prognostics[END_REF][START_REF] Baraldi | Investigation of uncertainty treatment capability of model-based and data-driven prognostic methods using simulated data[END_REF]].

The main types of uncertainty that inevitably influence RUL prediction are measurement uncertainty, model uncertainty and future uncertainty [START_REF] Gu | Uncertainty assessment of prognostics of electronics subject to random vibration[END_REF]]. Measurement uncertainty is due to sensor inaccuracy, model uncertainty is characterized by the difference between the degradation model and its real behavior, and future uncertainty is caused by unknown or uncertain environmental and operating conditions. Moreover, another source of uncertainty that is not negligible is the current state uncertainty at the prediction time k p as the future predictions depend on this current state estimation.

Model uncertainty.

A physics-based model of the degradation is used to estimate and forecast the degradation state. However, it is impossible to accurately represent a physical phenomena, the model that are used are therefore approximations. Usually, to take model uncertainty into account, a noise is added to the state equation. Model uncertainty can also include parameter uncertainty. Indeed, the degradation model can contain unknown or uncertain parameters to be estimated. Model uncertainty can be quantified as long as measurement are available until the prediction time k p , it is however more challenging to know their values in the prediction interval.

Measurement uncertainty. The measured output y k at time step k ≤ k p differs from the real value of the measured quantity due to sensor inaccuracy. It is difficult to quantify measurement uncertainty, however, in some cases, they can be known from the calibration of the measurement devices. Measurement uncertainty are generally taken into account by adding a noise to the measurement equation.

Current state uncertainty.

The predicted state at future time instants and thus the predicted RUL highly depend on the state estimate x kp at the prediction time. Indeed, the more accurate x kp is, the more accurate RUL predictions will be. The uncertainty in x kp is mainly due to the model and measurement uncertainty discussed above, and can be handled by choosing an appropriate estimation technique. It is important to have a precise estimation of x kp as the uncertainty at time k p will be propagated in time through the state equation.

Future uncertainty.

The main challenge in prognosis is to deal with future uncertainty. Indeed, the future loading and future operating conditions of the system are unknown or uncertain, making future state and RUL predictions very uncertain. Thus, taking these uncertainties into account while propagating the state equation until the failure threshold will improve RUL predictions.

To conclude, the choice of appropriate filtering and uncertainty propagation techniques is the core of prognosis, whose goal is to combine the effect of the uncertainty sources in order to measure the uncertainty in RUL prediction.

Usually, model-based prognosis methods represent and propagate uncertainty in a probabilistic framework [START_REF] Saha | Uncertainty management for diagnostics and prognostics of batteries using bayesian techniques[END_REF], in order to compute the RUL PDF. As stated in section 1.1.2, the main filtering techniques that are used to estimate the current degradation state with the measurements while taking uncertainty into account are based on particle filters and Kalman filters. As for the uncertainty propagation techniques, improved versions of PF-based methods were largely employed. However, the main drawback when choosing a probabilistic representation is that assumptions on the form of the PDF are made a priori, which can be inaccurate or misleading.

In order to avoid this problem, uncertainty representation within the set-membership framework could be adopted. In this case, any stochastic assumptions are required and RUL uncertainty is included in an interval. The main advantage of using setmembership techniques is that the predicted interval is guaranteed to contain the RUL. However, such approaches are limited to the cases where only a few number of parameters to be estimated is involved. In the field of model-based prognosis, set-membership methodologies are not yet very widespread. [Gucik-Derigny, Outbib, Ouladsine, and Basualdo 2011b] proposed a continuous interval observer for state and unknown input estimation. The reconstructed unknown input was the damage variable, but the authors did not perform damage prognosis.

In this thesis, stochastic filtering and unknown input observers are used to estimate the current degradation state with the measurements. Uncertainty propagation methods using stochastic filtering methods, reliability analysis methods and intervalbased methods are employed. The main requirements regarding the developed prognosis methods concern uncertainty quantification and on-line implementation.

Chapter 1. Prognosis and uncertainty propagation

Performance evaluation metrics

It is essential to evaluate the ability of prognosis methods to perform RUL prognosis. For this purpose, performance evaluation metrics can be used.

There is no strict agreement about which appropriate and acceptable set of metrics should be used in prognosis applications. However, it is widely admitted that accuracy and precision indicators are relevant to examine the performance of prognosis algorithms [START_REF] Edwards | Impact of input uncertainty on failure prognostic algorithms: Extending the remaining useful life of nonlinear systems[END_REF]. These metrics derive from the prediction error performed at time

k p which is expressed as e(k p ) = RU L true (k p ) -RU L(k p ) where RU L true (k p ) is the reference RUL value at time k p .
The three considered metrics are presented in the following:

• Accuracy is a measure of the degree of closeness of predicted failure time k f to the actual failure time k f T rue . This metric provides an exponential weight of the errors in RUL predictions over several experiments. The accuracy of a prognosis algorithm at a specific prediction time k p is defined here as in [START_REF] Vachtsevanos | Intelligent fault diagnosis and prognosis for engineering systems[END_REF]:

A(t p ) = 1 N s Ns ∑ n=1 exp ( - |e n (k p )| RU L(k p ) ) (1.4)
where N s is the number of experiments and e n (k p ) is the prediction error of the n th experiment. The range of the accuracy is between 0 and 1, where 1 gives the best accuracy.

• Precision is a measure of the narrowness of the interval in which the RUL predictions fall and is expressed as [START_REF] Butler | Prognostic Algorithms for Condition Monitoring and Remaining Useful Life Estimation[END_REF]]:

P (k p ) = exp ( - R R 0 ) (1.5)
where R is the width of the confidence interval of the prediction given by R = 2 × 3σ RU L where σ RU L is the standard deviation of the RUL pdf. R 0 is a normalizing factor. The precision value varies between 0 and 1, where 1 reflects the highest precision.

• Timeliness indicates the relative position of the predicted RUL pdf along the time axis with respect to the occurrence of the actual failure event. There are three cases (Fig. 1.2): (a) the failure occurs after the predicted failure time k f , (b) the failure occurs at the same time as the predicted failure time, and finally, (c) the failure occurs earlier than predicted. This last case must be absolutely avoided. To compute the timeliness metric, the following function is used [START_REF] Saxena | Metrics for evaluating performance of prognostic techniques[END_REF]:

1.4. Conclusion T = 1 N Ns ∑ n=1 T n (1.6) T n = ⎧ ⎨ ⎩ exp(-en(kp) R min ) -1 , if e n (k p ) ≤ 0 exp( en(kp) Rmax ) -1 , if e n (k p ) > 0 (1.7) [R min , R max ]

Conclusion

In this chapter, an overview of the different prognosis approaches that can be found in the literature was given. Prognosis approaches can be classified into data-driven approaches, model-based approaches and hybrid approaches which combines both. Data-driven methods mainly use statistical techniques (HMM, RVM/SVM, etc.) and artificial intelligence techniques (mainly NN). Model-based prognosis methods are most of the time based on particle filters which are powerful tools able to properly deal with non-Gaussian noises and nonlinear systems.

Model-based prognosis is divided in two main steps: (i) current degradation state estimation and (ii) future degradation state prediction and RUL computation. The first step is a filtering problem as it aims at estimating the current degradation state using measurements. The second step is an uncertainty propagation problem as the current degradation state is propagated in time until the failure threshold any new measurements. Thus, the main challenge while developing model-based prognosis methods is to include the different uncertainty sources in order to compute a measure of the uncertainty affecting the predicted RUL. The main uncertainty sources are model uncertainty, measurement uncertainty and future uncertainty (loading, operating conditions, etc.).

Chapter 2 Filtering methods for prognosis with stochastic assumptions

The primary focus of this chapter is to present in more detail the principle of model-based prognosis methods through an application example. As it was highlighted in section 1.2.1 of Chapter 1, there are two main steps in prognosis: (i) current degradation state estimation, (ii) future degradation state prediction and RUL computation.

However, these two steps are subject to uncertainty. Thus, the main challenge when developing a prognosis method is to take model, measurement and future loading uncertainty into account. In this chapter, only modeling and measurement uncertainty are addressed for clarity and ease of explanation.

To address this challenge, two stochastic filtering-based methods, namely an extended Kalman filter (EKF) and a particle filter (PF), are investigated. For the sake of illustration, the prognosis algorithms are applied to a nonlinear fatigue crack growth problem based on the Paris' law. The model used in this chapter are from [An, Choi, and Kim 2013b]. In addition, the performance evaluation metrics presented in Chapter 1 are used to compare the efficiency of the EKF and the PF in terms of accuracy, precision and timeliness. This first study has led to the publication of an international conference paper [START_REF] Robinson | Model-based prognosis algorithms with uncertainty propagation: application to fatigue crack growth[END_REF].

Recursive Bayesian filtering

In order to solve the dynamic degradation state estimation problem, two mathematical equations are needed: a first one describing the evolution of the degradation state (1.1) and a second one relating the state and the noisy measurements (1.2). Here, the input vector u is not considered, the following discrete state-space representation Chapter 2. Filtering methods for prognosis with stochastic assumptions is therefore adopted:

x k = f (x k-1 , θ k-1 , w k ) (2.1) y k = h(x k , v k ) (2.2)
As a reminder, x ∈ R n denotes the degradation state, y ∈ R p is the measured outputs and k ∈ N is a discrete time step. The functions f and h describe respectively the nonlinear evolution of the state and the measurements over time. The variables w and v are respectively the process and measurement noises which represent the model and measurements uncertainties.

Stochastic methods

The objective of stochastic filtering is to estimate the pdf p(x k |y 0:k ) which gives statistical information about the degradation state x k , based on the set of all measurements y

0:k = Y k = {y 0 , • • • , y k }.
The state equation (2.1) characterizes the state transition pdf p(x k+1 |x k ), whereas the measurement equation (2.2) describes the pdf p(y k |x k ) which is further related to the measurement noise model.

Concerning the unknown parameter vector θ, it is appended to the state vector x ∈ R n to form an augmented state vector X = Given p(X 0 ), p(X k+1 |X k ) and p(y k |X k ), the recursive Bayesian estimation presented below is used to solve the stochastic filtering problem.

Recursive Bayesian estimation principle

In Bayesian theory, the uncertainties are treated as random variables. Moreover, a recursive filtering approach means that received data can be processed sequentially rather than as a batch so that it is not necessary to store the complete data set nor to reprocess existing data if a new measurement becomes available.

First of all, it is assumed that: (i) the state vector is a first-order Markov process such that p(X k |X 0:k-1 ) = p(X k |X k-1 ) and (ii) the observations are independent of the states. Combining these two assumptions and the Bayes rules, the following equation is obtained [START_REF] Chen | Bayesian filtering: From Kalman filters to particle filters, and beyond[END_REF]:

p(X k |Y k ) = p(y k |X k )p(X k |Y k-1 ) p(y k |Y k-1 ) (2.3)
The posterior density p(X k |Y k ) is defined through the combination of three terms:

• The prior density p(X k |Y k-1 ) which is the prediction density of the state at

Extended Kalman filter

time k obtained via the Chapman-Kolmogorov equation:

p(X k |Y k-1 ) = ∫ p(X k |x k-1 )p(X k-1 |Y k-1 )dX k-1 (2.4)
where p(X k |x k-1 ) is the state transition density defined by the state equation (2.1);

• The Likelihood density p(y k |X k ) which is defined by equation (2.2);

• A normalizing constant p(y k |Y k-1
) which depends on the likelihood and the prior such that:

p(y k |Y k-1 ) = ∫ p(y k |X k )p(X k |Y k-1 )dX k . (2.5)
Bayesian filtering handles the computation or approximation of these three terms to deduce the pdf of the degradation state p(X k |Y k ). It is based on two steps: prediction and update. First, the required pdf p(X k-1 |Y k-1 ) is supposed to be available. During the prediction step, using the previous pdf and the system model (2.1), the prior p(X k |Y k-1 ) is approximated with equation (2.4). Then comes the update step at time k when a measurement y k becomes available, the likelihood pdf p(y k |X k ) is obtained with the measurement equation (2.2). Finally, the posterior density p(X k |Y k ) is deduced from equation (2.3).

Based on this Bayesian filtering algorithm, many types of filters have been developed [START_REF] Chen | Bayesian filtering: From Kalman filters to particle filters, and beyond[END_REF]]. In this chapter, an EKF and a PF are used, since they can handle nonlinear dynamical systems.

Extended Kalman filter

The EKF is an extended version of the original Kalman filter (KF) [START_REF] Kalman | New results in linear filtering and prediction theory[END_REF] developed for nonlinear systems. In this filtering approach, the state pdf p(X k |Y k ) is approximated by a Gaussian distribution. To achieve this, the nonlinear degradation model is linearized around the last predicted degradation state estimate and the conventional KF algorithm is applied to the linearized dynamics. To linearize the nonlinear functions f and h which are assumed to be differentiable, their respective Jacobian matrices F and H are computed at each time step with the predicted degradation state:

F k = ∂f ∂x ⏐ ⏐ ⏐ ⏐ ⏐ Xk-1
(2.6)

H k = ∂h ∂x ⏐ ⏐ ⏐ ⏐ ⏐ Xk

Current degradation state estimation with the EKF

To solve the prognosis problem, the first step consists in the estimation of the current degradation state, while measurements are available. Hence, the classical EKF algorithm [START_REF] Jazwinski | Stochastic Processes and Filtering Theory[END_REF]] expressed as follows is used:

Prediction step Xk|k-1 = f ( Xk-1|k-1 , w k = 0)
(2.7)

P k|k-1 = F k P k-1|k-1 F T k + Q k Update step ŷk = y k -H k Xk|k-1 (2.8) K k = P k|k-1 H T k (H k P k|k-1 H T k + R k ) P k|k = (I -K k H k )P k|k-1 Xk|k = Xk|k-1 + K k ŷk
where P , Q and R are covariance matrices, respectively of the estimation error, the process noise and the measurements noise. K is the Kalman gain.

Future degradation state estimation with the EKF

In this part, the previous estimation Xk p-1 |k p-1 is used as the initial degradation state while the prediction step remains the same. The update step is changed because the innovation term ŷk = y k -H k Xk|k-1 is no more available. Instead, the degradation state at step k ∈ {k p , • • • , k pf } is updated with the state transition model (2.1) using the standard deviation of the previous degradation state to approximate the distribution of the noise w. The algorithm for the estimation of the future degradation step from step k p to step k pf is given by the following equations:

Prediction step Xk|k-1 = f ( Xk-1|k-1 , w k = 0)
(2.9)

P k|k-1 = F k P k-1|k-1 F T k + Q Update step P k|k = P k|k-1 (2.10) ŵk ∼ N (0, P k|k ) Xk|k = f ( Xk|k-1 , ŵk )

Particle filter

The EKF generally provides satisfying results and is easy to implement. But as the desired state pdf is approximated by a Gaussian distribution, it may have significant deviation from the true distribution causing divergence in the case where the degradation model is highly nonlinear. In order to deal with more complex degradation models, particle filters can be used.

Particle filter

In the particle filter approach, the state pdf at time instant k is approximated by a set of

N part particles {X i k } Npart i=1
representing points in the state space, and a set of associated weights

{ω i k } Npart i=1
denoting discrete probability masses:

p(X k |y 0:k ) ≈ Npart ∑ i=1 ω i k δ(X k -X i k ) with Npart ∑ i=1 ω i k = 1 (2.11)
where δ is the Dirac delta function.

There exist several PF algorithms [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF]. One of the most used is the sequential importance resampling (SIR) particle filter. It is based on three main steps that are prediction, update and re-sampling, which are described in the following and represented in Figure 2.1:

0. Initialization at time instant k = 1 • Draw particles X i 0 ∼ p(X 0 ) • Compute the initial weights ω i k = 1

Npart

• Repeat the three following steps until the last available measurement 1. Prediction which are realizations of the predicted pdf p(X k |y 0:k-1 ).

Update

• Each sampled particle is assigned a weight based on the likelihood p(y k |X k ):

ω i k = ω i k-1 p(y k |X i k-1 ) = ω i k-1 p(y k |X i k )p(X i k |X i k-1 ) p(X i k |X i k-1 , y k )
(2.12)

• Normalize the weights:

ω i k = ω i k ⎛ ⎝ Npart ∑ i=1 ω i k ⎞ ⎠ -1
(2.13) Particles with small weights are eliminated so that the computational efforts are concentrated in those having large ones.

• Re-sampling condition: if the effective sample size N ef f is under some threshold N th , a re-sampling procedure is done. An estimate of

N ef f is Neff = ⎛ ⎝ Npart ∑ i=1 (ω i k ) 2 ⎞ ⎠ -1
(2.14)

In this work, systematic resampling is performed.

• Using the inverse cumulative distribution function method [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF] and the current set

{X k } Npart i=1 , a new set { Xk } Npart i=1
is drawn to replace the current one. Finally, with ωi k = N -1 part , the state is given by:

Xi k = Npart ∑ i=1 ωi k Xk i (2.15)

Current degradation state estimation with the PF

This basic SIR particle filter algorithm is applied during the first step of the prognosis process which consists in jointly estimating the model parameters and the current degradation state using the recorded data. This step is realized as long as measurements are available, i.e until the prediction time k p , from which a prediction of the future degradation state is performed.

Future degradation state estimation with the PF

To apply the PF algorithm and to obtain the degradation state pdf, the weight of each particle should be updated at every step. However, these weights depend on the acquisition of new measurements, which are inexistent in the prediction interval. To overcome this difficulty, the state is propagated only using the state model (2.1) while the current particle weights are propagated in time without any changes. In other words, only the prediction step is repeated until the threshold is reached, allowing to obtain the RUL pdf. Such a methodology can be applied as one can consider that the error induced while not changing the particle weights is negligible in comparison to the other sources of error involved in the prognosis process [START_REF] Orchard | A particle-filtering approach for on-line fault diagnosis and failure prognosis[END_REF].

Fatigue crack growth prognosis using an EKF and a PF

In this section, the two model-based prognosis methods presented are applied to estimate the evolution of a nonlinear fatigue crack growth process. The simulation results are shown and compared using performance metrics that are accuracy, precision and timeliness.

Fatigue crack growth model

In order to demonstrate the use of a PF and an EKF for prognosis, a crack evolving in a fuselage panel of an aircraft is considered as the system degradation. When a crack forms in a component, its size and propagation speed must be estimated in order to calculate the RUL of this component. For this purpose, the knowledge of the equation governing the crack growth is first required. A crack growing in the center of a fuselage panel is commonly treated as a through-thickness center crack in an infinite plate. Thus, the well-known Paris' law can be used to assess the crack size as a function of the cycle number [START_REF] Paris | A critical analysis of crack propagation laws[END_REF]:

da dN = C(∆K) m , ∆K = ∆σ √ πa (2.16)
where a is the crack size, N is the number of cycles, ∆K is the range of stress intensity factor and ∆σ is the stress range. C and m are the unknown model parameters to be estimated jointly with the crack length a.

Once a degradation model is available, the goal is now to estimate the state of the degradation and to compute the RUL using the PF and the EKF.

Chapter 2. Filtering methods for prognosis with stochastic assumptions

Simulation settings

Whether for crack growth prognosis with the PF or with the EKF, the following procedure is common to both methods.

Model discretization.

First of all, for simulation purpose, the crack growth model (B.9) is discretized using an Euler direct scheme:

a k = e C k (∆σ √ πa k-1 ) m k ∆N + a k-1 (2.17)
As the parameter C has a very small value, log(C) is considered instead for numerical tractability.

The objective here is to jointly estimate and forecast the crack length a and the two unknown parameters m and log(C). Therefore the augmented state vector is defined as:

X ⊤ k = [ a k log(C k ) m k ] ⊤ , (2.18)
and the evolution equations are:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ a k = a k-1 + e C k (∆σ √ πa k-1 ) m k ∆N log(C k ) = log(C k-1 ) m k = m k-1 (2.19)
Synthetic data generation. This model is used to create the reference crack size data. For this purpose, the values listed in Table 2.1 are chosen. The synthetic crack size measurements are then generated by adding a Gaussian noise with a standard deviation of 0.001. Filter tuning. The model is randomized in order to take the different uncertainty sources into account. Modeling and measurement uncertainty are integrated by adding zero-mean Gaussian noises to the state and measurement equations, that are respectively w ∼ N (0, σ 2 w ) and v ∼ N (0, σ 2 v ). Moreover, the unknown parameters m and C are considered as Gaussian random variables such that m ∼ N (µ m , σ 2 m ) and log(C) ∼ N (µ log(C) , σ 2 log(C) ).

Fatigue crack growth prognosis using an EKF and a PF

Therefore, the following process noise and observation noise covariance matrices, respectively denoted as Q and R, are obtained:

Q = ⎡ ⎢ ⎣ σ 2 w 0 0 0 σ 2 log(C) 0 0 0 σ 2 m ⎤ ⎥ ⎦ , R = σ 2 v (2.20)
The determination of these two covariance matrices is essential and refers to as filter tuning. It is the most challenging step when using a PF or an EKF.

In this work, filter tuning was performed off-line on the basis of a trial and error procedure. On-line filter tuning was not addressed as it is a complex problem and the main purpose here is to explain the prognosis process. In the literature, some references have treated this problem outside the context of prognosis. In [START_REF] Kontoroupi | Online noise identification for joint state and parameter estimation of nonlinear systems[END_REF], an online strategy to jointly estimate, through the formulation of an augmented vector, the state vector, the model parameters and the statistical parameters (covariance and mean) of the process and observation noises is proposed. They assume some conjugate prior distributions for both mean and covariance matrix of noises and perform Bayesian updating using an unscented Kalman filter (UKF). However, they point out the possible limitations of such a method due to the curse of dimensionality.

The variances of the stochastic variables that were used for the implementation of the PF and of the EKF algorithms are listed in Table 2.2.

σ 2 m σ 2 w σ 2 v Value 10 -2 10 -3 (10 -4 ) 2 (10 -3 ) 2
Initialization of the algorithms. The values of the initial state vector have to be determined for the initialization of the filtering algorithms. The initial crack length a 0 is supposed to be known because the prognosis module is launched only if a crack growth is detected in the monitored component. In this example, this value is a 0 = 0.01m. Concerning m and C, as they are material constants, their order of magnitude can be a priori known to initialize the algorithms. In this case, we have chosen m 0 = 3.5 and log(C) = -22.33.

Concerning the parameters that are specific to each filter, they are defined as follows.

Extended Kalman filter

The initial estimation error covariance matrix is chosen here as:

P 0 = ⎡ ⎢ ⎣ 10 -2 0 0 0 10 -2 0 0 0 10 -2 ⎤ ⎥ ⎦ (2.21)
This initial covariance matrix is chosen arbitrarily based on the information available about the initial state estimate. One should take P 0 large enough if the initial estimate is assumed to be not very close. On the contrary, if the initial state estimate is supposed to be close enough from its real value, a smaller P 0 can be defined. Here, the initial value of the state is supposed to be approximately known as the prognosis module is launched when a crack is detected. Therefore, a small value of P 0 was chosen.

Particle filter As for the particle filter, the number of particles should be provided. Here, a set of N part = 1000 particles is considered. The number of particles can be chosen based on a trial and error procedure, according to the user's requirements. Indeed, the higher N part is, the more precision can be gained. However, this is done at the expense of the computational time which is obviously increased. The influence of the number of particles is illustrated in the next subsection in Fig. 2.4, Fig. 2.5 and Fig. 2.6.

Simulation results

To simulate the experiment, crack length measurements were generated every ∆N = 50 cycles, from cycle 0 to cycle 1200. In other words, a crack length measurement is available every 50 cycles and the prediction time N p from which no more data is recorded and the forecasting of the crack length is performed is N p = 1200 cycles. From this time instant, the estimation of the crack size in the future without new measurements was realized until the threshold fixed at 0.0463 is reached (according to [An, Choi, and Kim 2013b]). Using the reference crack size dataset, for N p = 1200 cycles, the true RUL value is 1450 cycles.

Extended Kalman filter.

The measurements and the discretized model 2.4.2 are used in the EKF algorithm to estimate the current crack size in the observation interval from cycle 0 to cycle 1200. Then, the evolution of the crack length in the future is predicted according to the process described in Section 2.2.2. The crack growth prognosis with the EKF is illustrated in Fig. 2.2.

Once the threshold is reached, the failure time cycle is recorded, allowing to compute the RUL distribution. Indeed, the EKF supposes that Gaussian distributions are involved, and gives estimates of the mean and covariance of the RUL. Thus, the RUL distribution that was obtained is a Gaussian distribution N ∼ (µ RU L EKF , σ RU L EKF ) with µ RU L EKF = 950 cycles and σ RU L EKF = 1, 96 × 10 -3 cycles.

As is can be seen in Fig. 2.2, the predicted crack size and therefore the predicted RUL are far from the references. In the same way as for the EKF, the estimation of the current crack length with the measurements is realized from cycle 0 to cycle 1200. Then, the crack growth model is used to forecast the crack size until the threshold is reached.

As it was said before, the output of the EKF at each time step is a Gaussian distribution, and only the mean and covariance of the estimations are propagated. However, the PF propagates each particle through the state transition equation at each time step. This allows to obtain a number of N part crack length trajectories as illustrated in Fig. 2.3 in only 1 simulation run. In this way, the RUL distribution can be obtained by recording the failure time of each particle.

Different number of particles were tested in order to show the influence of N part . In Fig. 2.4, Fig. 2.5 and Fig. 2.6, N part respectively equals 100, 500 and 1000 particles. Only the median trajectories were plotted for the sake of clarity. These figures highlight the fact that the more particles are used, the closer to the reference the estimated crack length is. Thus, N part = 1000 particles was selected for the experiments as it is sufficiently small to keep a reasonable computational time.

The RUL distribution that was obtained with 1000 particles is given in Fig. 2.7. One can see that the computed RUL distribution gives satisfactory results with respect to the reference RUL value. Performance evaluation metrics. The capability of the EKF and of the PF to compute the RUL distribution was investigated in the previous results. However, one can have a more precise performance evaluation using metrics such as accuracy, precision and timeliness, whose expressions were given in Chapter 1. As for the parameters of the metrics, they were chosen here as R 0 = R min = R max = 700 cycles. This corresponds to the admissible prediction error, i.e the acceptable margin for maintenance decisions.

For this purpose, the simulation to compute the RUL distribution with both the EKF and the PF is run 100 times to have a large enough sample for the performance evaluation. The 100 crack length trajectories are shown in Fig. 2.8 for the EKF and in Fig. 2.9 for the PF. Only the median trajectory was plotted for each of the 100 simulation runs of the PF.

A distribution of the mean RUL values of each simulation run was obtained for the EKF (Fig. 2.10) and the PF (Fig. 2.11). Moreover, the numerical values of the performance metrics are depicted in The results show that the EKF based method is less performant than the PF in terms of accuracy and timeliness. The lower performance of the EKF for prognosis has already been reported in [START_REF] Saha | Comparison of prognostic algorithms for estimating remaining useful life of batteries[END_REF]. However, it was not quantified precisely with metrics.

The only disadvantage of the PF is that it requires a more complex implementation and has to propagate the entire distribution at each step, which tends to increase the computational time compared to the EKF algorithm.

Conclusion

Two stochastic methods for model-based prognosis that take modeling and measurement uncertainty have been presented. Both methods are able to generate a pdf associated to the RUL prediction. The two algorithms were compared using a classical fatigue crack growth example and synthetic data. It was shown that the PF method outperforms the EKF one in terms of accuracy and timeliness. This chapter has allowed to explain in more detail all the steps to compute the RUL pdf in a model-based framework. This work also shows the interest of using the metrics for evaluating prognosis performance.

In the next chapter, a more realistic fatigue crack growth prognosis scenario involving a highly nonlinear model is considered. Moreover, in addition to modeling and measurement uncertainty, future loading uncertainty is treated as well.

Chapter 3

Particle filter for fatigue crack growth prognosis in composite materials

Based on the model-based prognosis steps presented in Chapter 1, a more realistic fatigue crack growth problem is considered in this chapter. Fatigue crack growth prognosis in structural elements made of unidirectional fiber-reinforced composites is performed. The considered real test case is an interesting application as fiber reinforced composites are gaining importance for use in structural elements. Indeed, the components last longer because composite materials have a high strength and are corrosion-resistant. Moreover, in addition to modeling and measurement uncertainty, future loading uncertainty is treated as well.

First of all, a physics-based model is required to perform fatigue crack growth prognosis is such materials. Indeed, the Paris' law that was used in Chapter 1 is no longer valid. Thus, the effort was focused on deriving a crack growth model for fiber-reinforced titanium matrix composite materials. Various models are available in the literature, but they induce important computational time and efforts as Fredholm integrals or finite element method are involved. There was therefore a need to adapt these models to perform on-line fatigue crack growth prognosis. Although simplified, the model that was finally obtained is highly nonlinear. The model contains fixed model parameters that depend on material and loading parameters that may vary or not depending on the applied load. Thus, this has allowed to perform crack growth prognosis where model uncertainty and measurement uncertainty are included, but future loading uncertainty is taken into account as well.

Both cases of constant amplitude loading and variable amplitude loading (block loading) are treated in this chapter. Because of its ability to handle uncertainties, high nonlinearities and to perform joint parameter-state estimation (as shown in Chapter 2), a particle filter is used. In a first part, fatigue crack growth prognosis under constant amplitude loading is realized. The loading parameters are constant and known a priori, while the model parameters are jointly estimated along with the crack length. Then, in a second part, fatigue crack growth prognosis under variable amplitude loading is performed. This time, the loading parameters are unknown and change abruptly at unknown time-steps in accordance with the applied variable block loading. A two-sided cumulative sum (CUSUM) algorithm is implemented to detect abrupt load variations and help the particle filter to adapt and learn the new loading parameter values. With the combination of these two techniques, the prognosis module could be informed of the sudden crack length increase, and will correct the predicted remaining useful life. In both case studies, real data from fatigue tests on unidirectional fiber-reinforced titanium matrix composites is used.

The work presented in this chapter has led to the publication of a journal paper [Robinson, Marzat, and Raïssi 2018a] and an international conference paper [Robinson, Marzat, and Raïssi 2018b].

Future loading uncertainty in fatigue crack growth prognosis

The focus is placed on an application of particular interest which is the crack growth prognosis. Indeed, this is a common issue for critical systems as fatigue damage occurs when a structure is subject to repeatedly applied loads. Even if the applied stress values are relatively low, the fact that they are cyclically applied tends to weaken the material and cause its failure. In order to avoid catastrophic events, the crack propagation must be carefully monitored in order to constantly evaluate the remaining useful life of the considered structure. In the literature, the problem of crack growth prognosis was treated for different models and with various methods. A hybrid data-driven and model-based approach based on RVM and model fitting was implemented in [START_REF] Zio | Fatigue crack growth estimation by relevance vector machine[END_REF] and applied to an academic fatigue crack growth process based on the Paris' law. Using the same test case, the crack growth prognosis was solved with a particle-filter based algorithm in [START_REF] Zio | Particle filtering prognostic estimation of the remaining useful life of nonlinear components[END_REF]. [START_REF] Corbetta | Real-time prognosis of crack growth evolution using sequential monte carlo methods and statistical model parameters[END_REF] proposed a particle filter sequentially updated via a Metropolis-Hastings (MH) algorithm for crack growth prognosis on helicopter fuselage panels. A machine learning approach based on artificial neural networks was used to estimate the stress intensity factor (SIF) range, which is required to calculate the crack growth at each cycle. In [START_REF] Zárate | Bayesian model updating and prognosis of fatigue crack growth[END_REF], the SIF range was modeled by a polynomial equation with stochastic coefficients that were computed through Bayesian inference. The future crack length was then predicted using a Markov Chain Monte Carlo algorithm. The ability of these methods to predict the RUL of components subject to fatigue crack growth under measurement and modeling uncertainty has been proved through numerical examples. However, they have assumed known values of current and future 3.1. Future loading uncertainty in fatigue crack growth prognosis loading, which is usually not the case in real life scenario.

In order to address the estimation of the loading amplitude, other researchers have introduced crack growth prognosis methods based on structural health monitoring data. The main idea is to use real-time monitoring data to build models that characterize fatigue loading history, and then based on these models, the future loading can be predicted. In [START_REF] Ling | Integration of structural health monitoring and fatigue damage prognosis[END_REF], both flight parameters data related to acceleration and mass and data recorded from strain gauges were used to estimate and predict loading sequence through an autoregressive integrated moving average (ARIMA) modeling method and a Bayesian approach for the update process. The evaluation of the SIF range required for the crack growth calculation was made with a Gaussian process surrogate model that was previously trained with data from a finite element analysis. In [START_REF] Pais | Predicting fatigue crack growth under variable amplitude loadings with usage monitoring data[END_REF], usage monitoring data from an aircraft (acceleration, airspeed, angle of attack, fuel quantity and Mach number) were converted into a stress time history which was then transformed into a cyclic stress history via a rain-flow counting algorithm. The resulting cyclic stress history was used to consider the effects of variable amplitude loading in the determination of the crack growth direction. However, it was suggested that it could be used as the input into a prognosis method. What these proposed methodologies have in common is that real-time monitoring data related to the loading or flight parameters are required to build a model that characterizes the stress history. Such data may not be available, and if it is the case, the precision of the model depends on the quantity of data. Moreover, a finite element model was needed for the calculation of the SIF range or to train its surrogate model, which highly increases the computational time and efforts.

In this work, the SIF range ∆K is used to characterize the evolution of the crack length in the considered fiber-reinforced composite material. It is assumed that only crack length measurements are available, and no finite element model is used. Indeed, the effort was focused on the derivation of an analytical model of the SIF range. The obtained model contains model parameters that depend on material properties but also loading parameters that depend on external applied loads.

Therefore, fatigue crack growth prognosis under known value of applied load is first considered in order to estimate the model parameters. For this purpose, a joint parameter-state estimation is performed with a particle filter because of its ability to handle uncertainties and the high nonlinearities of the SIF range model. The use of an extended Kalman filter (EKF) to solve the same fatigue crack growth prognosis problem was also investigated as it is easy to implement and requires lower computational time [START_REF] Wang | Fatigue crack prognosis using Bayesian probabilistic modelling[END_REF]. The results were compared to those obtained with the particle filter using prognosis metrics.

In a second study, the problem of fatigue crack growth prognosis under unknown constant and then unknown variable amplitude loading is addressed. Variable amplitude loading can be divided into three categories [Laseure, Schepens, Micone,

Chapter 3. Particle filter for fatigue crack growth prognosis in composite materials

and De Waele 2015]: random loading, block loading and simple loading. In this work, block loading is considered, which means that loading parameters vary stepwise in time. Thus, using the model parameters estimated in the first study, the joint estimation of the crack length and of the unknown loading parameters was performed with a particle filter. However, if a significant abrupt load variation happens, the particle filter is not able to track the abrupt change in loading parameters values.

To solve this problem, a detection algorithm, namely the two-sided cumulative sum (CUSUM), is implemented in parallel with the particle filter to catch any sudden deviation of loading parameters. This way, the particle filter will receive the information and then reinitialize the previous estimates in order to learn the new loading parameter values. This detection step enables to adapt the prognosis algorithm to the changes in crack length evolution and therefore give a more accurate prediction of the RUL. Indeed, an increase in the amplitude loading can significantly accelerate the crack growth rate and lead to the critical crack length causing the failure of the component.

In both studies, real crack length trajectories provided by ONERA/DMAS (Département Matériaux et Structures/Materials and Structures Department) were used, and the prognosis metrics that were used to validate the ability of the particle filter to perform fatigue crack growth prognosis in such materials are accuracy, precision and timeliness (previously defined in Chapter 1).

Problem statement

The problem of fatigue crack growth prognosis in fiber-reinforced titanium matrix composite materials is addressed. This will allow to continuously assess the RUL of a structure made of such composite material in which a fatigue crack growth has occurred. As a model-based prognosis approach is used, the first requirement is to establish an analytical model to assess dynamic crack propagation in such materials. Before introducing this model, the experimental procedure to collect the real crack growth data used in this work is described.

Material and experimental procedures

A composite material is the combination of two or more different materials in order to create a superior material with different properties (stronger, lighter, etc.). Composites are mainly made up of two constituent materials: matrix and reinforcement. There are three main kinds of materials that are used for the matrix (polymer, metal and ceramic) and also different forms of reinforcement material (particles, fibers or laminates). In this work, the proposed methodology is applied to fatigue crack growth prognosis in titanium metal matrix composites with silicon carbide fibers used as reinforcement materials. The fatigue test data were provided by ONERA/DMAS, and were previously used in [START_REF] Maire | Multi-scale model of crack propagation in unidirectional metallic matrix composites[END_REF] to establish and validate a model to describe the fatigue crack growth in the concerned specimens, but the problem of RUL prognosis was not addressed. For consistency, the experimental protocol is summarized in what follows.

The titanium-matrix composite studied was SCS-6/Ti-6242. Ti-6242 is a near alpha titanium alloy with the composition Ti-6Al-2Sn-4ZR-2Mo (percent by weight), and this matrix was reinforced with 140 µm diameter SCS-6 Textron fibers. The fibers are regularly spaced in the matrix in such a way as to obtain a unidirectional composite. A parallelepiped notched specimen (Fig. 3.1) with a nominal width of 8 mm, thickness of 2.1 mm and length of 160 mm was machined from the composite material in a manner that the length of the specimen is parallel to fiber axis. An elliptic notch was drilled in the middle of the specimen in order to initiate the crack growth. During the experiment performed at ONERA/DMAS, crack growth trajectories of 7 samples of different notch lengths were recorded. Some specimens were tested under constant amplitude loading, while the others were subjected to load variations. The Table 3.1 gathers the information about the fatigue tests. For every specimen, a uniaxial cyclic loading oriented along the fiber direction was applied. The stress ratio R = σ min /σ max was equal to 0.1. The fatigue tests were performed at room temperature and at 400 • C under a frequency of 50 Hz. In order to measure the crack length, photographs of the crack extension were recorded by a digital camera monitored by a computer. The reader may refer to [START_REF] Thurner | Real-time detection and measurement of cracks in fatigue test applications[END_REF]] to be more informed about the real-time detection and measurement of cracks. As shown in Fig. 3.2, the typical cracking geometry involved the propagation of a crack in the matrix on each side of the notch, propagating perpendicularly to the loading direction.

The crack propagated through all the thickness of the specimens and across their entire width without causing the rupture of the composite material. This phenomena is due to unbroken fibers that have bridged the matrix crack. Indeed, as the fiber stress level did not exceed the value of the fiber strength, no fiber broke during the experiments. The constraints of these bridging fibers can be modeled by the distribution of a closure pressure P acting in the direction opposite to the applied 53 652, stress σ a in the bridged zone (Fig. 3.3). In this case, the bridged zone is equivalent to the crack length minus the notch length 2a 0 .

Chapter 3. Particle filter for fatigue crack growth prognosis in composite materials

Description of the crack propagation model

It is assumed that the composite has a linear elastic behavior except in a very small region at the crack tip, therefore the principles of linear elastic fracture mechanics can be applied. The following modified Paris' law that was previously used in [START_REF] Maire | Multi-scale model of crack propagation in unidirectional metallic matrix composites[END_REF] to study the same fatigue test dataset can then be used to model the Mode I crack propagation:

da dN = C(∆K m -K th ) m (3.1)
where a is the crack size, N is the number of cycles. Since the crack growth observed in the composite material was limited to the matrix cracking, the effective crackdriving force is assumed to be the SIF range ∆K m experienced by the matrix. The constants C, m and K th depend on matrix properties. K th is the threshold value of the SIF below which no crack growth occurs.

The calculation of the SIF range ∆K m is the main challenge in this modeling stage. A simple expression of the SIF range can be available. However, in more complex structures as composite materials, an analytical closed form of the SIF does not always exist or is too complicated to establish. In these cases, finite element simulations can be run to calculate the SIF values associated to different crack lengths. Then, based on the obtained database of SIFs, a regression model is used to allow the evaluation of the SIF for any crack length. This technique was used in [START_REF] Corbetta | Real-time prognosis of crack growth evolution using sequential monte carlo methods and statistical model parameters[END_REF] where a machine learning approach based on artificial neural networks was used to provide estimates of SIFs. In [START_REF] Neerukatti | Fatigue life prediction using hybrid prognosis for structural health monitoring[END_REF], two regression techniques were used, namely least absolute shrinkage and selection operator (LASSO) and relevance vector machine. Although very efficient to obtain good estimations of the SIF, the main drawback of such an approach is that the learning phase can be time consuming and a significant amount of data is required. In this work, an analytical expression of the stress intensity factor for bridged cracks in composite materials was established, based on different studies found in [START_REF] Johnson | Life prediction methodology for titanium matrix composites[END_REF].

The first step to determine ∆K m is to relate it to the continuum SIF range ∆K tip which is the homogenized composite stress intensity factor. Many discrete-continuum relationships were proposed, and three of them were compared in [START_REF] Bakuckas | Application of fiber bridging models to fatigue crack growth in unidirectional titanium matrix composites[END_REF]]. In their study, the one that gave the best results in the modeling of fiber-bridging effect on ∆K m was established by [START_REF] Mcmeeking | Matrix fatigue cracking in fiber composites[END_REF]:

∆K m = ∆K tip (3.2)
In order to calculate ∆K tip , the weight function technique proposed by [START_REF] Bueckner | Novel principle for the computation of stress intensity factors[END_REF]] that allows to calculate a stress intensity factor for arbitrary stress distributions is applied. In this case study, the applied stress range on the crack surface is

Chapter 3. Particle filter for fatigue crack growth prognosis in composite materials

∆σ br = ∆σ a -∆P (x), therefore we have:

∆K tip = 2 ∫ a 0 (∆σ a -∆P (x))G(x, a, c)dx (3.3)
where G(x, a, w) is a weight function that depends on the geometry of the specimen, c is the specimen width and 2a is the crack length.

For a center crack configuration in a finite width specimen, the following expression of ∆K tip is finally obtained [START_REF] Zheng | High temperature/high frequency fatigue crack growth in titanium metal matrix composites[END_REF]:

∆K tip = F ∆σ a √ πa -2F √ a π ∫ a a 0 ∆P (x) √ a 2 -x 2 dx (3.4)
where F is a geometric factor that depends on the specimen width (various expressions can be found in [START_REF] Tada | The stress analysis of cracks[END_REF]), 2a 0 is the length of the unbridged zone (i.e. notch length) and x is the distance from the crack center. It can be noted that this expression of ∆K tip takes into account both the contribution of the applied remote stress ∆σ a and the contribution of the bridging stress ∆σ br produced by unbroken fibers in the crack wake.

Finally, the model that describes crack growth propagation in the studied specimens is given by:

da dN = C ( F ∆σ a √ πa -2F √ a π ∫ a a 0 ∆P (x) √ a 2 -x 2 dx -K th ) m (3.5)
The determination of the change in closure pressure ∆P (x) is the critical issue of the fatigue crack propagation problem in fiber reinforced composites. In the literature, two types of analytical models have been widely used to calculate ∆P (x): (i) the shear-lag model and (ii) the fiber pressure model. For more details about these two approaches, the reader may refer to [START_REF] Ghosn | Modeling of crack bridging in a unidirectional metal matrix composite[END_REF]. As an alternative to these methods, [START_REF] Davidson | The micromechanics of fatigue crack growth at 25 C in Ti-6al-4v reinforced with SCS-6 fibers[END_REF]] have used a uniform closure pressure over the entire bridged zone and obtained satisfying results for SCS-6/Ti-6Al-4V composite materials. This approach was considered in this paper as it allows to avoid the integration of ∆P (x) to calculate ∆K tip , which highly decreases the computational time. Indeed, this is an important criteria for on-line RUL prognosis.

Adopting this constant approximation of the closure pressure leads to the following model:

da dN = C ( F ∆σ a √ πa -2F √ a π ∆P ∫ a a 0 1 √ a 2 -x 2 dx -K th ) m (3.6)
where

∫ a a 0 1 √ a 2 -x 2 dx = π 2 -arcsin ( a 0 a ) .
(3.7)

Fatigue crack growth prognosis with a particle filter for joint parameter-state estimation

Finally, the model that is used for crack growth prognosis in the studied fiberreinforced titanium matrix composite material is:

da dN = C ( F ∆σ a √ πa -2F √ a π ∆P ( π 2 -arcsin ( a 0 a )) -K th ) m (3.8)
The model (3.8) that was obtained involves two different types of parameters. One needs to distinguish between model parameters and loading parameters:

• m and C are constant model parameters which depend on the material under study and various testing conditions (stress ratio R, temperature, etc.). They can be estimated for a given specimen and considered to have the same values for the other similar test specimens;

• ∆σ a and ∆P are loading parameters directly related to the external loading applied to each specimen. Their values may therefore vary according to the loading conditions. The estimation of these two loading parameters is an interesting feature as most of the existing prognosis methods have assumed known values of current and future loading conditions, which is not the case in real life scenarios. Their real-time monitoring would allow the prognosis system to catch any abrupt load variation and update the RUL estimation. Indeed, depending on the load variation suffered by the structure, the crack growth rate may rapidly increase and cause its sudden failure.

Because of their different nature (constant or time-varying), these two sets of parameters will not be estimated at the same time for the fatigue crack prognosis in the studied composite material. Firstly, a dataset of crack growth trajectory from a specimen tested under a known constant amplitude loading will be used to estimate the model parameters m and C and perform fatigue crack growth prognosis. Then, using these estimated values of m and C, the case of fatigue crack growth prognosis under variable amplitude loading where ∆σ a and ∆P are unknown will be considered.

In both studies, the parameters are jointly estimated along with the crack length with a particle filter. However, in the second case study, ∆σ a and ∆P are time-varying parameters that change abruptly with the applied load. A CUSUM algorithm is then used to detect load variations and help the particle filter to adapt and learn the new loading parameters values.

Fatigue crack growth prognosis with a particle filter for joint parameter-state estimation

A particle filter was chosen to perform the joint estimation of the unknown parameters along with the crack length in order to calculate the RUL of the composite 57

Chapter 3. Particle filter for fatigue crack growth prognosis in composite materials

specimens subjected to fatigue crack growth. The ability of the particle filter to perform uncertainty propagation but also to deal with the high nonlinearities of the fatigue crack growth analytical model (3.8) has motivated this choice. Indeed, the performances of the particle filter were investigated in Chapter 2. Moreover, as the particle filter is a Bayesian filtering method, the RUL predictions are updated and improved as new crack length data is collected.

When specimens were tested under variable amplitude loading (block loading), loading parameters are piecewise constant functions of time. In order to detect load variations and to keep an accurate estimate of these loading parameters, a two-sided CUSUM algorithm is added to the particle filter so it can keep on tracking the time-varying parameters. The idea is to jointly estimate the unknown loading parameters along with the crack length using the particle filter, while monitoring their estimations with the two-sided CUSUM algorithm to detect abrupt load variations. If a load variation happens, the CUSUM algorithm informs the particle filter which has then to reinitialize its previous estimates of the loading parameters and learn their new values.

The principle of the joint parameter-state estimation using the particle filter was exposed in Chapter 2. In this section, the CUSUM algorithm is given, then, it will be explained how the particle filter and a detection algorithm can be combined to correctly estimate the crack length in a structure subjected to variable amplitude loading.

Two-sided CUSUM algorithm

A two-sided CUSUM algorithm is used for the detection of load variation. It is the combination of two CUSUM algorithms, where one is for the detection of an increase in the mean of the monitored variable, and the other one to detect a decrease in the mean. The general idea is to calculate a cumulative sum S k that depends on the monitored process ∆x, on its initial mean value µ 0 and on the minimal size of change to detect denoted by ν. And when the value of the sum exceeds a defined threshold value S th , a change in the mean value is detected. Therefore, the two-sided CUSUM algorithm is based on the following equations:

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ S + k = max ( 0, S + k-1 + ∆x k -µ 0 -ν 2 ) S - k = max ( 0, S - k-1 -∆x k -µ 0 -ν 2 ) N detect = min{k : S + k ≥ S th ∪ S - k ≥ S th } (3.9)
where N detect is the time at which the detection is made. There are two parameters that have to be chosen in this algorithm: S th and ν. This choice depends on how the signal to process looks like, and for a Gaussian distribution, one can set S th = 2 σ ∆x ν where σ ∆x is the standard deviation of ∆x. Further information about the two-sided CUSUM algorithm can be found in [START_REF] Blanke | Diagnosis and fault-tolerant control[END_REF].

Prognosis methodology with the particle filter and the detection algorithm

The general methodology that was presented above is adopted when the constant model parameters m and C are jointly estimated along with the crack length. In the case where the crack is propagating under unknown variable amplitude loading (block loading), the loading parameters are unknown and time-varying (piecewise constant functions of time). However, the SIR particle filter itself cannot handle sudden change in loading parameters. Indeed, between two load variations, the loading parameters have constant values. The marginal posterior parameter distributions become increasingly concentrated around the true parameters values as more data is available. Thus, if some abrupt load variation occurs, the new loading parameters values may fall far into the tails of the parameter distributions, leading to particles with negligible weights. The particle filter is not able to estimate the new parameters values and a reinitialization of the particle system is then needed. The two-sided CUSUM algorithm is thus implemented to detect load variations and trigger the re-initialization of the particle filter for the learning of the new loading parameters values.

The constant monitoring of loading parameters with the detection algorithm allows not only to detect abrupt load variations but also to reinitialize the parameters values right after the alarm. This would help the particle filter to converge more quickly and more easily to the actual crack length. This procedure leads to the introduction of a new uncertain parameter which is the initial loading parameters values after load variation. In most cases this value is unknown, but the user may have an order of magnitude of it. Indeed, depending on the monitored system, the critical input value that can accelerate the degradation state evolution can be obtained from expert knowledge. Therefore, because of the uncertainty associated to this value, it is included in an interval L = [L, L].

Numerical applications

The aim of this section is to show how on-line fatigue crack growth prognosis with a SIR particle filter can be assessed in fiber-reinforced titanium matrix composite materials. The case of a fatigue crack prognosis with unknown model parameters m and C but known loading parameters ∆σ a and ∆P is first presented. For this purpose, dataset #4 of Table 3.1 where the specimen was tested under constant

Chapter 3. Particle filter for fatigue crack growth prognosis in composite materials

amplitude loading is used. In order to justify the choice of a particle filter rather than an extended Kalman filter (EKF) which is easier to implement and is less time consuming, the performance of these two stochastic filters will be compared. Then, in a second part, a more complex and realistic case of fatigue crack growth prognosis where the loading parameters ∆σ a and ∆P are unknown is considered. In this second example, dataset #3 of Table 3.1 from a specimen tested under variable amplitude loading is used.

In order to evaluate the performance of the particle filter and of the EKF in estimating the RUL of the considered composite materials, the three prognosis metrics (accuracy, precision and timeliness) presented in Chapter 1 are calculated. Accuracy measures the degree of closeness of the predicted RUL to the actual RUL, and its values are between 0 and 1 where 1 gives the best accuracy. Precision evaluates the narrowness of the interval in which the RUL predictions fall, and ranges between 0 and 1 which reflects the highest precision. Finally, timeliness indicates the relative position of the predicted RUL pdf along the time axis with respect to the occurrence of the actual failure event. There are three cases: (i) the failure occurs after the predicted failure time, (ii) the failure occurs at the same time as the predicted failure time, and finally, (iii) the failure occurs earlier than predicted. This last case must be absolutely avoided, that is why the timeliness function allows to penalize late predictions. Timeliness has positive values and 0 is the best score.

For both case studies, the input parameters required for the implementation of the algorithms are provided, then the simulation results are shown. Finally, the prognosis methods are compared using the performance metrics that have been defined above. The following discrete-time form of the crack growth model (3.8) is used:

a k+1 = a k + C∆N [ F ∆σ a k √ πa k -2F √ a π ∆P k ( π 2 -arcsin ( a 0 a k )) -K th ] m
(3.10)

Unknown model parameters m and C

In this first case study, a crack growth trajectory dataset where the loading amplitude is maintained constant is used (dataset #4). The loading parameters ∆σ a and ∆P are then constant and equal to 196, 2 MPa. The objective here is to jointly estimate and forecast the crack length and the two unknown material parameters m and C, therefore the augmented state vector is defined as

X ⊤ k = [ a k C k m k ] .
Filter tuning. The variances of the stochastic variables that were used for the implementation of the PF and the EKF algorithms are listed in Table 3.2. Filter tuning was performed off-line on the basis of a trial and error procedure, as it was described in Chapter 2. The process noise w is assumed to be additive and the measurements noise v is unknown. (10 -2 ) 2 (10 -3 ) 2 (10 -3 ) 2 (10 -3 ) 2 Initialization of the algorithms. The values of the initial state vector have to be determined for the initialization of the filtering algorithms. The initial crack length a 0 is supposed to be known because the prognosis module is launched only if a crack growth is detected in the monitored component. In the tested specimen, this value is a 0 = 3.1mm. Concerning m and C, as they are material constants, their order of magnitude can be a priori known to initialize the algorithms. In this case, we have chosen m 0 = 2.4 and log(C) = -23.65. Finally, crack length measurements were collected every ∆N = 1172 cycles and a total number of 500 particles was used by the particle filter. This number of particles was chosen after a trial and error procedure, it was the best compromise between computational time and precision of the estimation. In the figures presented in the next subsections, the plotted curves correspond to the median trajectory.

The crack growth prognosis has been done for different prediction time cycles N p , which means that the number of available measurements in the observation interval has been increased over cycles, as in on-line prognosis. From prediction time cycle N p , the forecasting of the future crack length was performed without any measurements until the threshold was reached in order to deduce the RUL value. This critical crack length was fixed at 7.5 mm.

Performance evalutation metrics.

As it was done in the case study of Chapter 2, 100 simulation runs of the filtering algorithms were performed in order to evaluate the performance of the algorithms in terms of accuracy, precision and timeliness. Concerning the choice of the performance metrics parameters [I min ,I max ] and I 0 defined in Chapter 1, they depend on the material and user requirements. They can be derived from prior knowledge about the fatigue behavior of the material under study characterized by its S -N curve [START_REF] Kawai | A failure-mode based anisomorphic constant life diagram for a unidirectional carbon/epoxy laminate under off-axis fatigue loading at room temperature[END_REF]. Based on these considerations, the user can define an allowable margin as a decision criterion. In this case, we have chosen I 0 = 2 × 10 5 cycles. As for the interval [I min ,I max ], its width equals I 0 = 2 × 10 5 cycles.

The evolution of the crack growth prognosis for different values of N p are shown from Fig. 3.4 to Fig. 3.6 for the EKF and from Fig. 3.7 to Fig. 3.9 for the PF. The performance metrics values for both methods are listed in Table 3.3. The results show that the EKF did not manage to predict correctly the future crack length until time cycle N p = 2.794 × 10 5 cycles. Indeed, the EKF gives reliable results only if the system is locally linear and is not highly nonlinear, or when there is no joint

Chapter 3. Particle filter for fatigue crack growth prognosis in composite materials

parameter-state estimation as in [START_REF] Wang | Fatigue crack prognosis using Bayesian probabilistic modelling[END_REF]. The only way to obtain satisfactory results with the EKF was to decrease the standard deviation of the exponent parameter m to σ m = 10 -6 , however it is not a realistic value because the parameter m generally varies by 10% from one specimen to another. As regards the particle filter, even at only N p = 1.198 × 10 5 cycles, the RUL prognosis gives appropriate results in terms of accuracy, precision and timeliness, and the metric values improve as more measurements are available. The high performance of the particle filter for the fiber-bridged crack prognosis was highlighted in the results and confirms the results obtained in Chapter 2 on an academic test case with a Paris' law. Moreover, the choice of the constant value of the closure pressure is appropriate in this case. Thus, the proposed model for fatigue crack growth propagation in fiber-reinforced composites can be validated.

Unknown loading parameters ∆σ a and ∆P

This second case study is based on dataset #3 where the cracked specimen was subjected to variable loading amplitude. A joint input-state estimation will allow the constant monitoring of external loads while estimating the crack length. The augmented state vector is then defined as

X ⊤ k = [ a k ∆σ a k ∆P k ] .
Filter tuning. The variances of the stochastic variables that were used for the implementation of the particle filter algorithm are listed in Table 3.4. As in the previous subsection, filter tuning is based on a trial and error procedure. Moreover, the process noise w is supposed to be additive and the measurements noise v is unknown.

Table 3.4: Distributions of random parameters for the estimation of ∆σ a and ∆P Parameter ∆σ a ∆P w v Variance (0.07) 2 (0.07) 2 (10 -3 ) 2 (5 × 10 -3 ) 2

Initialization of the algorithm. The initial crack length a 0 is supposed to be known because the prognosis module is launched only if a crack growth is detected in the monitored component. In the tested specimen, this value is a 0 = 0.7mm. Concerning the loading parameters ∆σ a and ∆P , their initial values are assumed to be included in the interval L = [300, 400] MPa.

In this part, the material parameters m and C are assumed to be constant variables as the focus is placed on the estimation of the unknown loading parameters. They have already been calculated in the previous case study using dataset #4, and we have m = 2.4 and C = 5.4 × 10 -11 . Finally, a total number of 500 particles was used and crack length measurements were collected every ∆N = 1172 cycles.

Fatigue crack growth under unknown constant amplitude loading

The first step of prognosis consists in jointly estimating the current crack length a and the unknown loading parameters ∆σ a and ∆P using the crack growth model and the collected crack length data at each time step k corresponding to a cycle number N . Then, from a prediction time cycle N p , the forecasting of the future state vector is performed without any new measurements.

First of all, in order to demonstrate the robustness of the proposed particle filter to the initial unknown loading parameters values, the prognosis results for different initial values of ∆σ a and ∆P in L 0 have been plotted in Figure 3 Moreover, one can notice that the crack length does not increase and the estimated unknown loading parameters remain constant after the transient state (Figure 3.11). Indeed, the crack growth stops after some extension because of the fiber bridging phenomenon mentioned in Section 3.2. Therefore, the crack growth reaches a steadystate and the failure threshold will not be reached.

In the following, a case of crack growth under variable amplitude loading is treated to study the ability of the proposed methodology to jointly estimate the crack length and the unknown loading parameters even in these circumstances.

Fatigue crack growth under unknown variable amplitude loading

In this subsection, crack growth prognosis under unknown variable amplitude loading is considered. The crack length measurements are from the real dataset obtained during fatigue test on fiber-reinforced titanium matrix composite materials presented in Section 3.2.

As a first step, the estimation of the crack length and the unknown loading parameters is performed using the available measurements. Then, the forecasting step is realized from the prediction time cycle N p . Before load variation, the prognosis results are the same as described in the previous subsection. After load variation, the forecasting of the future crack growth without new measurements is made at different prediction time cycles N p . The evolution of the future crack length predictions are depicted in Figure 3.12 and Figure 3.13. The first part of the crack length evolution is not shown as it is the same as in Figure 3.10.

The higher amplitude of the external applied load has led to an increase in the crack length, which indicates that the critical stress has been attained. The results show that the particle filter has some difficulties to converge even after several time cycles of estimation using the measurements. Indeed, the particle filter needs more data in order to estimate the unknown loading parameters whose values have significantly increased after this abrupt load variation. This problem might be addressed by increasing the number of particles. In this study, 500 particles were used, and a test with 3000 particles was realized. The results were almost the same, and this number cannot be further increased because the computational time would be too important, which is not suitable for online applications.

In order to circumvent this issue, the two-sided CUSUM algorithm is integrated to the particle filter. The unknown loading parameters ∆σ a and ∆P are constantly monitored to detect the sudden load variation. Indeed, as these parameters are related to the applied load, when the loading amplitude changes, their values change as well (Figure 3.14). It can be seen that this monitoring must start after the transient state to avoid any false detection.

The parameter µ 0 of the CUSUM algorithm is the mean value of the estimated parameters in the transient state. The minimal size of change to detect in the unknown loading parameters variables was fixed to ν = 2 MPa. This value was chosen based on the parameter estimations shown in Fig. 3.11. The two-sided CUSUM algorithm has detected the variation at N detect = 10.881 × 10 6 cycles while the actual load variation time is N load = 10.868 × 10 6 cycles.

After the detection of load variation, the values of ∆σ a and ∆P are reinitialized. The choice of the interval L is based on the a priori knowledge of the necessary load amplitude that may cause a rapid propagation of the crack in the considered specimen. The values ∆σ a crit and ∆P crit associated to the critical load are known to be around 450 MPa. Therefore, in order to take the uncertainty associated to this [400,500]. This has allowed to estimate and then forecast an interval that contains the estimated crack length after load variation. The bounds of this interval are derived from two extreme loading cases. The evolution of this interval for different prediction cycles N p is shown in Figure 3.15 and in Figure 3.16.

First of all, the gain of the detection algorithm and the reinitialization of the loading parameters values is highlighted in these figures. It can be seen that even if the loading parameter values after the load variation are uncertain, the interval that was derived from the interval L still gives an accurate prediction of the future crack length. Moreover, the predictions improve as more data is available.

Performance evaluation metrics. Even if the plotted results

give an idea about the efficiency of the proposed methodology, a more precise performance evaluation using metrics such as accuracy, precision and timeliness is performed. The metrics were calculated for several reinitialization values of ∆σ a and ∆P included in L. Indeed, it allows to take the uncertainty related to these values into account. Thus, for each value of ∆σ a and ∆P included in L, 100 RUL pdf computations were simulated to calculate the metrics. Then, the mean values and the standard deviation of each metrics in the predicted interval are given in Table 3.5. The calculated metrics confirm that the proposed methodology gives satisfactory results in terms of accuracy, precision and timeliness despite the high uncertainty on the reinitialized values of the unknown loading parameters after load variation. Moreover, the obtained values show the usefulness of the timeliness metric. Indeed, we can see that the highest value of timeliness is at the prediction cycle N p = 12.337×10 6 . This is explained by the fact that a significant part of the predicted RUL have fallen after the actual RUL value. Thus, even if the accuracy is the highest one, the timeliness value must be taken into account carefully to optimize maintenance decisions and avoid catastrophic events. As for the case of the particle filter without the detection algorithm, the metrics could not even be calculated. Indeed, the particle filter did not manage to estimate the increase in crack length after the load variation. Therefore, the crack growth was so slow that the critical crack length could never be reached, making the calculation of the metrics impossible.

Chapter 3. Particle filter for fatigue crack growth prognosis in composite materials

1 1.1 1.2 1.3 1.4 1.5 x 10

Conclusions

In this chapter, fatigue crack growth prognosis in unidirectional fiber-reinforced composites using real data from fatigue tests on titanium matrix composites has been considered. First of all, the model that describes fatigue crack growth in such materials has been derived. For this purpose, the model that was previously developed in [START_REF] Maire | Multi-scale model of crack propagation in unidirectional metallic matrix composites[END_REF] has been simplified. Indeed, the problem in model-based prognosis is to make a compromise between a very accurate model that induces very high computational time, and a simplified model that will reduce computational time and effort. This simplified model was obtained by making assumptions on the physics-based model, and the resulting modeling uncertainties were represented with uncertain parameters and additive noise.

Then, using this nonlinear model and real crack growth data, two cases of fatigue crack prognosis in fiber-reinforced composite materials were treated. In a first study, the joint estimation of model parameters that depend on material in conjunction with the crack length was performed while assuming known value of the applied load. The performances of a particle filter and an extended Kalman filter (EKF) for the crack growth prognosis were compared using three prognosis metrics (accuracy, precision

Chapter 3. Particle filter for fatigue crack growth prognosis in composite materials

and timeliness). The results highlighted the high performance of the particle filter and the difficulty of the EKF to predict the future crack growth in presence of high nonlinearities. It was also shown that the more measurements are collected, the more the performance metrics of both filters were improved. However, the particle filter needed much less measurements to obtain a predicted RUL value close to the real one, which also proves the efficiency of the particle filter over the EKF in the studied crack growth prognosis.

In a second case study, a real crack growth dataset recorded under variable amplitude loading was used to perform fatigue crack growth prognosis in unidirectional fiber-reinforced composites under unknown applied load. A model-based on-line prognosis method that is able to estimate and forecast unknown loading parameters was proposed. A particle filter was used for the joint parameter-state estimation and a two-sided CUSUM algorithm was integrated to detect load variations. Indeed, it was noted that after an abrupt load variation, the particle filter had some difficulties to converge to the degradation state whose trajectory had suddenly changed. Therefore, the role of the CUSUM algorithm is to monitor the unknown loading parameters values, and to give an alert when a load variation is detected. Once the particle filter has received this alert, the estimated values of the unknown loading parameters are reinitialized in an interval that is chosen from a priori knowledge. The association of these two algorithms have enabled to keep the RUL predictions accurate even after load variation.

Finally, the simplified crack growth model that was proposed for fatigue crack growth prognosis in unidirectional fiber-reinforced composite materials can be validated. Indeed, even when considering that the absolute value of the closure pressure is equal to the applied load, the crack growth prognosis results were satisfactory. That is an important aspect in on-line prognosis because the computational time must be minimized.

Chapter 4

Set-membership framework for model-based prognosis

In this chapter, two set-membership state estimation techniques based on interval analysis are applied to the problem of RUL prognosis. In this context, uncertainty propagation is performed without no stochastic assumption and only the bounds of the state and measurement noises are required. This allows the estimation of the degradation state in a guaranteed way.

In the first approach, a set-membership methodology using constraint satisfaction is adopted to estimate and then forecast the degradation state. For this purpose, the feasible domain of the model parameters is first estimated based on the available measurements.Then, using the natural inclusion function of the degradation model and the projection of the outer approximation of the feasible domain of the parameters, the future degradation behavior is predicted and the interval that contains the RUL of the system is finally computed. The proposed methodology is illustrated using the same nonlinear crack growth example as in Chapter 2.

The second method belongs to the category of set-membership observers. The case where the degradation cannot be directly measured (inaccessible or difficult to measure) is considered. However, it can be detected through the drift of a model parameter. This parameter drift is viewed as an unknown input, whose reconstruction allows the estimation of the degradation signal by solving an inverse problem. Thus, in order to estimate and forecast the bounds of the degradation state, the bounds of the degradation signal should be reconstructed first. The RUL calculation is composed of three main steps. First, the bounds of the state and the unknown input are estimated, while taking the model and measurement noises into account. For this purpose, a novel methodology to synthesize an interval unknown input observer (UIO) for linear time-invariant (LTI) discrete-time systems is proposed. The bounds of the degradation signal are then deduced by solving an inverse problem based on interval arithmetics. The second step consists in estimating the feasible domain of the parameters of the degradation model based on a set-membership methodology using constraint satisfaction. Finally, as a third step, the natural inclusion function of the degradation model is used to compute the bounds of the future degradation state until failure threshold is reached, allowing the deduction of RUL interval. In order to demonstrate the efficiency of the proposed model-based prognosis methodology, the degradation of a suspension system is studied. A crack that cannot be directly measured is evolving in the suspension spring, and is affecting its stiffness constant, which is a parameter of the suspension dynamical model. This work has led to the publication of an international conference paper [Robinson, Marzat, and Raïssi 2017b].

Set-membership methodology using constraint satisfaction

In this section, a set-membership methodology is proposed based on a guaranteed estimation and prediction of models described by (2.1), (2.2) where only the bounds of the noises and disturbances are available without any additional stochastic assumption. Disturbances w and noises v satisfy |v k |≤ V , |w k |≤ W for some positive bounds V and W . Since w and v belong to intervals, the parameter vector θ and the degradation state x cannot take single values but they also belong to some compact domains. The proposed set-membership methodology to compute the RUL is based on two steps:

1. The available measurements over the time step interval [k 0 , k p ] are used to estimate the feasible domain of the parameter vector θ given by

Θ = ⎧ ⎪ ⎨ ⎪ ⎩ θ ∈ R q | x k = f (x k-1 , θ, w k ), h(x k , v k ) ∈ [y m k -V, y m k + V ], ∀k ∈ [k 0 , k p ], ∀x k-1 ∈ [x k-1 ], ∀|w k |≤ W ⎫ ⎪ ⎬ ⎪ ⎭ (4.1)
where y m k are the measurements and [y m k -V, y m k + V ] is the domain of the output taking the noises into account.

2. The estimated feasible parameter domain is used to predict the degradation behavior for time steps k > k p in order to estimate the remaining useful life of the system.

The estimation and prediction steps are based on interval tools to deal with the bounded uncertainties.

Interval techniques

Interval analysis techniques represent a powerful tool to tackle uncertainty propagation without no stochastic assumption. Indeed, the evaluation of the whole set of possible model outputs could be performed using only one interval evaluation. A real interval [a] = [a, a] is a connected and closed subset of R. The set of all real intervals of R is denoted by IR. Real arithmetic operations are extended to intervals (see [START_REF] Moore | Interval analysis[END_REF]). Consider an operation • ∈ {+; -; * ; /} and [a], [b] two intervals, then:

[a]

• [b] = {x • y | x ∈ [a], y ∈ [b]}
The width of an interval [a] is defined by w[a] = aa and its midpoint by mid[a] = (a + a)/2. The midpoint represents a point estimation of a variable and the radius is the uncertainty.

Let f : R n → R m be a vector function and [x] an interval vector (also called a box) of R n . An inclusion function, denoted by [f ], aims at computing an interval [f ] ([x]) that contains the images of f over the box [x], and can be defined as:

[f ]([x]) ⊂ {f (x) | x ∈ [x]}, (4.2)
or equivalently, an interval function [f ] : IR n → IR m is an inclusion function for f if:

∀[x] ∈ IR n , f ([x]) ⊆ [f ]([x]) (4.3)
An inclusion function of f could 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 constraint satisfaction problem (CSP) is defined by a set of n variables x = x 1 , x 2 , . . . , x n and a set of m constraints C 1 , C 2 , . . . , C m . Each variable x i has an initial nonempty domain D i of possible values. Each constraint C i involves a subset of the variables and specifies the possible combinations of values for such subset. A state of the problem is defined by an assignment of values to some or all of the variables, x i = v i , . . . , x j = v j . An assignment that does not violate any constraint is called a consistent assignment. A complete assignment is one in which every variable is mentioned, and a solution to a CSP is a complete assignment that satisfies all the constraints. In contrast to conventional techniques, interval methods do not suffer from local convergence and the computed set is guaranteed to contain the global solutions. In addition, an empty set is returned if the CSP has no solution in the initial searching domain. The goal of propagation techniques is to reduce as much as possible the domains for the variables without losing any solution. The most popular approach is based on the Waltz filtering algorithm [START_REF] Waltz | Generating semantic descriptions from drawings of scenes with shadows[END_REF]] which has initially been proposed to reduce the combinatory associated with line labeling of three-dimensional scenes. It has proved its effectiveness in solving some control problems such as identification, filtering and robust control.

A contractor associated to a set X is an operator C which associates to a box [x] ∈ R n another box C ([x]) ∈ R n such that the two following properties are always satisfied [START_REF] Jaulin | Applied interval analysis[END_REF]:

{ C([x]) ⊂ [x] (contractance property) C([x]) ∩ X = [x] ∩ X (completeness property)

Damage estimation and RUL prediction

In the sequel, we propose to use interval techniques for estimating the parameters of system (4.1). Thus, the following CSP is formulated:

C : { C k : x k = f (x k-1 , θ, w k ), y k = h(x k , u k ), k = k 0 , . . . k p θ ∈ [θ, θ] (4.4)
The characterization of the solution set Θ defined by the above CSP is based on set inversion techniques and is performed here with a Matlab R ⃝ toolbox called VSIVIA.

In the following, the general methodology for damage estimation and RUL prediction is given, as well as the basic notions of set inversion and VSIVIA.

Methodology.

During the first step which consists in estimating the parameters of the current degradation state with the measurements, inner and outer approximations of the parameter set Θ defined by (4.1) are computed. They are respectively denoted as Θ in and Θ out and satisfy:

Θ in ⊆ Θ ⊆ Θ out (4.5)

In the second step, the degradation prediction is computed with the natural inclusion function of the degradation state equation (2.1), where [θ] is computed through a projection of Θ out on the axes θ i . This inclusion function is given by

[x k ] = [f ]([x k-1 ], [θ], [-W, W ]) (4.6)
where [-W, W ] is the feasible domain of the disturbances w. A consistency check of the predicted intervals [x k ] (k > k p ) and the degradation threshold is used to estimate the remaining useful life. Due to uncertainty propagation, it is not possible to compute a reliable point estimation of the RUL. In the following, we propose to define lower and upper bounds of the RUL (i.e.

k RU L ∈ [k RU L , k RU L ]) as: { k RU L = k f -k p k RU L = k f -k p (4.7)
Similarly to the stochastic case where the RUL is characterized by a probability density function, in the set-membership context, the RUL can be considered as a random variable with an uniform pdf within the bounds [k RU L , k RU L ].

Set inversion and VSIVIA. Set inversion is one of the most used technique that was developed in the context of set characterization [START_REF] Moore | Interval analysis[END_REF][START_REF] Jaulin | Set inversion via interval analysis for nonlinear bounded-error estimation[END_REF].

Set inversion of a given box Y by a vector function f aims at characterizing the set defined by:

X = f -1 (Y) = {x ∈ R n | f (x) ∈ Y} (4.8)
Usually, the Set Inverter via Interval Analysis (SIVIA) [START_REF] Jaulin | Set inversion via interval analysis for nonlinear bounded-error estimation[END_REF]] algorithm is used to approximate the inner and outer approximations X in and X out of the set X defined by equation (4.8), given a set Y ∈ R m , an inclusion function [f ](•) : IR n → IR m for f , and an inclusion test [t](•) defined by:

[t](x) = ⎧ ⎪ ⎨ ⎪ ⎩ true if [f ]([x]) ⊂ Y, f alse if [f ]([x]) ∩ Y = ∅, undecided otherwise. (4.9)
The characterization of the set X with SIVIA is based on three main steps: (i) choice of an initial box [x] 0 supposed to contain at least a part of X, (ii) inclusion tests and (iii) bisection.

The inclusion test aims at checking if the considered box [x] 0 belongs or not to X, in function of the output of the inclusion test [t], and classifies [x] 0 into one of these three categories:

• [x] 0 is feasible if [t](x) = true, • [x] 0 is infeasible if [t](x) = f alse, • [x] 0 is indeterminate if [t](x) = undecided.
At the end of this inclusion test, the interval box [x] 0 is stored into a list. In SIVIA, the list that contains the feasible boxes is denoted S, the infeasible boxes are stored into N . As for the indeterminate boxes, they are bisected, stored into a list L and are subject to inclusion test and bisections until their size gets smaller than ϵ, which defines the precision of the outer approximation X out . They are then stored into the list of ambiguous boxes E such that X out = X in ∪ E. This procedure is repeated until the list L of indeterminate boxes whose size is bigger than ϵ becomes empty. The computational complexity of SIVIA algorithm exponentially increases with the parameter vector dimension. In order to reduce the number of bisections, SIVIA is used with interval constraint propagation techniques.

In this work, a Matlab R ⃝ toolbox called VSIVIA was used instead of SIVIA. VSIVIA is a vector implementation of SIVIA, which was designed so as to benefit from the advantages of Matlab R ⃝ in terms of vectorization. Therefore, the basic principle of VSIVIA is the same as SIVIA, except that VSIVIA performs the calculations vectorially. Indeed, VSIVIA evaluates all the boxes of the list L in a vector way instead of processing them one by one. Thus, the vector inclusion test [t](L) is defined by:

[t](L) = ⎧ ⎪ ⎨ ⎪ ⎩ in if [f ](L) ⊂ Y out if [f ](L) ∩ Y = ∅ undecided ¬in ∧ ¬out (4.10)
The reader may refer to [START_REF] Herrero | An efficient implementation of the sivia algorithm in a high-level numerical programming language[END_REF] to have a detailed algorithm of VSIVIA as well as application examples.

Each time VSIVIA is used in this chapter, the set of inputs required by the algorithm will be defined as:

Inputs = {[x] 0 , Y, ϵ} (4.11)
where [x] 0 is the initial box that contains the set to characterize, Y is the set to be inverted and ϵ defines the precision of the algorithm.

In the next section, the proposed methodology is applied to the academic test case presented in Chapter 2 where the prognosis of a fatigue crack growth with the Paris' law is performed. The same performance evaluation metrics are used to evaluate the algorithm and to compare it to the EKF and PF presented in Chapter 2.

Simulation results

In this section, the degradation that is considered is a fatigue crack propagation in an infinite rectangular plate. The case study is exactly the same as the one presented in Chapter 2. Indeed, the same measurement data and the same following discretized model are used:

a k = e C k (∆σ √ πa k-1 ) m k ∆N + a k-1 (4.12)
where a is the crack size, N is the number of cycles and ∆σ is the stress range whose value is known. C and m are the unknown model parameters. the user may choose between the two methods, knowing that the interval method is easier to implement despite its lower precision.

Conclusion

The proposed methodology based on interval techniques allows a guaranteed estimation and prediction of the degradation state. Its efficiency was demonstrated with an illustration example, which is a nonlinear fatigue crack growth. This method can be applied directly in the case where the degradation state is measured. However, this is not always the case. Thus, an interval based methodology that allows damage prognosis when the degradation state is not directly measurable is presented in the next section.

Unknown input interval observer for prognosis of discrete-time linear systems

In some cases, the degradation of a system may not be directly measurable. In this context, [Gucik-Derigny et al. 2011a] has introduced a model-based prognosis methodology that consists in considering the hidden degradation as a disturbance that affects the system dynamics, and is therefore referred to as an unknown input. The authors have used a nonlinear high gain unknown input observer to estimate the state and the unknown input which is the degradation state. Then, the parameters of the degradation model were identified with a least mean square algorithm. Finally, based on a known future input model, the RUL of the system was deduced when the failure threshold was reached. The authors have illustrated their methodology on an electromechanical oscillator. However, uncertainties were not considered.

Based on this previous study, this section considers the case where the degradation cannot be directly measured but only detectable through the drift of a model parameter. This parameter drift is viewed as an unknown input, whose reconstruction allows the estimation of the degradation signal by solving an inverse problem. Moreover, model and measurement uncertainties are considered in a set-membership framework. The RUL calculation is composed of three main steps. First, the bounds of the state and the unknown input are estimated, while taking the model and measurement noises into account. For this purpose, a novel methodology to synthesize an interval unknown input observer (UIO) for linear time-invariant (LTI) discrete-time systems is proposed. The bounds of the degradation signal are then deduced by solving an inverse problem based on interval arithmetics. The second step consists in estimating the feasible domain of the parameters of the degradation model based on a set-membership methodology using constraint satisfaction. Finally, as a third step, interval propagation is performed to estimate the bounds that include the future degradation state until failure threshold is reached, allowing the deduction of RUL interval. In order to demonstrate the efficiency of the proposed model-based prognosis methodology, the degradation of a suspension system is studied. A crack that cannot be directly measured is evolving in the suspension spring, and is affecting its stiffness constant, which is a parameter of the suspension dynamical model.

Unknown input interval observer for joint state and unknown input estimation

In this section, the methodology to synthesize an interval UIO for LTI discrete systems is presented.

Chapter 4. Set-membership framework for model-based prognosis

Consider the following system:

⎧ ⎨ ⎩ x(k + 1) = Ax(k) + Bu(k) + F d(k) + ω(k) y(k) = Cx(k) + δ(k) (4.14)
where x ∈ R n , u ∈ R m and y ∈ R p are respectively the state, the input and the measurement vectors; d ∈ R q is the unknown input vector which does not affect the outputs. A, B, C and F are constant matrices of appropriate dimensions. Finally, ω ∈ R n and δ ∈ R p are the state and measurement noises which are assumed to be bounded with a priori known bounds |ω|≤ ω and |δ|≤ δ where ω ∈ R n and δ ∈ R p are constant component-wise positive vectors and |•| is the component-wise absolute value for vectors. Moreover, it is assumed that n ≥ q and p ≥ q.

In order to solve the problem of state and unknown input estimation, unknown input observers (UIO) have been developed since 1970's [START_REF] Meditch | Observers for systems with unknown and inaccessible inputs[END_REF][START_REF] Hostetter | Observers for systems with unknown, unmeasurable inputs[END_REF][START_REF] Wang | Observing the states of systems with unmeasurable disturbances[END_REF], and are mostly used in the field of fault detection [START_REF] Frank | Survey of robust residual generation and evaluation methods in observer-based fault detection systems[END_REF][START_REF] Chen | Design of unknown input observers and robust fault detection filters[END_REF]]. However, in presence of measurement noise or uncertain parameters, classical Luenberger-based observers face some limitations. This is why interval observers have been proposed, to cope with uncertainties by evaluating the set of admissible values of the state at each time instant and computing the lower and upper bounds.

The problem of state estimation without unknown inputs using interval observers has been widely studied. In [START_REF] Gouzé | Interval observers for uncertain biological systems[END_REF], it was shown that in the particular case of asymptotically stable and cooperative systems (i.e systems where the Jacobian matrix of the state vector field has non-negative off-diagonal elements), interval observers can be designed directly. This assumption of cooperativity is the main limitation of interval observers as most of the systems are not cooperative. However, in the case of LTI systems this hypothesis can be relaxed by using a time-varying change of coordinates [START_REF] Mazenc | Interval observers for linear time-invariant systems with disturbances[END_REF] or by time-invariant ones [START_REF] Raïssi | Interval state estimation for a class of nonlinear systems[END_REF]].

In the set-membership framework, joint state and unknown input estimation has been considered but only for continuous-time systems [Gucik-Derigny et al. 2016b]. Standard discrete-time UIOs have already been used for state estimation in the presence of unknown inputs [START_REF] Valcher | State observers for discrete-time linear systems with unknown inputs[END_REF][START_REF] Darouach | Functional observers for systems with unknown inputs[END_REF], and also for unknown input estimation [START_REF] Maquin | Estimation of unknown inputs in linear systems[END_REF]]. Based on these works, the goal of this section is to establish an interval observer to jointly estimate the state and the unknown input of systems described by (4.14). The idea is to use the results from [START_REF] Maquin | Estimation of unknown inputs in linear systems[END_REF] to obtain an unknown input-free subsystem. A second state transformation using a time-invariant change of coordinates is performed in order to ensure the cooperativity property of the observation error in the new coordinates. Then, an interval observer is designed in the new coordinates and allows to deduce lower and upper bounds for the state in the original basis. Finally, the bounds of the unknown input are computed. An extension of these results to the case where the known and unknown inputs also affect the measurement equation is available in Appendix A.

Preliminaries

Notations and definitions

• Given a matrix M ∈ R n×m , define M + = max(0, M ), M -= M -M + , M * = [M + M -] and M = [ M 0 0 M ] . • For two vectors x 1 , x 2 ∈ R n or matrices M 1 , M 2 ∈ R n×n , the operators ≤, <, >
are understood component-wise.

• For x, x ∈ R n with x ≤ x, define X T = [ x x ] and X T = [ x x
] .

• A matrix M ∈ R n×n is called Schur stable if its spectral radius is less than one.

• A matrix M ∈ R n×n is called nonnegative if all its elements are nonnegative.

Lemma 1 [START_REF] Hirsch | Monotone maps: a review[END_REF]; Consider the linear system

x(k + 1) = Ax(k) + ω(k) (4.15)
where ω ∈ R n + and A is a nonnegative matrix. Then, ∀k > 0, x(k) ≥ 0 provided that x(0) ≥ 0. Such dynamical systems are called cooperative.

The following lemma was stated and proven in the continuous-time case in [START_REF] Gouzé | Interval observers for uncertain biological systems[END_REF] and is derived here in the discrete-time context.

Lemma 2 Consider the system described by (4.15) and suppose that the following assumptions are fulfilled:

• The matrix A is Schur stable and nonnegative;

• ω is bounded by a fixed positive vector Ω;

• x(0) ≤ x(0) ≤ x(0).
Then, the state vector x is asymptotically lower than the positive vector

ρ x = (I -A) -1 Ω (4.16)
Proof Starting from (4.15), it is shown that

x(k) = A k x(0) + k-1 ∑ i=0 A i ω(k -1 -i) (4.17)
Then, as ω ≤ Ω, we can deduce that ∀k ∈ Z + :

x(k) ≤ A k x(0) + k-1 ∑ i=0 A i Ω (4.18)
A is Schur stable, therefore the term A k x(0) converges to 0. Moreover, using classical results on the convergence of geometric series, the fact that A is Schur stable allows to deduce the convergence of

k-1 ∑ i=0 A i to ∞ ∑ i=0 A i = (I -A) -1 (4.19)
Finally, we can conclude that ∀k ∈ Z + :

x(k) ≤ (I -A) -1 Ω. (4.20)

□

Interval observers and cooperative discrete-time linear systems

Consider the following system without unknown inputs:

⎧ ⎨ ⎩ x(k + 1) = Ax(k) + Bu(k) y(k) = Cx(k) (4.21)
The aim of an interval observer is to compute two trajectories x(k) and x(k) such that x(k) ≤ x(k) ≤ x(k) for all k ∈ Z + with an initial condition verifying x(0) ≤ x(0) ≤ x(0). The upper and lower bounds could be obtained with Luenbergerbased observers defined by:

⎧ ⎨ ⎩ x(k + 1) = (A -LC)x(k) + Bu(k) + Ly(k) x(k + 1) = (A -LC)x(k) + Bu(k) + Ly(k) (4.22)
The dynamics of the errors e(k) = x(k)x(k) and e(k) = x(k)x(k) are given

Unknown input interval observer for prognosis of discrete-time linear systems

by:

⎧ ⎨ ⎩ e(k + 1) = (A -LC)e(k) e(k + 1) = (A -LC)e(k) (4.23)
Based on Lemma 1, the observation errors (4.23) are always positive and bounded if and only if the matrix (A -LC) is Schur stable and nonnegative [START_REF] Smith | Monotone dynamical systems: an introduction to the theory of competitive and cooperative systems[END_REF]; Efimov et al. [2013b]).

In the following sections, interval observers whose structures are similar to (??) will be designed to compute lower and upper bounds for the state vector even in the presence of noise and disturbances. It will also be shown that the assumption on the non-negativity of (A -LC) can be relaxed by using some changes of coordinates. Finally, an original approach will be presented to estimate the bounds of the unknown inputs.

Problem statement

The methodology proposed to jointly estimate the state and the unknown inputs of a LTI discrete-time system is splitted in two steps. First, two bounds x k ,x k ∈ R n for the state are estimated. Then, a technique to build the lower and upper bounds d k ,d k ∈ R q for the unknown input is described.

The unknown input interval observer developed here is based on the UIO proposed in [START_REF] Maquin | Estimation of unknown inputs in linear systems[END_REF] whose methodology is extended to systems with bounded disturbances and noise. Following a change of coordinates, the state is divided into two subsytems, one affected by the unknown input and the second one is unknown input-free. This allows to design an interval observer in the new coordinate basis to estimate the upper and lower bounds of the state. Then, by returning into the initial coordinates, the upper and lower bounds for the unknown input can be computed. First of all, the following assumption is required.

Assumption 1 C is a full row rank matrix and F is a full column rank matrix.

Under Assumption 1, there exists an orthogonal matrix H ∈ R n×n and matrices R 0 ∈ R q×q and K ∈ R q×q such that:

F = H [ R 0 0 ] K ⊤ (4.24)
This leads to the transformation of the system (4.14) into an equivalent one:

Chapter 4. Set-membership framework for model-based prognosis

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ z(k + 1) = Ãz(k) + Bu(k) + ⎡ ⎣ R 0 0 ⎤ ⎦ d(k) + ω(k) y(k) = Cz(k) + δ(k) (4.25)
where:

H = [ H 11 H 12 H 21 H 22 ] , Ã = H ⊤ AH = [ Ã11 Ã12 Ã21 Ã22 ] , B = H ⊤ B = [ B1 B2 ] , C = CH = [ C1 C2 ] z(k) = H ⊤ x(k) = [ z 1 (k) z 2 (k) ] , d(k) = K ⊤ d(k), ω(k) = H ⊤ ω = [ ω1 (k) ω2 (k)
] H ⊤ is supposed to be bounded, therefore |ω|≤ ω where ω is a constant positive vector. The system (4.25) is decomposed into an unknown input depending subsystem and an unknown input-free subsytem described by:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ z 1 (k + 1) = Ã11 z 1 (k) + Ã12 z 2 (k) + B1 u(k) + R 0 d(k) + ω1 (k) z 2 (k + 1) = Ã21 z 1 (k) + Ã22 z 2 (k) + B2 u(k) + ω2 (k) y(k) = C1 z 1 (k) + C2 z 2 (k) + δ(k) (4.26)
where C1 ∈ R p×q and C2 ∈ R p× (n-q) . C1 is supposed to be a full column rank matrix [START_REF] Hou | Design of observers for linear systems with unknown inputs[END_REF] and can be decomposed as: (p-q) ) and ỹ(k) = H ⊤ 1 y(k); the measurements equation can be decomposed as

C1 = H 1 [ R 1 0 ] K ⊤ 1 (4.27) with H 1 = [ H 011 H 012 ] (H 011 ∈ R p×q and H 012 ∈ R p×
⎧ ⎨ ⎩ ỹ1 (k) = R 1 K ⊤ 1 z 1 (k) + H ⊤ 011 C2 z 2 (k) + H ⊤ 011 δ(k) ỹ2 (k) = H ⊤ 012 C2 z 2 (k) + H ⊤ 012 δ(k) = C 2 z 2 (k) + H ⊤ 012 δ(k) (4.28) As ỹ1 (k) = G ⊤ s ỹ(k) with G ⊤ s = [ I q O q×(p-q) ]
, the expression of z 1 is extracted from (4.28):

z 1 (k) = E(y(k) -C2 z 2 (k) -δ(k)) (4.29) with E = K 1 R -1 1 G ⊤ s H ⊤ 1 .

State estimation

The state estimation is first performed in the coordinates r 2 . In the sequel, we define ∆

⊤ = [ δ -δ ] , ∆ ⊤ = [ -δ δ ] and Ω ⊤ = [ ω -ω ] , Ω ⊤ = [ -ω ω ] .
The following theorem allows to carry out an interval state estimation in the coordinates r 2 .

Theorem 1 Assume that r 2 (0) ≤ r 2 (0) ≤ r 2 (0). Then, for all k ∈ Z + the estimates r 2 (k) and r 2 (k) given by

⎧ ⎨ ⎩ r 2 (k + 1) = Rr 2 (k) + P B 2 u(k) + M y(k) + (-M ) * ∆ + P * Ω2 r 2 (k + 1) = Rr 2 (k) + P B 2 u(k) + M y(k) + (-M ) * ∆ + P * Ω2 (4.33)
are bounded and verify

r 2 (k) ≤ r 2 (k) ≤ r 2 (k) (4.34)
In addition, if the gain L is chosen such that (A 2 -LC 2 ) is Schur stable, then r 2 and r 2 are bounded.

Proof There are two results to prove:

1. ∀k ∈ Z + , r 2 (k) ≤ r 2 (k) ≤ r 2 (k)
i.e the upper and lower observation errors are positive.

2. Stability of the interval observer.

Step 1: Positivity of the observation errors.

The upper and lower observation errors are defined as

⎧ ⎨ ⎩ e r 2 (k) = r 2 (k) -r 2 (k) e r 2 (k) = r 2 (k) -r 2 (k)
The dynamics of the errors are given by

⎧ ⎨ ⎩ e r 2 (k + 1) = Re r 2 (k) + V 1 (k) e r 2 (k + 1) = Re r 2 (k) + V 2 (k)
where

V 1 (k) = (-M ) + (δ -δ(k)) -(-M ) -(δ + δ(k)) + P + (ω 2 -ω2 (k)) -P -(ω 2 + ω2 (k))
and

V 2 (k) = (-M ) + (δ + δ(k)) -(-M ) -(δ -δ(k)) + P + (ω 2 + ω2 (k)) -P -(ω 2 -ω2 (k))
According to Lemma 3, the matrix R = P (A 2 -LC 2 )P -1 is nonnegative.

For a given matrix M ∈ R n×m we have defined in Section 3.2 that M = M + +M -, therefore M + and -M -are nonnegative matrices. Moreover, |δ|≤ δ and |ω 2 |≤ ω2 . Therefore, as V 1 and V 2 are the sums of positive terms, we can deduce that ∀k ∈

Z + , V 1 (k) ≥ 0 and V 2 (k) ≥ 0.
Adding to that, as r 2 (0) ≤ r 2 (0) ≤ r 2 (0), we have e r 2 (0) ≥ 0 and e r 2 (0) ≥ 0. Therefore Lemma 1 allows to conclude that the upper and lower observation errors are always positive.

Step 2: Stability and convergence of the interval observer.

Defining E ⊤ r 2 = [ e r 2 e r 2 ]
leads to the following dynamics:

E r 2 (k + 1) = RE r 2 (k) + V (k) (4.35)
where R is defined in Section ?? as

R = P (A 2 -LC 2 )P -1 and V ⊤ = [ V 1 V 2 ] . As |δ|≤ δ and |ω 2 |≤ ω2 , it is straightforward to prove that V 1 ≤ V and V 2 ≤ V with: V = 2((-M ) + -(-M ) -)δ + 2(P + -P -)ω 2
Moreover, Lemma 3 states that R is Schur stable and nonnegative, therefore we can ensure the stability of the observation error dynamics (4.35).

Finally, with Lemma 2, we can deduce that E r 2 is asymptotically lower than the non-negative vector

ρ r 2 = (I -R) -1 V . (4.36)
□ Furthermore, since r 2 = P z 2 , the bounds of z 2 (k) are given by the following corollary.

Corollary 1 Under the conditions of Theorem 1, we have z

2 (k) ≤ z 2 (k) ≤ z 2 (k) with ⎧ ⎨ ⎩ z 2 (k) = (P -1 ) + r 2 (k) + (P -1 ) -r 2 (k) z 2 (k) = (P -1 ) + r 2 (k) + (P -1 ) -r 2 (k) (4.37) Proof We have P -1 z 2 (k) = ((P -1 ) + + (P -1 ) -)z 2 (k). If ∀k ∈ Z + , r 2 (k) ≤ r 2 (k) ≤ r 2 (k), then ((P -1 ) + r 2 (k) + (P -1 ) -r 2 (k) ≤ (P -1 )r 2 (k) ≤ (P -1 ) + r 2 (k) + (P -1 ) -r 2 (k) which is equivalent to z 2 (k) ≤ z 2 (k) ≤ z 2 (k).
Finally, as r 2 (k) and r 2 (k) are bounded, z 2 (k) and z 2 (k) are bounded as well.

□ The last step consists in computing the bounds for the whole state in the original coordinates

x ⊤ = [ x 1 x 2 ] and x ⊤ = [ x 1 x 2 ]
. Based on Theorem 1 and Corollary 1, the following theorem ensures the interval estimation of the state x.

Theorem 2 Assume that the conditions of Theorem 1 are satisfied and x(0) ≤ x(0) ≤ x(0). Then, for all k ∈ Z + the estimates x(k) and x(k) given by

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ x 1 (k) = H 11 Ey + (H 12 ) * Z 2 (k) + (-E 1 ) * Z 2 (k) + (-H 11 E) * ∆ x 1 (k) = H 11 Ey + (H 12 ) * Z 2 (k) + (-E 1 ) * Z 2 (k) + (-H 11 E) * ∆ x 2 (k) = H 21 Ey + (H 22 ) * Z 2 (k) + (-E 2 ) * Z 2 (k) + (-H 21 E) * ∆ x 2 (k) = H 21 Ey + (H 22 ) * Z 2 (k) + (-E 2 ) * Z 2 (k) + (-H 21 E) * ∆ (4.38) are bounded and verify x(k) ≤ x(k) ≤ x(k) (4.39) with E 1 = H 11 E C2 and E 2 = H 21 E C2 .
Proof We have x = Hz:

[ x 1 x 2 ] = [ H 11 H 12 H 21 H 22 ] [ E(y(k) -C2 z 2 (k) -δ(k)) z 2 (k) ] = [ H 11 Ey(k) + H 12 z 2 (k) -E 1 z 2 (k) -H 11 Eδ(k) H 21 Ey(k) + H 22 z 2 (k) -E 2 z 2 (k) -H 21 Eδ(k) ]
The observation errors relative to the state x are given by:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ e x 1 (k) = x 1 (k) -x 1 (k) e x 1 (k) = x 1 (k) -x 1 (k) e x 2 (k) = x 2 (k) -x 2 (k) e x 2 (k) = x 2 (k) -x 2 (k) 92 4.

Unknown input interval observer for prognosis of discrete-time linear systems

and in a developed form we get:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ e x 1 (k) = (H + 12 + (-E 1 ) + )e z 2 (k) -(H - 12 + (-E 1 ) -)e z 2 (k) +(-H 11 E) + (δ -δ(k)) -(-H 11 E) -(δ + δ(k)) e x 1 (k) = (H + 12 + (-E 1 ) + )e z 2 (k) -(H - 12 + (-E 1 ) -)e z 2 (k) +(-H 11 E) + (δ + δ(k)) -(-H 11 E) -(δ -δ(k)) e x 2 (k) = (H + 22 + (-E 2 ) + )e z 2 (k) -(H - 22 + (-E 2 ) -)e z 2 (k) +(-H 21 E) + (δ -δ(k)) -(-H 21 E) -(δ + δ(k)) e x 2 (k) = (H + 22 + (-E 2 ) + )e z 2 (k) -(H - 22 + (-E 2 ) -)e z 2 (k) +(-H 21 E) + (δ + δ(k)) -(-H 21 E) -(δ -δ(k)) (4.40)
With the same reasoning as in the proof of Theorem 1, it can be deduced from (4.40) that the state observation errors are positive. Therefore,

x(k) ≤ x(k) ≤ x(k) , ∀k ≥ k 0 .
By defining the compact error E x as

E x = [ e ⊤ x 1 e ⊤ x 1 e ⊤ x 2 e ⊤ x 2
] ⊤ and the matrices

J 1 = ⎡ ⎢ ⎢ ⎢ ⎣ H + 12 + (-E 1 ) + H - 12 + (-E 1 ) - H - 12 + (-E 1 ) -H + 12 + (-E 1 ) + H + 22 + (-E 2 ) + H - 22 + (-E 2 ) - H - 22 + (-E 2 ) -H + 22 + (-E 2 ) + ⎤ ⎥ ⎥ ⎥ ⎦ J 2 = ⎡ ⎢ ⎢ ⎢ ⎣ (-H 11 E) + -(-H 11 E) - -(-H 11 E) -+ (-H 11 E) + (-H 21 E) + -(-H 21 E) - -(-H 21 E) -+ (-H 21 E) + ⎤ ⎥ ⎥ ⎥ ⎦
it can be deduced that

E x ≤ J 1 E z 2 + 2J 2 δ.
Moreover, from Corollary 1, we have

E z 2 = RE r 2 with R = [ (P -1 ) + -(P -1 ) - -(P -1 ) -(P -1 ) + ] . It follows that E x ≤ J 1 RE r 2 + 2J 2 δ.
Finally, it is shown that E x is asymptotically element-wise lower than the nonnegative vector ρ x = J 1 Rρ r 2 + 2J 2 δ. □

Unknown input estimation

In this subsection, the upper and lower bounds of the unknown input d will be estimated. The expression of d is expressed from the first equation of (4.26):

(4.41) d(k) = KR -1 0 [z 1 (k + 1) -Ã11 z 1 (k) -Ã12 z 2 (k) -B1 u(k) -ω1 (k)]
By replacing z 1 with its expression in (4.29), equation (4.41) becomes:

(4.42) d(k) = KR -1 0 [Ey(k + 1) -E C2 z 2 (k + 1) -Eδ(k + 1) -Ã11 (Ey(k) -E C2 z 2 (k) -Eδ(k)) -Ã12 z 2 (k) -B1 u(k) -ω1 (k)]
The following theorem ensures the interval estimation of the unknown input d.

Theorem 3 Assume that the conditions of Theorem 1 are satisfied. Then, for all

k ∈ Z + the estimates d(k) and d(k) given by ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ d(k) = QEy(k + 1) -Q Ã11 Ey(k) -Q B1 u(k) + G * 1 Z 2 (k + 1) + G * 2 Z 2 (k) +G * 3 ∆ + G * 4 ∆ + G * 5 Ω1 d(k) = QEy(k + 1) -Q Ã11 Ey(k) -Q B1 u(k) + G * 1 Z 2 (k + 1) + G * 2 Z 2 (k) +G * 3 ∆ + G * 4 ∆ + G * 5 Ω1 (4.43)
are bounded and verify

d(k) ≤ d(k) ≤ d(k) (4.44) With Q = KR -1 0 , G 1 = -QE C2 , G 2 = Q( Ã11 E C2 -Ã12 ), G 3 = -QE, G 4 = Q Ã11 E and G 5 = -Q.
Proof The lower and upper observation errors of the unknown input d are:

⎧ ⎨ ⎩ e d (k) = d(k) -d(k) e d (k) = d(k) -d(k)
By developing these expressions we obtain:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ e d (k) = G + 1 e z 2 (k + 1) -G - 1 e z 2 (k + 1) + G + 2 e z 2 (k) -G - 2 e z 2 (k) +G + 3 (δ -δ(k + 1)) -G - 3 (δ + δ(k + 1)) + G + 4 (δ -δ(k)) -G - 4 (δ + δ(k)) +G + 5 (ω 1 -ω1 (k)) -G - 5 (ω 1 + ω1 (k)) e d (k) = G + 1 e z 2 (k + 1) -G - 1 e z 2 (k + 1) + G + 2 e z 2 (k) -G - 2 e z 2 (k) +G + 3 (δ + δ(k + 1)) -G - 3 (δ -δ(k + 1)) + G + 4 (δ + δ(k)) -G - 4 (δ -δ(k)) +G + 5 (ω 1 + ω1 (k)) -G - 5 (ω 1 -ω1 (k)) (4.45)
With the same reasoning as in the proof of Theorem 1, it is deduced from (4.45) that the unknown input observation errors are positive.

As e z 2 , e z 2 , δ and ω1 are bounded, then d and d are bounded as well.

Unknown input interval observer for prognosis of discrete-time linear systems

If we define E

⊤ d = [ e d e d
] and:

N 1 = [ G + 1 -G - 1 -G - 1 G + 1 ] , N 2 = [ G + 2 -G - 2 -G - 2 G + 2 ] , N 3 = G + 3 -G - 3 , N 4 = G + 4 -G - 4 , N 5 = G + 5 -G - 5 ,
then it can be deduced that

E d ≤ (N 1 + N 2 )E z 2 + 2(N 3 + N 4 ) [ δ δ ] + 2N 5 [ ω1 ω1
]

.

As E z 2 = RE r 2 , we obtain

E d ≤ (N 1 + N 2 )RE r 2 + +2(N 3 + N 4 ) [ δ δ ] + 2N 5 [ ω1 ω1
] .

Using the upper bound ρ r 2 of the observation error E r 2 , we can deduce that E d is asymptotically lower than the non-negative vector:

ρ d = (N 1 + N 2 )Rρ r 2 + 2(N 3 + N 4 ) [ δ δ ] + 2N 5 [ ω1 ω1
] . □

Model-based prognosis using an unknown input interval observer

In this section, a prognosis methodology using the previous unknown input interval observer is presented. This methodology could be applied to various linear systems for which a model is available. Here, remaining useful life prognosis is performed on a suspension system.

Presentation of the model

A passive suspension system is composed of springs that absorb the shocks by oscillating, and dampers that limit excessive suspension movement by reducing the oscillations of the springs. The suspension model that is used in this work is represented in Figure 4.4, and the corresponding mathematical model that describes the system dynamics was obtained from the Newton's second law: During its lifetime, a suspension spring is subject to repeated bending and sometimes overloading that tend to weaken it and eventually cause cracks. A crack evolving in the suspension spring may cause the system failure if a critical crack length is reached. Therefore, estimating and forecasting the crack length will allow to calculate the RUL of the system. The problem is that this crack length is not directly measurable. However, the propagation of this crack affects the value of the suspension stiffness k 1 , which tends to decrease when the length of the crack increases. Therefore, the degradation induced by the crack propagation can be tracked by estimating the evolution of the decrease of the suspension stiffness k 1 .

⎧ ⎨ ⎩ m 1 ẍ1 = -b 1 ẋ1 + b 1 ẋ2 -k 1 x 1 + k 1 x 2 m 2 ẍ2 = b 1 ẋ1 -b 1 ẋ2 + k 1 x 1 -k 1 x 2 -b 2 ẋ2 -k 2 x 2 + b 2 ẋr + k 2 x r (4.46) m 1 m 2 x 1 x 2 k 1 k 2 b 1 b 2
In order to take the degradation of k 1 into account, the following equation is introduced:

k 1 = k N 1 + δk 1 (4.47)
where k N 1 is the nominal value of k 1 and δk 1 is the decrease in the suspension stiffness value.

After replacing this expression of k 1 in (4.46), the new equations that consider 96 4.2. Unknown input interval observer for prognosis of discrete-time linear systems the evolution of k 1 are:

⎧ ⎨ ⎩ m 1 ẍ1 = -b 1 ẋ1 + b 1 ẋ2 -k N 1 x 1 + k N 1 x 2 -δk 1 (x 1 -x 2 ) m 2 ẍ2 = b 1 ẋ1 -b 1 ẋ2 + k N 1 x 1 -k N 1 x 2 -b 2 ẋ2 -k 2 x 2 + b 2 ẋr + k 2 x r + δk 1 (x 1 -x 2 ) (4.48)
A discrete-time state-space representation of the system is adopted. The state vector is defined as

x = [ x 1 x 2 ẋ1 ẋ2 ] ⊤ and the input vector is u = [ ẋr x r ] ⊤ .
Moreover, modeling and measurement uncertainties are taken into account in the form of additive noises that are respectively denoted as ω and δ. Therefore, (4.48) can be rewritten in the following state-space form:

⎧ ⎨ ⎩ x(k + 1) = Ax(k) + Bu(k) + F d(k) + ω(k) y(k) = Cx(k) + δ(k) (4.49)
where

A = ⎡ ⎢ ⎢ ⎢ ⎣ 0 0 1 0 0 0 0 1 -k N 1 /m 1 k N 1 /m 1 -b 1 /m 1 b 1 /m 1 k N 1 /m 2 -(k N 1 + k 2 )/m 2 b 1 /m 2 -(b 1 + b 2 )/m 2 ⎤ ⎥ ⎥ ⎥ ⎦ , B = ⎡ ⎢ ⎢ ⎢ ⎣ 0 0 0 0 0 0 b 2 /m 2 k 2 /m 2 ⎤ ⎥ ⎥ ⎥ ⎦ C = ⎡ ⎢ ⎣ 1 -1 0 0 0 0 1 0 0 0 0 1 ⎤ ⎥ ⎦ , F = ⎡ ⎢ ⎢ ⎢ ⎣ 0 0 -1/m 1 1/m 2 ⎤ ⎥ ⎥ ⎥ ⎦
and the unknown input d is given by:

d(k) = δk 1 (k)(x 1 (k) -x 2 (k)) (4.50)

Estimation of the RUL of the suspension system

The degradation of the suspension system is caused by the crack evolving in the spring of stiffness k 1 . Therefore, in order to estimate the RUL of the suspension system, the evolution of this crack should be estimated and then predicted. The crack length is not directly measurable, however it is related to the stiffness k 1 by a polynomial expression:

k 1 (ξ) = deg ∑ i=0 α i ξ i (4.51)
where ξ ∈ [0, 1] is the degradation measure, which is a quantity derived from the normalization of the crack length, and is described by a Paris' law. This polynomial allows to reflect the stiffness degradation which evolves slowly at the crack initiation and then at a faster rate when the crack approaches its critical length. The procedure to estimate the degree and the coefficients of this polynomial is described in [START_REF] Sobczyk | Stochastic dynamics with fatigue-induced stiffness degradation[END_REF].

Therefore, the goal is to track and predict the evolution of the degradation measure ξ. For this purpose, k 1 should be estimated. This can be done using:

k 1 = k N 1 + δk 1 (4.52)
where the expression of δk 1 is derived from equation (4.50):

δk 1 (k) = d(k) x 1 (k) -x 2 (k) (4.53)
Thus, the first step consists in estimating the unknown input d and the states x 1 and x 2 . The linear discrete-time unknown input interval observer synthesized in Section ?? is used for this. Moreover, it is considered that x 1 ̸ = x 2 and that the slowly evolving variable δk 1 is independent of x 1 and x 2 in order to keep a linear framework.

Estimation of x and d with the UIO

The numerical values of the model parameters were taken from [START_REF] Alvarez-Sanchez | Nonlinear observer-estimator application: quarter-car suspension system under road disturbances[END_REF]. They are listed in Table 4 which reflects the road excitation is modeled as a stochastic process (Figure 4.5). The state and measurement noises are assumed to be uniformly distributed and bounded with δ = 10 -3 [1 1 1] ⊤ and ω = 10 -4 [1 1 1] ⊤ .

Finally, the bounds of the initial state vector are chosen as

x 0 = [-1 -1 -1] ⊤ and x 0 = [1 1 1] ⊤ .
State estimation. Assumption 1 is verified since D is a full column rank matrix and C is a full row rank matrix. In addition, Assumption 2 holds as the pair (A 2 , C 2 ) is observable. The gain L is chosen with a Schur stable pole assignment {0.94 , 0.93 , 0.95} and is given by L T = ⎡ ⎢ ⎣ -0.0004 -0.0246 0.0572 0.3237 0.0015 0.0786

⎤ ⎥ ⎦.
Finally, Lemma 1 is satisfied with a transformation matrix P given by: allowing to obtain a Schur stable and nonnegative matrix R = P (A 2 -LC 2 )P -1 .

P = ⎡ ⎢ ⎣ -2.
The results for state x 1 and x 2 bounds estimation are depicted in Fig. 4.6. The results show that the proposed interval UIO is well suited for state estimation in discrete-time with the presence of bounded noises and unknown input.

In the following, the results obtained for the computation of the lower and upper bounds of the unknown input are presented.

Unknown input estimation. The estimation of the bounds of the unknown input d are required to estimate the stiffness degradation δk 1 . These bounds are computed using (4.43) and the results are shown in Figure 4.7, where a zoom is done in order to highlight the fact that the bounds actually contain the reference. It can be seen that the older the suspension system gets, the greater is the amplitude of the unknown input d. This indicates that a degradation is actually evolving. Once the bounds of the unknown input are obtained, the stiffness degradation δk 1 is estimated in the following.

Estimation of the stiffness degradation δk 1 and k 1

Interval arithmetics is applied to compute the bounds of δk 1 from the bounds of d and equation (4.53). This calculation is done in 3 steps:

1. Denote x 12 = x 1x 2 , then:

⎧ ⎨ ⎩ x 12 = x 1 -x 2 x 12 = x 1 -x 2 (4.55) 2. Denote x inv 12 = 1 x 1 -x 2 , then: ⎧ ⎨ ⎩ x inv 12 = min( 1 x 12 , 1 x 12 ) x inv 12 = max( 1 x 12 , 1 x 12 ) (4.56)
3. Finally we have:

⎧ ⎨
The computation of δk 1 and δk 1 requires numerical treatment. Indeed, their values will diverge when x 12 → 0 or x 12 → 0, which happens frequently as x 1 and x 2 are the vertical displacements of the suspension system. To overcome this problem, it is verified at each time step if 0 ∈ [x 12 , x 12 ], and if it is the case, the point is eliminated.

Then, the bounds of the stiffness k 1 are deduced with k 1 = k N 1 + δk 1 as the degradation measure ξ depends on k 1 and not on δk 1 . Figure 4.8 shows the evolution of k 1 .

A zoom of Figure 4.8 is displayed in Figure 4.9 to show that the reference is included in the estimated bounds. The effect of the numerical treatment to avoid the divergence of the bounds clearly appears.

Reconstruction of the degradation signal ξ

Using the reconstructed signal of k 1 , the bounds of the degradation signal ξ can be estimated with a polynomial whose expression was adapted from [START_REF] Luo | Model-based prognostic techniques applied to a suspension system[END_REF]:

k 1 (ξ) = k N 1 -8.84 × ξ + 4.34 × ξ 2 -1.66 × ξ 3 + 2.53 × ξ 4 -1.38 × ξ 5 (4.58)
where

Inputs = {[x] 0 = [0, 1] , Y = [k 1 (k), k 1 (k)] , ϵ = 0, 1}
The resulting bounds of the degradation signal are shown in Figure 4.10.

Estimation and prediction of the bounds of the degradation state

Once the bounds that contain the degradation signal are available, the methodology presented in Section ?? to estimate and then forecast the bounds of the degradation state is finally applied.

For this purpose, a model that describes its evolution is required. The model that is used for the fatigue crack growth is a Paris' law, where ξ is the normalized crack length. Therefore, the discretized model is: Therefore, the set of input parameters of the VSIVIA algorithm required to characterize the parameter set Θ is:

ξ k+1 = ∆N e C k (∆σ √ πξ k ) m k + ξ k (4.
Inputs = {[x] 0 = [2, 5] × [-27 -19] , Y = {[ξ k , ξ k ]∀k ∈ [k 0 , k p ]} , ϵ = 0, 1}
In the same way as in Section 4.1, the bounds of the degradation state ξ are estimated with VSIVIA until the prediction time N p (Figure 4.11

[ξ k+1 ] = ∆N e [C k ] (∆σ √ π[ξ k ]) [m k ] + [ξ k ] (4.60)
where [m] and [log(C)] are derived from the projection of Θ out on the axes m and log(C).

The prediction of the degradation state ξ was performed for different prediction time cycles N p in order to highlight the fact that the more data is available, the more the predicted interval tends to shrink. The results are depicted from Figure 4.12 to Figure 4.14.

In can be seen that the estimated and predicted bounds of the degradation state ξ include the reference degradation state, and the interval is shrinking when more data is available.

The values of the interval bounds containing the RUL are depicted in Table 4.4. It can also be concluded that the bounds actually include the reference RUL value, and they tend to shrink over time. The performance of the proposed methodology is evaluated next in terms of accuracy, precision and timeliness.

Performance evaluation metrics

In this case, we have chosen R 0 = R min = R max = 0.5 × 10 6 cycles. The obtained metrics values are shown in Table 4.5. The same conclusion as before can be drawn concerning the fact that the more measurements are available, the better the metrics values get.

Conclusion

In this chapter, two set-membership methodologies adapted to model-based prognosis were presented.

In the first section, the case were the degradation state is directly measured was considered. The available measurements, the degradation model and an interval technique named set inversion were used to estimate the current degradation state. Then, from prediction time, the inclusion function and the estimated feasible parameter domain were used to prediction the future degradation state. The same academic fatigue crack growth case study presented in Chapter 2 was used to demonstrate the efficiency of the method.

In the second section, an interval observer for LTI discrete-time systems was proposed for the estimation of state and unknown input in the presence of bounded process and measurement noises. All the steps to design the interval UIO were described and the proofs of the main results were given. Finally, the use of the proposed interval UIO in the context of model-based prognosis was demonstrated on a suspension system subject to a fatigue crack growth. The obtained results have highlighted the efficiency of the method to perform on-line prognosis as the computational time is reasonable and the data received incrementally allow to update the predictions which get more and more precise.

To conclude, the ability of set-membership methods to perform on-line modelbased prognosis while taking uncertainty into account was demonstrated in this chapter.

Chapter 5 Inverse First-Order Reliability method for prognosis

In this chapter, another kind of uncertainty propagation method based on probability theory is presented. This method is an alternative to sampling-based methods and is called the inverse first-order reliability method (Inverse FORM). It is an analytical method originally developed in the structural reliability analysis field, which aims at evaluating the failure probability of a structure. Here, a method to compute the RUL probability density function (pdf) based on the Inverse FORM algorithm is presented. This methodology can be applied if an expression of the RUL or an algorithm whose output is the RUL is available.

In a first step, a filter is employed for the estimation of the current degradation state with the available measurements. Then, the Inverse FORM algorithm is used to find for which parameter values the system fails, with reference to a specified failure probability level. Finally, these parameter values and the explicit expression of the RUL allow to calculate the desired percentiles of the RUL pdf in terms of the cumulative distribution function according to the user's requirements. Indeed, one may compute only some points of the RUL pdf and then reconstruct the entire pdf by interpolation if needed. In the cases where only the mean of the RUL pdf and the 95% probability bounds are required, the Inverse FORM can be very efficient.

Therefore, the proposed methodology does not require huge computational efforts and is not time consuming. Indeed, the computation of the entire pdf of the RUL is not required, and there is no need to propagate the model equation step by step until the threshold is reached. Only the desired RUL pdf points can be computed, for example the mean value, the 5% bounds and so on.

For the sake of illustration, the crack growth model using the Paris' law presented in Chapter 2 is used again. This work has led to the publication of an international conference paper [Robinson, Marzat, and Raïssi 2017a].

Failure probability estimation. If f X is the joint pdf of X, then computing the failure probability of the system is equivalent to evaluating the probability that the output Y exceeds the threshold y th , or also to evaluating the probability that the realizations x of the input vector X are in the failure domain:

P f = P (Y ≥ y th ) = P (G(X) ≤ 0) = ∫ G(x)≤0 f X (x) dx.
(5.2)

Therefore, for a given threshold value y th , the corresponding value of P f can be computed. Moreover, the relationship between P f and the CDF of Y , denoted by F Y , evaluated at a given y th , is given by:

F Y (y th ) = P (Y ≤ y th ) = 1 -P f (5.3)
In this way, with several threshold values y th , the entire CDF F Y can finally be reconstructed.

If required, the pdf of Y (denoted f Y ) can be obtained with the following relationship:

d dy F Y (y) = f Y (y) (5.4)
The integral in equation (5.2) can be calculated analytically only in some simple academic cases, therefore numerical methods have been developed to approximate it. Available techniques for failure probability estimation can be split into two categories: one the one hand, the simulation methods, relying on the use of Monte Carlo simulations to compute the integral. Among others, one can cite the crude Monte Carlo sampling, importance sampling, subset sampling [START_REF] Robert | Monte Carlo statistical methods[END_REF]. On the other hand, a class of approximation methods, more specific to structural reliability analysis, has been proposed. They are known as the first and second-order reliability methods (FORM/SORM) [START_REF] Ditlevsen | Structural reliability methods[END_REF]. In this case, the function G(•) is approximated by its first/second-order Taylor series expansion so as to simplify the estimation problem and to reduce it to a geometric problem.

In this work, the first-order reliability method (FORM) is of interest for RUL calculation. Before presenting the methodology, the relationship between solving a reliability analysis problem and RUL computation is highlighted in the next subsection.

RUL calculation: an uncertainty propagation problem

In this section, the problem of RUL calculation is formulated as an uncertainty propagation problem in order to make a connection with reliability analysis.

During the current state estimation, sensor data are available for a specific observation interval whose size depends on the prediction time step k p . Then, from

Inverse First-Order Reliability Method for RUL calculation

First-order reliability method (FORM)

FORM is one of the most commonly used approximation methods. Indeed, its ease of implementation and low computation time make it a valuable tool.

Main steps of the FORM algorithm. One first needs to transform the uncertain parameters vector X into a vector of standard normal variables U and to work in this new standard normal space.

Then, the idea is to proceed to the linearization of the function G using a firstorder Taylor series approximation. This linearization is done around the so-called Most Probable Point (MPP), which is the point on the limit-state surface closest to the origin in the standard normal space. The MPP is found by solving an optimization problem aiming at minimizing the distance between the origin and the limit-state surface.

Once the MPP is identified, the reliability index β which is equal to the distance from the origin of the standard normal space to the MPP is computed.

Finally, the failure probability is obtained with:

P f ≈ P f,F ORM = Φ(-β) (5.7)
where Φ(•) denotes the standard normal CDF. If the true limit-state function G is linear, then this equation is exact. In the case where the true G function is nonlinear, the failure probability obtained from the above formula is only an approximation.

In the following, the use of the FORM algorithm to compute the RUL CDF/PDF is described.

Computation of the RUL CDF.

As it was explained in Section 5.1.1, one needs to apply the following steps to calculate the failure probability P f of a system and to deduce the CDF/PDF of the system output Y : 1. General case: define the computational model M that allows to calculate the quantity of interest Y . ⇒ RUL calculation: here, the quantity of interest is the RUL R and the model M may be an analytical expression of the RUL or a black-box function.

2. General case: identify the vector of random variables X which contains the model input uncertainties (material properties, geometry, loading, etc.). ⇒ RUL calculation: the uncertain inputs involved in the RUL calculation problem at the prediction time k p are mainly the current degradation state estimate and the model parameter estimates at time k p , future loading, future process and measurement noises.

3. General case: choose a failure scenario, i.e define a threshold value y th . This allows to define the limit-state function G such that G(X) = y th -M(X) = y th -Y .

⇒ RUL calculation: in this case, the threshold value is denoted r th , then the limit-state function is defined as:

G(X) = r th -M(X) = r th -R (5.8)
4. General case: calculate P f,F ORM and deduce the corresponding CDF value F Y (y th ). ⇒ RUL calculation: the above FORM algorithm is applied to compute one value of P f,F ORM for a particular value r th :

P f,F ORM (r th ) = P (G(X) ≤ 0) = P (R ≥ r th ) (5.9)
Then, the corresponding RUL CDF value F R (r th ) is deduced according to equation (5.3):

F R (r th ) = P (R ≤ r th ) = 1 -P f,F ORM (r th ) (5.10)
Thus, the FORM algorithm should be applied for multiple values r (instead of one particular value r th ) to compute the corresponding values of P f,F ORM (r). At this point, one should notice that the values r under consideration can be seen either as varying threshold values, or as possible realizations of the random variable R (i.e. the unknown RUL). The values of F R (r) are then deduced using equation (5.10) to finally obtain the entire RUL CDF F R after a certain number of repetitions of FORM.

However, choosing the values of r for the calculation of P f,F ORM (r) is difficult as the RUL distribution is not known a priori. Yet, the values of P f,F ORM (r) are assumed to be included in the interval [0, 1] by definition, and so is η = 1 -P f,F ORM (r). Thus, an inverse strategy can be applied (Fig. 5.1): starting from different probability values of P f,F ORM (or CDF values η = F R (r) = 1 -P f,F ORM (r)), one can find the corresponding r values. This strategy is known in reliability literature as Inverse FORM algorithm. Inverse FORM is used to solve an inverse reliability problem. Indeed, reliability analysis aims at computing the failure probability P f of a system. In the case of inverse reliability analysis, the goal is to find the vector of uncertain inputs X associated to a given failure probability P f .

Inverse First-Order Reliability Method for RUL calculation

The main steps of both forward and inverse reliability problems in the context of RUL PDF/CDF computation are illustrated in Fig. 5.2 for clarity: variables we have:

RUL CDF computation with FORM -Generation of one realization x ∈ R s of X at time k p -Choice of a threshold value r FORM Computation of the RUL CDF value F R (r) = 1 -P f,F ORM (r) x ∈ R s , r Pf,FORM (r) RUL CDF computation with Inverse FORM Choice of a value η = F R (r) ∈ [0, 1] Inverse FORM -Computation of r = M(x) -r is the RUL CDF value such that F -1 R (η) = r Pf = 1 -η x ∈ R s
u j i = x j i -µ i σ i (5.12)
where µ and σ are respectively the mean and the standard deviation of the uncertain variables. These quantities can be derived from the estimation step, as an output of the filtering algorithm (e.g. a Kalman filter).

3. Calculation of the gradient vector of G:

α i = ∂G ∂u i = ∂G ∂x i × ∂x i ∂u i .
(5.13) 4. Calculation of the next point u j+1 :

u j+1 = -β α ||α|| (5.14)
where β = -Φ -1 (P f ).

5. Transformation into original space to compute x j+1 .

The steps from 3 to 5 are repeated until the following convergence criteria within thresholds δ 1 and δ 2 are satisfied:

(i) The solution belongs to the limit-state surface:

|G(x j )r th |≤ δ 1 .

(5.15) (ii) The solution is almost constant between two iterations:

|x j+1 -x j |≤ δ 2 .
(5.16)

The values δ 1 and δ 2 are very small and are fixed by the user. Here, they are equal to 10 -4 , which is the value that is usually chosen.

The above procedure can be repeated for different values of P f . According to the needs of the user, the entire CDF can be computed with P f = {0.1, 0.2, . . . , 0.9}, or only the 95% probability bounds and the mean can be computed with P f = {0.05, 0.5, 0.95}. [START_REF] Sankararaman | Analytical algorithms to quantify the uncertainty in remaining useful life prediction[END_REF] have proved that the pdf computed using Inverse FORM is is included in confidence bounds generated with 1000 Monte Carlo simulations.

Numerical results

In this section, the academic fatigue crack growth case study presented in Chapter 2 is used again. The filtering algorithm applied to estimate the current crack length with measurements is the EKF algorithm. Then, the proposed Inverse FORM methodology is used to compute the RUL at the prediction time cycle when there is no more crack length data.

Explicit model of the RUL

The value of the RUL calculated at a prediction time cycle N p is given by the following expression: R(N p ) = N f -N p (5.17)

where N f is the failure time cycle i.e the cycle at which the crack size reaches its critical value a f . Thus, in order to obtain an analytical expression of the RUL, the expression of N f is required.

First of all, the degradation model of the system is needed. In this case, the degradation is the fatigue crack growth and the Paris' law presented in Chapter 2 is used:

da dN = C(∆K) m , ∆K = ∆σ √ πa (5.18)
where a is the crack size, N is the number of cycles, ∆K is the range of stress intensity factor and ∆σ is the stress range. C and m are the unknown model parameters.

The RUL is calculated from the prediction time cycle N p where the crack length is a p until the failure time cycle N f where the crack length is a f . Therefore, the differential equation 5.18 is solved in the intervals [a p , a f ] and [N p , N f ] to obtain the expression of N f : This differential equation is solved to find the expression of the crack length a with respect to the cycle number N :

C(∆σ √ π) m dN = 1 ( √ a) m da C(∆σ √ π) m dN = a -m 2 da ∫ N f Np C(∆σ √ π) m dN = ∫ a f ap a -m 2 da C(∆σ √ π) m (N f -N p ) = ( 2 2 -m ) ( a ( 2-m 2 ) f -a ( 2-m 2 ) p ) .
This yields:

N f = ( 2 2 -m ) ⎛ ⎜ ⎝ a ( 2-m 2 ) f -a ( 2-m 2 ) p C(∆σ √ π) m ⎞ ⎟ ⎠ + N p . (5.19)
Finally, the following expression of the RUL calculated at the prediction time N p is obtained:

R(N p ) = ( 2 2 -m ) ⎛ ⎜ ⎝ a ( 2-m 2 ) f -a ( 2-m 2 ) p C(∆σ √ π) m ⎞ ⎟ ⎠
(5.20)

The next step consists in identifying the random variables of the problem. Here, it is considered that the values of the parameters m and C at time cycle N p are uncertain. Moreover, the estimated crack length at time cycle N p is also a source of uncertainty. Therefore, the uncertain parameters vector is

X = [a p , C p , m p ].
As explained in the previous section, the Inverse FORM is used to find for which parameter values the system fails, with reference to a specified failure probability level.

Simulation results

Synthetic data generation. The parameters used to generate data to simulate crack length measurements are exactly the same as those presented in Section 2.4.2. Thus, a new crack length data is obtained every ∆N = 50 cycles.

Estimation of the current crack length. The first step consists in estimating the current crack size with the available measurements until the prediction time cycle N p . Here, an EKF is used because of the low computational time and effort required for the implementation. Moreover, its efficiency as long as there are measurements was demonstrated in Section 2.4.3. The variances used for the tuning of the EKF are the same as in Section 2.4.2.

RUL computation.

From time cycle N p , the Inverse FORM algorithm described in Section 5.2.2 is applied to compute the RUL:

1. Setting j to 0 and initial guess of the realizations of the uncertain parameter vector x j = [x j 0 , . . . , x j s ] where s is the number of uncertain parameters. ⇒ The vector X ⊤ =

[ a C m ] ⊤ estimated by the EKF at time cycle N p is used as the initial guess.

2. Transformation into the standard normal space. In the case of Gaussian variables we have:

u j i = x j i -µ i σ i (5.21)
where µ and σ are respectively the mean and the standard deviation of the uncertain variables. ⇒ These quantities are derived from the values obtained with the EKF at time cycle N p .

3. Calculation of the gradient vector of G: .22) 4. Calculation of the next point u j+1 :

α i = ∂G ∂u i = ∂G ∂x i × ∂x i ∂u i . ( 5
u j+1 = -β target α ||α|| (5.23)
where β target = -Φ -1 (P f ).

⇒ The user specifies the value of P f and the value of β target is deduced.

5. Transformation into original space to compute x j+1 .

6. Steps 3 to 5 are repeated until the convergence criteria defined in Section 5.2.2 are satisfied.

This whole algorithm is repeated for each value of P f needed by the user. For example, one may require only the mean RUL value associated to P f = 0. The 1% and 99% bounds are large at the beginning, which means that the uncertainty in the RUL is high initially. However, these bounds tend to shrink as more data is available and the estimation becomes more and more accurate.

The use of the Inverse FORM to compute the mean RUL value and its confidence interval is very concluding because the computational time is very low and there was no need to compute the entire RUL pdf as with the particle filter. Moreover, with the particle filter, the crack propagation equation was propagated step by step from the prediction time cycle N p until the critical crack length was reached, increasing the computational time and effort.

Performance metrics are calculated to evaluate the accuracy, precision and timeliness of the Inverse FORM.

Numerical results

Performance evaluation metrics. The simulation to compute the mean RUL value associated to P f = 0.5 at N p = 1200 cycles is run 100 times to have a large enough sample for the performance evaluation and to compare with the results obtained with the PF and the EKF in Section 2.4.3 with the same test case.

The distribution of the mean RUL values of each simulation run is shown in FIg. 5.7. Moreover, the numerical values of the performance metrics are depicted in Table 5.1. putation in terms of accuracy and timeliness. However, the Inverse FORM is less successful than the interval method concerning the three metrics.

Nevertheless, although the EKF was used for the estimation of the current degradation state, the performance metrics of the Inverse FORM remain acceptable, especially if a high precision is not required.

Even if the particle filter has the best results, its main drawback is its complex implementation and its higher computational time as the entire pdf of the RUL must be computed. Moreover, a large number of the whole trajectories of the degradation state until the threshold is reached have been generated to deduce the failure time and then compute the RUL pdf (Fig. 2.3).

The Inverse FORM only needs to compute a few values of RUL associated to to P f ={0.01,0.5,0.99}. This is sufficient to obtain the mean value of the RUL with the 99% probability bounds. Therefore, the computational time with the Inverse FORM is much lower, and the performance metrics are satisfying and although lower than the particle filter and the interval method.

Conclusion

In this chapter, a general methodology to compute the RUL with the Inverse FORM was presented, then evaluated and compared to filter-based methods through different performance metrics. The results have highlighted that in terms of accuracy and timeliness the Inverse FORM gives better results than the EKF. Even if the Inverse FORM does not outperform the interval method and the PF, the calculated metrics remain reasonable. Moreover, it presents a significant advantage in terms of computational time.

Therefore, the Inverse FORM can be an alternative to PF-based methods and interval methods for RUL prognosis in the case where low computational cost algorithms are needed, and a high precision is not required. However, if highly nonlinear and complex models are involved, the particle filter remains the best approach to choose.

Conclusion and perspectives

Conclusion

In this thesis, different methods for model-based prognosis were developed. The main purpose was the inclusion of uncertainty sources in order to measure the uncertainty in the predicted RUL. The main uncertainty sources that were considered are model uncertainty, measurement uncertainty and future uncertainty related to loading conditions. In order to deal with these uncertainty sources, probabilistic methods, interval methods and reliability analysis methods were extended to modelbased prognosis. In order to compare the proposed methods with performance evaluation metrics (accuracy, precision and timeliness), they were each time applied to nonlinear fatigue crack growth prognosis based on the Paris' law, which is referred to as the Paris case study.

A literature review was first conducted in Chapter 1 in order to be aware of the existing model-based approaches, and to better identify the main challenges while developing model-based prognosis methods. It was concluded that filtering and uncertainty propagation techniques capable of quantifying model, measurement and future uncertainty into account are required to perform an accurate RUL prediction. Particle filter-based methods are the most widespread because of their ability to tackle different uncertainty sources, nonlinear systems and non-Gaussian noises. Moreover, performance evaluation metrics, namely accuracy, precision and timeliness, aiming at evaluating and comparing the performances of prognosis methods were presented.

In Chapter 2, all the steps of RUL prognosis using an extended Kalman filter (EKF) and a particle filter (PF) were described. Joint state-parameter estimation was used to estimate the model parameters in conjunction with the degradation state. The Paris case study was used as an illustration example. Synthetic data generated with the Paris' law were considered, and only model and measurement uncertainty were included in the study. The performance evaluation metrics showed the superiority of the PF over the EKF.

Chapter 3 considered a more realistic case study in order to validate the performances of the PF for prognosis. Indeed, fatigue crack growth prognosis under variable Conclusion amplitude loading (block loading) in composite materials was performed. Thus, in addition to model and measurement uncertainty, future loading uncertainty was also considered. The proposed methodology combined a PF with a detection algorithm, namely the two-sided cumulative sum (CUSUM). The proposed methodology was illustrated using real data from fatigue tests on composite materials realized at ONERA/DMAS. The performance evaluation metrics results showed the ability of the PF to handle highly nonlinear models and the different uncertainty sources.

Another framework for uncertainty representation and propagation was proposed in Chapter 4. This time, set-membership methods based on interval analysis were used. Only the bounds of the state and measurement noises are required without any stochastic assumption. This allows the estimation of the degradation state in a guaranteed way. Two different case studies were considered. In the first approach, a set-membership methodology using constraint satisfaction and setinversion techniques was adopted to estimate and then forecast the degradation state, based on the Waltz filtering algorithm and the Matlab R ⃝ toolbox VSIVIA. The method was applied to the Paris case study and compared to the results obtained with the PF and EKF in terms of accuracy, prediction and timeliness. The PF-based method remains the most efficient one, however, the interval method outperforms the EKF algorithm in terms of accuracy and timeliness. As for the precision, which measures the width of the interval that contains the predicted RUL, the results are lower than the EKF and the PF. Indeed, the interval computed with the setmembership method contains all the possible solutions, whereas the distributions obtained with the stochastic methods were truncated to the range (-3σ, 3σ).

In the second study of Chapter 4, the case of a hidden degradation state which cannot be directly measured was considered. A discrete-time linear time-invariant interval UIO was synthesized to jointly estimate the bounds of the system and of the degradation states in the context of bounded process and measurement noises. Then, the bounds of the future degradation state were predicted using the same interval propagation method as in the first study of Chapter 4. A suspension system subject to a crack growing in the suspension spring was used to illustrate the proposed method which has proven to be efficient.

Finally, the use of a reliability analysis method, namely the Inverse first-order reliability method (FORM) was investigated in Chapter 5. Inverse FORM allowed the calculation of the RUL cumulative distribution function (CDF) at desired points (median, probability bounds) at a low computational time and effort. The Inverse FORM was applied to the Paris case study and the performance evaluation metrics were computed. The Inverse FORM provided lower performance than the PF and the interval methods. However, it can be an alternative to PF-based methods and interval methods for RUL prognosis in the case where low computational cost algorithms are needed, and a high precision is not required. This thesis has focused on the development of model-based prognosis methods that take uncertainty into account and compute a measure of the RUL uncertainty. It was shown that uncertainty propagation can be tackled within various uncertainty representation frameworks. The choice of the method to use is mainly driven by the user requirements who can perform on-line RUL prognosis in the probabilistic, in the set-membership or in the reliability analysis framework. However, when highly nonlinear degradation models are considered, the PF-based methods have been found to be the most efficient ones for the application considered.

Perspectives

Based on the work presented in this thesis, some future investigations could be proposed to improve the presented model-based algorithms.

A key area for future work concerns the tuning of the filtering techniques. In this thesis, trial and error procedures were used to set the initial noise covariances. However, such method is not realistically feasible when data are not a priori available. Therefore, on-line estimation of the noise covariances would be a significant upgrade of the filtering method used in this work.

Although the particle filter has proved to be a powerful tool to handle model, measurement and future loading uncertainty, the various improvements (e.g. choice of proposal distribution) that can be found in the literature review of Chapter 1 can be applied to obtain more accurate state pdf estimations and thus more accurate RUL predictions.

In the case where future loading uncertainty is considered, more advanced detection algorithm such as the leave-out sign-dominant correlation regions (LSCR) algorithm [START_REF] Weyer | A randomised subsampling method for change detectio[END_REF]. With this algorithm, a region in parameter space which has a guaranteed probability of containing the true parameters can be identified.

The proposed interval UIO could be adapted to the case of nonlinear systems in order to extend its application domain as real systems are often nonlinear. [Gucik-Derigny, Outbib, Ouladsine, and Basualdo 2011b] has developed a continuous nonlinear interval observer for state and unknown input estimation. The authors suggested that it could be applied to prognosis but did not presented this part. However, the proposed observer is not adapted to the kind of nonlinear models presented in this thesis. Therefore, there is a need to develop an observer for state and unknown input estimation which is discrete-time, interval-based, and able to handle the type of nonlinearities involved in the models used for fatigue or structural monitoring in particular.

Stochastic and set-membership methods for prognosis were proposed and compared. The advantages of each kind of method have been highlighted. Thus, based on [START_REF] Abdallah | Box particle filtering for nonlinear state estimation using interval analysis[END_REF], where an extension of the particle filter Conclusion algorithm using interval analysis and constraint satisfaction techniques was proposed, a similar prognosis method could be developed.

With the emergence of technologies increasing the availability of data from the degrading systems, the proposed model-based prognosis methods could be combined with data-driven methods in order to obtain more accurate models and thus improve the RUL predictions.

Finally, some possible extensions of this work concern the interpretation and translation of the prognosis results into actions. Indeed, prognosis information could be very useful for maintenance decisions or control reconfiguration to ensure mission success and users safety. from (A.7): .8) where

z 1 (k) = R K [ỹ 2 (k) -Cy2 z 2 (k) -D2 u(k) -δ2 (k)] (A
R K = K -1 C R -1 C G ⊤ y H ⊤ C
. Finally, we have the following dynamical system:

⎧ ⎨ ⎩ z 2 (k + 1) = Ā2 z 2 (k) + B2 u(k) + Ē2 ỹ2 (k) -Ē2 δ2 (k) + ω2 (k) ỹ22 (k) = H ⊤ C2 Cy2 z 2 (k) + H ⊤ C2 D2 u(k) + H ⊤ C2 δ2 (k) (A.9)
where Ē2 = Ã21 R K , Ā2 = Ã22 -Ē2 Cy2 , B2 = B2 -Ē2 D2 ,

Introduction

Les systèmes critiques sont des systèmes complexes pour lesquels les anomalies de fonctionnement et la dégradation des composants peuvent avoir un impact considérable sur la sécurité des utilisateurs, la réussite des missions et les coûts liés à la maintenance. Les procédures de surveillance actuelles intègrent des outils de diagnostic pour la détection, la localisation et l'identification des défauts survenus sur le système. Selon leur nature, ces défauts altèrent lentement ou subitement le fonctionnement normal du système jusqu'à entraîner sa défaillance partielle ou totale (panne). Il est donc primordial de pouvoir suivre l'évolution de ces défauts et agir sur le système (reconfigurer les lois de commande, entreprendre une action de maintenance) avant sa défaillance.

De ce besoin est né il y a une dizaine d'années un domaine de recherche, le pronostic, qui consiste à calculer le temps qu'il reste avant que le système ou le composant où se trouve le défaut naissant ne soit défaillant. Ce temps est appelé durée de vie résiduelle du système (RU L : Remaining Useful Life). Ainsi, le pronostic permettra de développelr une stratégie de maintenance prédictive, qui consiste à surveiller et prédire l'état du système afin d'agir au plus près de son moment de défaillance. Cette maintenance prédictive permettra la disponibilité maximale du système et la minimisation des coûts qui peuvent être élevés dans le cas des maintenances classiques, tout en garantissant un système fiable et sécurisé.

Les méthodes de pronostic peuvent être classées en trois catégories principales [START_REF] Liu | A review on fault prognostics in integrated health management[END_REF]] : les approches à base de données, les approches à base de modèles et enfin les approches hybrides combinant les deux précédentes. Pour le pronostic à base de données, des données en ligne (température, pression, courant, ...) ainsi que des connaissances a priori sont utilisés pour extraire des caractéristiques permettant d'évaluer et de prédire l'état de santé du système. Pour cela, des méthodes basées sur l'intelligence artificielle et l'apprentissage statistique sont utilisées. Le pronostic B.2. Pronostic de propagation de fissure dans un matériau composite soumis à un chargement inconnu d'amplitude variable avoir une évaluation numérique plus précise. L'exactitude mesure la proximité de la RUL prédite par rapport à la vraie RUL. La précision mesure la largeur de l'intervalle contenant les RUL prédites. L'opportunité mesure la position relative de la RUL prédite par rapport à la vraie RUL, c'est-à-dire si c'est trop tôt ou trop tard. L'exactitude et la précision sont compris entre 0 et 1, où 1 reflète la meilleure performance. Plus l'opportunité est grande, plus on s'éloigne de la vraie distribution de la RUL, et les prédictions en retard ont un poids plus fort afin de les pénaliser. L'évaluation des performances a été réalisée sur 100 simulations. Les valeurs de ces métriques sont résumées dans le Tableau B.3.

B.2 Pronostic de propagation de fissure dans un matériau composite soumis à un chargement inconnu d'amplitude variable

Dans ce chapitre, le matériau dans lequel se propage la fissure est un matériau composite renforcé avec des fibres. Le chargement est cette fois inconnu et la structure est soumise à des variations brusques de chargement, par exemple dues à un choc. Un cas plus réaliste qui n'avait pas encore été traité dans la littérature a donc été considéré. De plus, des données réelles issues de tests de fatigue réalisés au département matériaux et structures de l'Onera (DMAS) ont été utilisées.

A Modèle de propagation de la fissure

Comme c'est un problème de pronostic qui n'avait pas encore été traité dans la litérature, il faut tout d'abord établir le modèle de propagation de la fissure.

Une loi de Paris modifiée a été utilisée : Ce modèle est fortement non linéaire et contient deux types de paramètres : (i) les paramètres de modèle constants qui sont m et C qui dépendent du matériau, et (ii) les paramètres de chargement ∆P et ∆σ qui dépendent du chargement appliqué à la structure et peuvent donc varier. Ainsi, dans un premier temps, un jeu de données a été utilisé pour estimer les paramètres constants m et C, puis d'autres jeux de données ont été considérés pour le pronostic de propagation de fissure avec les paramètres ∆P et ∆σ à estimer conjointement avec la longueur de la fissure.

da dN = C(∆K -K th ) m (B.

B Pronostic de propagation de fissure avec un chargement constant

Dans un premier temps, le chargement appliqué au matériau composite est constant. La En effet, le filtre particulaire diverge lorsqu'il y a une variation brusque de paramètres. On pourrait envisager dans un premier temps de laisser le filtre particulaire corriger ses estimations avec plus de données. Cependant ce n'est pas possible car on risque d'atteindre le seuil de dégradation. On pourrait également augmenter le nombre de particules, ce qui n'est pas possible car cela va considérablement augmenter le temps de calcul, ce qui n'est pas approprié au pronostic en temps réel. Finalement, un algorithme de détection a été ajouté afin de réinitialiser le filtre lorsque le chargement brusque a été détecté pour le faire converger vers les valeurs réelles des paramètres. (B.17) où µ 0 est la valeur initiale de la moyenne de ∆x k et ν correspond à l'amplitude du changement de moyenne à détecter.

Ainis, lorsqu'un changement est détecté, les valeurs des paramètres estimés est réinitialisée pour que le filtre particulaire converge vers les vraies valeurs après la variation. Les valeurs de réinitialisation sont déterminées à partir de la connaissance de la valeur de chargement critique ∆σ crit au-delà de laquelle la longueur de la fissure risque d'augmenter plus rapidement. Dans ce cas, ∆σ crit = 450 MPa. Afin de tenir compte de l'incertitude sur cette valeur, ∆σ a et ∆P après la variation brusque sont réinitialisés dans un intervalle L = [400, 500] MPa contenant ∆σ crit .

Les résultats d'estimation et de prédiction de la longueur de la fissure après la variation brusque de chargement avec le filtre particulaire combiné à l'algorithme de détection two-sided CUSUM sont représentés sur la Figure B.7.

Un intervalle contenant les paramètres réinitialisés après la variation de chargement a été considéré, ainsi, un intervalle contenant la longueur de fissure en utilisant les valeurs extrêmes des paramètres a ou être obtenu. Plus il y a de données, plus les estimations et prédictions se rapprochent de la longueur de fissure de référence.

Ainsi, le filtre particulaire combiné à l'algorithme two-sided CUSUM a permis de traiter le cas où les chargement extérieurs sont inconnus et peuvent varier brusquement.

B.3 Conclusion

Dans cette thèse, différentes méthodes de pronostic à base de modèle ont été développées. Le principal enjeu concernait l'intégration et le traitement des différentes sources d'incertitudes afin d'évaluer l'incertitude autour de la RUL prédite. Il y a principalement les incertitudes de modèle, les incertitudes de mesures et les incertitudes liées aux futures conditions d'opération du système. Afin de gérer ces incertitudes, des méthodes de pronostic probabilistes et ensemblistes ont été développées. Afin d'illustrer le fonctionnement des méthodes, le pronostic de propagation de fissure a été considéré. Les performances de ces méthodes ont ensuite été évaluées et comparées à travers trois métriques : l'exactitude, la précision et l'opportunité.

Dans un premier temps, un filtre de Kalman étendu ainsi qu'un filtre particulaire sont appliqués au pronostic de propagation de fissure, en utilisant la loi de Paris et des données synthétiques. Puis, une méthode combinant un filtre particulaire et un algorithme de détection (algorithme des sommes cumulatives) a été développée

B.3. Conclusion

puis appliquée au pronostic de propagation de fissure dans un matériau composite soumis à un chargement variable. Cette fois, en plus des incertitudes de modèle et de mesures, les incertitudes liées aux futures conditions d'opération du système ont aussi été considérées. De plus, des données réelles ont été utilisées. Ensuite, deux méthodes de pronostic sont développées dans un cadre ensembliste où les erreurs sont considérées comme étant bornées. Elles utilisent notamment des méthodes d'inversion ensembliste et un observateur par intervalles pour des systèmes linéaires à temps discret. Enfin, l'application d'une méthode issue du domaine de l'analyse de fiabilité des systèmes au pronostic à base de modèles est présentée. Il s'agit de la méthode Inverse First-Order Reliability Method (Inverse FORM).

Les résultats obtenus ont permis de montrer que la propagation d'incertitude peut être effectuée avec différentes méthodes de représentation d'incertitude. Le choix de la méthode à utiliser dépend principalement du domaine d'étude de l'utilisateur qui peut effectuer le pronostic avec des méthodes probabilistes à base de filtres stochastiques, des méthodes ensemblistes basées sur l'analyse par intervalles ou encore des méthodes issues du domaine de la fiabilité des systèmes. Cependant, lorsque des modèles de dégradation fortement non-linéaires sont en jeu, les méthodes à base de filtre particulaire restent les mieux adaptées pour les cas qui ont été traités dans cette thèse.

Filtering and uncertainty propagation techniques for model-based prognosis

Elinirina Iréna ROBINSON

  Figure 1.1: Prognosis process scheme

  represents the interval around the ground-truth RUL (Fig. 1.2). The values of this timeliness function are in the interval [0; +∞] and the perfect score for timeliness is 0.

Figure

  Figure 1.2: Timeliness metric

  perform the identification of the model parameters in conjunction with state estimation.

  Figure 2.1: Principle of the SIR particle filter

Figure 2

 2 Figure 2.2: Results obtained with the EKF at N p = 1200 cycles

  Figure 2.3: 1000 particle trajectories for N p = 1200 cycles

Figure 2

 2 Figure 2.4: Results obtained with the PF with N part = 100 particles and N p = 1200 cycles

Figure 2

 2 Figure 2.5: Results obtained with the PF with N part = 500 particles and N p = cycles

Figure 2

 2 Figure 2.6: Results obtained with the PF with N part = 1000 particles and N p = cycles

Figure 2

 2 Figure 2.8: Results of 100 experiments obtained with an EKF

Figure

  Figure 2.9: Results of 100 experiments obtained with a PF

  Figure 3.1: Schematic of the specimen showing the dimensions and the loading axis parallel to fiber direction

Figure 3 . 2 :

 32 Figure 3.2: Photography of the specimen with a crack perpendicular to fiber direction

Figure 3 Figure 3 . 5 :Figure 3 . 6 :

 33536 Figure 3.4: Crack growth prognosis with the EKF at N p = 1.198 × 10 5 cycles

Figure 3 Figure 3

 33 Figure 3.10: Crack length evolution at N p = 7.702 × 10 6 cycles

Figure 3

 3 Figure 3.12: Evolution of the crack length at different prediction cycles N p without the detection algorithm

Figure 3 Figure 3

 33 Figure 3.15: Evolution of the crack length at different prediction cycles N p with the detection algorithm

Figure 4 Figure 4

 44 Figure 4.1: Estimation of the feasible domain of m and log(C) with VSIVIA

Figure 4

 4 Figure 4.3: Crack size bounds prediction and estimation

Figure 4 . 4 :

 44 Figure 4.4: Quarter-car suspension model

  Figure 4.5: Road surface excitation

  Figure 4.6: State estimation (vertical displacements x 1 and x 2 )

Figure 4

 4 Figure 4.7: Estimation of the bounds of the unknown input d

Figure 4 . 8 :

 48 Figure 4.8: Estimation of the bounds of the stiffness degradation k 1

Figure 4 . 9 :

 49 Figure 4.9: Zoom on k 1

Figure 4

 4 Figure 4.10: Estimation of the bounds of the degradation ξ

  59) where ∆N = 50 cycles and ∆σ = 78 MPa as in Section 4.1. The first step consists in estimating the bounds of the current degradation state. The feasible domain of the parameter vector θ = [m C] is estimated based on the bounds of the degradation signal computed with the measurements recorded until the prediction time k p or equivalently time cycle N p . It can be considered that m ∈ [2, 5] and log(C) ∈ [-27 -19] as the spring is a metallic component.

Figure 4 .

 4 Figure 4.11: Estimation of ξ bounds with the measurements

Figure 4 .Figure 4

 44 Figure 4.12: Degradation ξ bounds estimation and prediction at time cycle N p = 0, 20 × 10 6 cycles

Figure 4 .

 4 Figure 4.14: Degradation ξ bounds estimation and prediction at time cycle N p = 2, 99 × 10 6

Figure 5

 5 Figure 5.1: CDF and Inverse CDF

Figure 5 . 2 :

 52 Figure 5.2: Computation of a RUL CDF value with both FORM and Inverse FORM

  5 or to choose several values of P f to reconstruct the entire RUL CDF. Examples of the RUL CDF calculated at N p = 1200 cycles obtained with different values of P f are depicted from Fig. 5.3 to Fig. 5.5.

Figure 5 Figure 5

 55 Figure 5.3: RUL CDF generated with values of P f from 0.001 to 0.99 every 0.001 steps

Figure 5

 5 Figure 5.7: Results of 100 experiments obtained with Inverse FORM

  Figure B.5 représente les estimations et prédictions de la longueur de la fissure et des paramètres ∆P et ∆σ.

B. 2 .

 2 Figure B.5: Estimation et prédiction de la longueur de fissure et des paramètres de chargement

  Pour la détection, l'algorithme CUSUM (CUmulative SUm) a été utilisé. C'est un algorithme très utilisé pour la détection de changement de moyenne d'une distribution gaussienne. Afin de détecter le changement que ce soit dans le cas d'une diminution ou d'une augmentation, un algorithme CUSUM bilatéral (two-sided CUSUM) est considéré. L'idée générale est de calculer le temps d'arrêt N detect de l'algorithme qui correspond au moment où le changement de moyenne du paramètre ∆x k surveillé est 160B.2. Pronostic de propagation de fissure dans un matériau composite soumis à un chargement inconnu d'amplitude variabledétecté :N detect = min{k : S + k ≥ S th ∪ S - k ≥ S th } (B.16) où S + k et S - k sont respectivement les fonctions de décision correspondant aux tests d'augmentation et de diminution de la moyenne de ∆x k , et S th est un seuil de détection. Les expressions des sommes cumulatives S + k et S - k sont les suivantes : - k-1 -∆x kµ 0 -ν 2 )

  Figure B.7: Estimation et prédiction de la longueur de fissure et des paramètres de chargement

  Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Recursive Bayesian filtering . . . . . . . . . . . . . . . . . . . . . . . Schematic of the specimen showing the dimensions and the loading axis parallel to fiber direction . . . . . . . . . . . . . . . . . . . . . . 3.2 Photography of the specimen with a crack perpendicular to fiber direction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 The closure pressure in the crack wake . . . . . . . . . . . . . . . . .

	2 Filtering methods for prognosis with stochastic assumptions	31
	2.1

nos actions, et nous dirigeons le reste.» Nicolas Machiavel iii Introduction 3 1 Prognosis and uncertainty propagation 9 1.1 Prognosis approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.1.1 Data-driven approaches . . . . . . . . . . . . . . . . . . . . . 1.1.2 Model-based approaches . . . . . . . . . . . . . . . . . . . . . 1.1.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Model-based prognosis . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Principle of model-based prognosis . . . . . . . . . . . . . . . 1.2.2 Uncertainty in model-based prognosis . . . . . . . . . . . . . . 1.3 Performance evaluation metrics . . . . . . . . . . . . . . . . . . . . . 1.4 2.1.1 Stochastic methods . . . . . . . . . . . . . . . . . . . . . . . . 2.1.2 Recursive Bayesian estimation principle . . . . . . . . . . . . . 2.2 Extended Kalman filter . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures 1.1 Prognosis process scheme . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Timeliness metric . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Principle of the SIR particle filter . . . . . . . . . . . . . . . . . . . . 2.2 Results obtained with the EKF at N p = 1200 cycles . . . . . . . . . . 2.3 1000 particle trajectories for N p = 1200 cycles . . . . . . . . . . . . . 2.4 Results obtained with the PF with N part = 100 particles and N p = 1200 cycles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Results obtained with the PF with N part = 500 particles and N p = 1200 cycles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Results obtained with the PF with N part = 1000 particles and N p = 1200 cycles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 RUL PDF obtained with the PF with N part = 1000 particles and N p = 1200 cycles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Results of 100 experiments obtained with an EKF . . . . . . . . . . 2.9 Results of 100 experiments obtained with a PF . . . . . . . . . . . . 2.10 EKF RUL PDF for 100 experiments . . . . . . . . . . . . . . . . . . . 2.11 PF RUL PDF for 100 experiments . . . . . . . . . . . . . . . . . . . 3.1 3.4 Crack growth prognosis with the EKF at N p = 1.198 × 10 5 cycles . . 3.5 Crack growth prognosis with the EKF at N p = 2.794 × 10 5 cycles . . 3.6 Crack growth prognosis with the EKF at N p = 3.542 × 10 5 cycles . . xi

Table 2

 2

	.1: True model parameters to generate the synthetic data
	∆σ ∆N a true log(C true ) m true
	78	50 0.01	-22.62	3.8

Table 2 .

 2 2: Initial distribution of the random variables

	Variance σ 2 log(C)

  Table 2.3.

				EKF 100 runs for N p =1200 cycles
		0.05
	Crack size (m)	0.02 0.03 0.04
						Reference crack size
						Measurements
		0.01				Estimated crack size
						Predicted crack size
						Threshold
		0 0	500	1000	1500	2000	2500	3000
					N (cycle number)

Table 2

 2

		.3: Performance evaluation results
	Method Accuracy Precision Timeliness
	EKF	0.7150	0.7630	1.0084
	PF	0.9511	0.7030	0.1132

Table 3 .

 3 1: Crack growth trajectories dataset provided by ONERA/DMAS

	Specimen Notch length (mm) Number of data Loading ∆σ a (MPa)
	1	1,175	228	Constant	326,
				Block 1	194, 4
	2	2,88	5910	Block 2 Block 3	243 303,
				Block 4	379,
	3	0,75	2314	Block 1 Block 2	360 450
	4	2, 86	584	Constant	196, 2
	5	0, 73	386	Constant	278, 1
	6	0, 76	565	Constant	221, 4
				Block 1	199, 8
				Block 2	249,
				Block 3	312,
	7	0,75	1740	Block 4	389,
				Block 5	487,
				Block 6	560,
				Block 7

Table 3 .

 3 2: Distributions of random parameters for the estimation of m and C

	Parameter log(C)	m	w	v
	Variance

Table 3

 3

	Prediction time Method Accuracy Precision Timeliness
	1.198 × 10 5	PF EKF	0.88 0.67	0.72 0.93	0.77 4.83
	2.370 × 10 5	PF EKF	0.95 0.71	0.82 0.95	0.17 2.04
	3.542 × 10 5	PF EKF	0.96 0.72	0.85 0.95	0.08 0.91

.3: Performance evaluation results for the estimation of m and C

Table 3 .

 3 5: Performance evaluation results for the estimation of ∆σ a and ∆P

	Prediction time	Accuracy	Precision	Timeliness
	11.058 × 10 6 11.197 × 10 6 11.335 × 10 6 12.337 × 10 6	0.59 (±0.04) 0.62 (±0.02) 0.65 (±0.04) 0.75 (±0.02) 0.81 (±0.08) 0.86 (±0.04) 0.87 (±0.15) 0.59 (±0.07) 207.12 (±486.04) 4.56 (±1.05) 2.74 (±0.73) 0.92 (±0.46)

Table 4

 4

		.2: Performance evaluation results
	Method Accuracy Precision Timeliness
	Interval	0.9334	0.4493	0.1536
	EKF	0.7150	0.7630	1.0084
	PF	0.9511	0.7030	0.1132

  .3.

	Parameter	Value Unit
	Sprung mass (m 1 )	315	kg
	Unsprung mass (m 2 )	51	kg
	Suspension stiffness (k 1 )	43.3	kN/m
	Tire stiffness (k 2 )	210	kN/m
	Suspension damping coefficient (b 1 ) 3.9	kN•s/m
	Tire damping coefficient (b 2 )	1.1	kN•s/m
	Table 4.3: Values of the model parameters
	The input of the system u =

[ ẋr x r ] ⊤

  k N 1 = 43.3 × 10 3 N/m in this work.Here, the bounds of ξ are computed at each time step k with VSIVIA. The inputs 1024.2. Unknown input interval observer for prognosis of discrete-time linear systems

			x 10 4		Degradation of the stifness k 1
		4.3
		4.28
	k 1 (N/m)	4.24 4.26
		4.22
			Upper bound
			Lower bound
		4.2	Reference
			0.6	0.8	1	1.2	1.4	1.6
					Time (cycles)	x 10 6

  ).

				Degradation level estimation with the measurements
		0.25
			Measurements
			Estimated Lower bound
			Estimated Upper bound
		0.2
	ξ	0.15
		0.1
		0 0.05	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2
						Time (cycles)				x 10 6

Table 4

 4

		.4: Calculated RUL intervals
	Prediction time RUL lower bound k RU L True RUL RUL upper bound k RU L
	0, 20 × 10 6 1, 90 × 10 6 2, 99 × 10 6	2, 29 × 10 6 1, 56 × 10 6 1, 29 × 10 6	4, 22 × 10 6 2, 52 × 10 6 1, 43 × 10 6	5, 15 × 10 6 3, 25 × 10 6 1, 59 × 10 6

Table 4

 4

	.5: Performance evaluation results
	Prediction time Accuracy Precision Timeliness
	0, 20 × 10 6 1, 90 × 10 6 2, 99 × 10 6	0.8891 0.9544 0.9936	0.0572 0.1847 0.7467	1.6950 0.2657 0.0185

Table 5

 5

		.1: Performance evaluation results
	Method Accuracy Precision Timeliness
	IFORM	0.8989	0.3251	0.2716
	PF	0.9511	0.7030	0.1132
	Interval	0.9334	0.4493	0.1536
	EKF	0.7150	0.7630	1.0084
	The results show that the Inverse FORM outperforms the EKF for RUL com-
				123

Reliability method for prognosis

Table B

 B Le filtre particulaire donne de meilleurs résultats que l'EKF en terme de précision et d'exactitude. Ceci peut être principalement dû aux approximations linéaires effectuées par l'EKF.

		.3: Evaluation des performances
	Method Accuracy Precision Timeliness
	EKF	0.7150	0.7630	1.0084
	PF	0.9511	0.7030	0.1132

  12)où K th est la valeur de ∆K à partir de laquelle la fissure commence à se propager. L'expression de ∆K fait intervenir un terme qui permet de modéliser la tension exercée par les fibres : La principale difficulté dans ce modèle concerne l'estimation de ∆P qui se fait généralement par la méthode des éléments finis. Cependant ce n'est pas adapté au pronostic en ligne à cause des temps de calcul beaucoup trop importants. Dans la littérature, il a été proposé de considérer delta P comme étant uniforme. Cette hypothèse a été retenue pour simplifier le modèle et l'adapter au pronostic.Finalement, le modèle de propagation de fissure utilisé dans ces travaux est le suivant :

				∆K = F ∆σ a	√	πa -2F	√ a π	∫ a a 0	∆P (x) √ a 2 -x 2 dx	(B.13)
	où ∆P représente la force exercée par les fibres.
	Ainsi, le modèle de propagation de fissure est :
		da dN	= C	( F ∆σ a	√	πa -2F	√ a π	∫ a a 0	∆P (x) √ a 2 -x 2 dx -K th	) m	(B.14)
	da dN	= C	( F ∆σ a	√	πa -2F	√ a π	∆P	( π 2	-arcsin	( a 0 a	))	-K th	) m	(B.15)

1.1. Prognosis approaches

Durant les essais, le chargement appliqué au matériau composite est augmenté brusquement, ce qui a pour effet d'accélérer la propagation de la fissure. Cependant, le filtre particulaire n'arrive pas à estimer et prédire la longueur de la fissure après la variation brusque même avec plus de données(Figure B.6).

Remerciements

Estimation of the current degradation state using available measurements

In this case, the degradation state is the crack length a and it is directly measured. The crack size measurements from initial time cycle 0 to prediction time cycle N p = 1200 are used to estimate the current degradation state. First, the feasible domain Θ of the parameter vector θ =

[ log(C) m

] is estimated using a set inversion technique. Then, the bounds of the current crack size are computed with the natural inclusion function and the estimated set Θ. For this purpose, VSIVIA is used.

Estimation of Θ using VSIVIA. To estimate Θ, an initial box containing the parameter set should be specified first. The values of the parameters m and C may differ from one specimen to another due to material and operating conditions variability. However, they have known range values that depend on the considered material. Here, aluminium plates are considered, therefore, according to [START_REF] Steadman | On the form of fatigue crack growth formulae[END_REF], m ∈ [2, 5] and log(C) ∈ [-27, -19].

Thus, the set of input parameters of the VSIVIA algorithm is:

The obtained set is shown in Figure 4.1 where Θ int appears in red, Θ ext in white and the part that do not belong to the solution set is in blue.

Estimation of the bounds of the current crack length using VSIVIA.

Then, using the estimated Θ, the bounds of the current crack length are computed. The estimated bounds are shown in Figure 4.2, and it can be seen that the measurements are all contained in these bounds.

Estimation of the future degradation state

As it was explained in Section 4.1.2, the natural inclusion function of the state equation is used to compute the prediction of the degradation bounds:

where [m] and [log(C)] are derived from the projection of Θ out on the axes m and log(C).

Figure 4.3 shows the estimated bounds from time cycle 0 to prediction time cycle 1200, and the predicted bounds until the threshold is reached.

In both steps (estimation and prediction), the estimated bounds include the actual reference crack size. Indeed, the interval analysis techniques allow a guaranteed estimation.

The results of the estimated bounds of the RUL, as well as the true RUL value are depicted in Table 4.1. Performance evaluation metrics were used to evaluate the efficiency of the proposed algorithm and to compare it to the results obtained in Chapter 2 with the EKF and the PF.

Performance evaluation metrics

The same metrics (accuracy, precision and timeliness) are calculated here, using the same parameter values as in Chapter 2, R 0 = R min = R max = 700 cycles. The metrics values that are obtained are listed in Table 4.2 and compared to those obtained previously with the PF and the EKF.

As expected by looking at Figure 4.3, the precision of the interval-based method is lower than the PF and EKF. However, the estimated interval is guaranteed to contain all the uncertainty associated to the predicted RUL, while the precision metrics obtained with the stochastic filters are based on a truncated Gaussian distribution. Here, the distributions were truncated to the range (-3σ, 3σ), and with other values the precision would have not been the same.

Nevertheless, in terms of accuracy and timeliness, the interval method can be compared to the PF method, and largely outperforms the EKF.

The computational time of the interval method remains reasonable as long as the precision parameter ϵ is not chosen too small. Thus, based on the required precision,

Unknown input interval observer for prognosis of discrete-time linear systems

By replacing this expression of z 1 (k) in the second equation of (4.26) we obtain:

Finally we have the following dynamical system:

where

and F 2 = Ã21 E. In order to be able to design an interval observer for the discrete-time system (4.31), the following assumption which is standard in the field of observer design is required [START_REF] Hou | Design of observers for linear systems with unknown inputs[END_REF].

Based on Assumption 2 the following lemma allows to transform (4.31) into a suitable form for interval observer design [Efimov et al. 2013a].

Lemma 3 There exists a gain L ∈ R (n-q)×(p-q) and a transformation matrix P of appropriate dimensions such that (A 2 -LC 2 ) is Schur stable and R = P (A 2 -LC 2 )P -1 is nonnegative.

Such a transformation always exists, and in the case where the eigenvalues of (A 2 -LC 2 ) are real, R can be chosen as diagonal or as Jordan form of A 2 -LC 2 [Efimov et al. 2013a]. After the change of coordinates r 2 = P z 2 , the system (4.31) is described in the new coordinates by:

where M = P (F 2 + LH ⊤ 012 ).

Interval observer design for state and unknown input estimation

In this section, an interval observer is designed for state estimation and then for unknown input estimation.

From reliability analysis to prognosis

Basics of reliability analysis

During its lifetime, a system is subject to various solicitations that can degrade its performance and in time, can cause its failure. In order to ensure the safety, serviceability and durability of the system, reliability analysis was introduced to evaluate the probability that the system does not behave as desired during a given period of time, whilst taking uncertainties into account. This probability is referred to as probability of failure or failure probability of the system. The basic principle of reliability analysis which is required for the understanding of this chapter is briefly described in what follows, but the reader may refer to [START_REF] Sudret | Uncertainty propagation and sensitivity analysis in mechanical models-Contributions to structural reliability and stochastic spectral methods[END_REF][START_REF] Lemaire | Structural reliability[END_REF] to have more details.

Uncertainty propagation in reliability analysis. Reliability analysis aims at studying the impact of the input uncertainties (material properties, geometry, loading conditions, etc.) on the behavior of a system with respect to a certain failure scenario.

The input uncertainties are usually described within a probabilistic framework by a s-dimensional continuous random vector X containing the input variables affected by uncertainties.

The system is modeled by a function M : R s → R, X ↦ → Y = M(X) whose output is supposed to be scalar. This function can be either analytical or a black-box computer code (e.g., a finite element code) which can be only known pointwise. Due to the presence of uncertainties affecting the input variables, the output Y ∈ R is a random variable too. Thus, the goal is to compute the pdf of Y through uncertainty propagation from the input to the output of the model.

Finally, for reliability assessment, the quantity of interest is the relative behavior of the output Y with respect to a failure criterion (e.g., a threshold overrun). To characterize this failure criterion, a so-called limit-state function G is introduced:

where y th ∈ R is a safety threshold value characterizing the behavior of the system. If the realization y of Y equals or exceeds the threshold, the system is considered in a failure state. In this way, if x denotes the realizations of the input vector X, this limit-state function G separates the variable space in two domains:

• D s = {x : G(x) > 0} defines the safety domain;

this time instant, the forecasting of the degradation state in the future is carried out for time instants k > k p without new measurements. The main challenge in this prediction step lies in the fact that future operational conditions of the system are unknown. Therefore, the forecasting of the degradation state must be performed by taking uncertainties into account.

The future state is predicted until the failure threshold is reached, giving the predicted failure time k f . Finally, the RUL at time k p which is denoted by R(k p ) can be calculated as:

(5.5)

For better readability, the dependency in k p will be omitted in the following and the RUL simply denoted by R.

However, as the future state prediction is uncertain, the predicted failure time k f is uncertain, making the RUL a random variable that depends on:

If the vector X contains all of the uncertain quantities mentioned earlier i.e X = [x kp , u kp , u kp+1 , . . . ] ∈ R s where s is the number of uncertain parameters, and M is a function M : R s → R, X ↦ → M(X) that allows to compute the RUL, then we have: R = M(X).

(5.6)

The model M can be a known function defined by a mathematical expression of the RUL, or may be a black-box function that takes input values and provides a result.

Thus, RUL calculation can be considered as an uncertainty propagation problem. Indeed, computing the RUL pdf f R is equivalent to propagating the uncertainty in X through M.

Therefore, uncertainty propagation methods issued from the field of reliability analysis can be adapted to the problem of RUL computation.

Chapter 5. Reliability method for prognosis

Inverse FORM for RUL calculation

In this section, a methodology to compute the RUL with the Inverse FORM is presented.

The methodology can be applied if any type of solution to calculate the RUL from the dynamical degradation model is found. The solution may be analytical, and in the case where no analytical expression is available, a numerical or an approximate solution can be used.

Therefore, the first step to perform RUL prediction with the Inverse FORM is to establish a computational model M whose output is a RUL value. For example, in [Bressel, Hilairet, Hissel, and Bouamama 2016a], RUL estimation of a Proton Exchange Membrane Fuel Cell (PEMFC) has been performed. A simple mathematical expression of the RUL in function of the state of health (SoH) of the PEMFC and its degradation speed was obtained using experimental data. In their work, the random state variable X contained the SoH and the degradation speed variables. In [Sankararaman, Daigle, and Goebel 2014a], the RUL was computed using a computational code, i.e an implicit model M. Indeed, the random state variable X contained the future loading (considered constant but random) and especially the process noise generated from the prediction time step k p to the time k f at which the degradation threshold is reached. Thus, the size of X depends on k f , making the RUL computation model implicit. Here, the RUL computational model is derived from a physics-based model of the degradation.

Once this model is available, the Inverse FORM algorithm is applied to find the input parameters vector X for a given failure probability

For the numerical search of the uncertain parameter vector X, the following constraints must be satisfied:

where β is the reliability index as defined in Section 5.2.1, u is a realization of U which is the random vector X expressed in the standard normal space and ||•|| is the Euclidean norm. To find the optimum solution to (5.11.c), a numerical search is required (see Der [START_REF] Der Kiureghian | Inverse reliability problem[END_REF] for more details). The steps that are followed for each iteration j of the algorithm are:

1. Setting j to 0 and initial guess of the realizations of the uncertain parameter vector x j = [x j 0 , . . . , x j s ] where s is the number of uncertain parameters.

2. Transformation into the standard normal space. In the case of Gaussian

Appendices

Appendix of Chapter 4 Extension of the interval UIO

Consider the following LTI discrete-time system:

Step 1: Unknown input-free state vector

As it was done in Section 4.2.1, the state vector in eq. (A.1) is decomposed into an unknown input depending subsystem and an unknown input-free subsystem using the transformation x = Hz:

Step 2: Unknown input-free measurement vector

Similarly to the procedure applied in Section 4.2.1 to decompose the state vector, the measurement vector is divided into an unknown input depending subsystem and an unknown input-free subsystem.

Assumption 3 C is a full row rank and G is a full column rank matrix.

Under Assumption 3, there exists an orthogonal matrix

and matrices R 0y ∈ R q×q and K y ∈ R q×q such that:

This leads to the transformation of the measurement equation of system (A.1) into an equivalent one: .4) where:

Finally, with x = Hz, we obtain the following unknown input-free subsystem:

where Cy1 ∈ R (p-q)×q and Cy2 ∈ R (p-q)×(n-q) .

Step 3: Eliminate z 1 from system (A.5)

Cy1 is supposed to be a full column rank matrix and can be decomposed as:

×q and H C2 ∈ R q×q ) and ỹ2 = H ⊤ C ỹ2 , we can decompose ỹ2 (k) as:

, the expression of z 1 can be extracted à base de modèles fait appel à un modèle physique qui décrit l'évolution de la dégradation du système. Enfin, les méthodes hybrides combinent à la fois les données disponibles sur le système et le modèle physique de la dégradation. Les travaux de cette thèse s'inscrivent dans le cadre des méthodes de pronostic à base de modèles. Cette catégorie de méthodes peut s'avérer difficile à mettre en place car un modèle physique de la dégradation n'est pas toujours disponible. Cependant, lorsque c'est le cas, les méthodes de pronostic à base de modèles sont très précises car elles intègrent non seulement la connaissance des conditions d'opération du système mais aussi l'effet de leur variabilité sur son comportement.

Le pronostic à base de modèle est composé de deux étapes principales (Figure B.1) : (i) estimation de l'état actuel de la dégradation et (ii) prédiction de l'état futur de la dégradation. La première étape, qui est une étape de filtrage, est réalisée à partir du modèle et des mesures disponibles. La seconde étape consiste à faire de la propagation d'incertitudes. Le principal enjeu du pronostic concerne la prise en compte des différentes sources d'incertitude pour obtenir une mesure de l'incertitude associée à la RUL prédite [START_REF] Sankararaman | Significance , interpretation , and quantification of uncertainty in prognostics and remaining useful life prediction[END_REF]]. Les principales sources d'incertitude sont les incertitudes de modèle, les incertitudes de mesures et les incertitudes liées aux futures conditions d'opération du système. Afin de gérer ces incertitudes et de les intégrer au pronostic, il est nécessaire de choisir des méthodes de représentation et de propagation d'incertitudes appropriées. Parmis les approches existantes, il y a par exemple celles basées sur l'analyse par intervalles, la théorie des ensembles flous, ou encore les approches probabilistes. In the literature, modeling and measurement uncertainty are mainly handled with probabilistic methods. The random variables are represented using probability distributions. Particle filtering-based techniques are commonly used as they can deal with nonlinearities. However, Monte-Carlo simulation-based approaches require huge computational efforts and are not therefore practicable for online RUL prediction.

As regards future loading conditions uncertainty, they are not always taken into account. When it is the case, high computational time methods such as finite element models are involved. Moreover, in order to reduce computational time, the models are often linearized or simplified through the use of a surrogate model. However, these procedures tend to increase modeling uncertainty.

B.1 Filtre de Kalman étendu et filtre particulaire pour le pronostic

Ce chapitre a pour objectif de présenter plus en détail le principe et les principaux enjeux du pronostic à base de modèles à travers un exemple d'application. Le pronostic à base de modèles se compose de deux grandes étapes : (i) l'estimation de l'état de degradation actuel et (ii) la prédiction de l'état futur de la dégradation.

Cependant, ces deux étapes sont entachées d'incertitudes. Ainsi, le principal enjeu du pronostic est de tenir compte des incertitudes de modèle, de mesures et sur les futures conditions d'opération du système. Dans ce chapitre, seules les incertitudes de modèle et de mesures sont considérées.

Afin de traiter ces incertitudes, deux méthodes basées sur le filtrage stochastique sont utilisées : le filtre de Kalman étendu et le filtre particulaire. Les méthodes proposées sont appliquées au pronostic de propagation de fissure utilisant la loi de Paris.

A Filtre de Kalman étendu et filtre particulaire

Le filtrage stochastique a pour objectif d'estimer la densité de probabilité p(x k |y 0:k ) qui donne des informations statistiques sur l'état de dégradation x k , basé sur les mesures y 0:k = Y k = {y 0 , • • • , y k }. L'estimation conjointe de l'état et des paramètres du modèle est effectuée en formant un vecteur d'état augmenté X = Le filtre de Kalman étendu (EKF) est une adaptation du filtre de Kalman classique aux systèmes non linéaires. La densité de probabilité (pdf) de l'état de dégradation p(X k |Y k ) est approximée par une distribution Gaussienne.

En ce qui concerne le filtre particulaire (PF), la pdf de l'état de dégradation à l'instant k est approchée par un ensemble de N part particules qui ont chacune un poids associé

où δ est la fonction de Dirac.

Il existe de nombreux algorithmes de PF [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF], et celui qui a été considéré ici est le filtre à rééchantillonnage par importance (SIR : sequential importance resampling).

A.1 Estimation de l'état actuel de dégradation avec l'EKF

Pour l'étape de filtrage avec le filtre de Kalman étendu, on retrouve les équations d'un filtre de Kalman étendu classique avec les étapes de linéarisation, prédiction et correction à partir des mesures relevées. A.2 Prédiction de l'état futur de dégradation avec l'EKF Pour l'estimation de l'état futur de la dégradation, plus aucune mesure n'est disponible. Ainsi pour l'étape de correction, le bruit d'état est généré à partir de la matrice de covariance à l'instant k puis est injecté dans l'équation du modèle. Cette étape est réalisée jusqu'à atteindre le seuil de dégradation.

Prédiction

Prédiction

Xk|k-1 = f ( Xk-1|k-1 , w k = 0) (B.6)

A.3 Estimation de l'état actuel de dégradation avec le PF

Pour l'étape de filtrage avec le filtre particulaire, un algorithme de filtre particulaire classique avec une étape de ré-échantillonnage a été utilisé, le poids des particules est mis à jour à chaque nouvelle mesure :

A.4 Prédiction de l'état futur de dégradation avec le PF

Puis pour l'étape de propagation d'incertitudes, comme il n'y a plus de mesures disponibles, il a été considéré que le poids des particules était invariant à partir de l'instant de prédiction jusqu'à atteindre le seuil de dégradation.

Cette approximation peut être considérée car l'erreur induite est négligeable [START_REF] Orchard | A particle-filtering approach for on-line fault diagnosis and failure prognosis[END_REF].

B Application au pronostic de propagation de fissure avec la loi de Paris

Dans un premier temps, ces deux méthodes ont été appliquées au pronostic de propagation de fissure en utilisant la loi de Paris.

B.1 Modèle de propagation de fissure

L'objectif était de calculer le temps de vie restant d'une structure dans laquelle une fissure se propage, c'est-à-dire calculer le temps qu'il reste avant que la fissure n'atteigne une longueur critiques. La loi de Paris-Erdogan qui est très utilisée dans la littérature pour le pronostic de propagation de fissure a été utilisée :

Ce modèle décrit l'évolution de la longueur de fissure a en fonction du nombre de cycle N . Le chargement appliqué δσ est connu et vaut 78 MPa et les paramètres du modèle m et C sont inconnus.

B.2 Paramètres de la simulation

Les étapes sont les mêmes pour l'EKF et pour le PF.

Discrétisation du modèle Le modèle est discrétisé de la manière suivante :

Comme le paramètre C a une valeur très faible, log(C) est considéré ici.

Le vecteur d'état augmenté permet l'estimation conjointe de la longueur de fissure et des paramètres m et C :

Génération des données de mesure synthétiques Le modèle a été utilisé pour générer la longueur de fissure de référence ainsi que les mesures en y ajoutant un bruit. Les valeurs suivantes ont été utilisées :

Un nombre total de 1000 particules a été considéré pour le PF.

B.3 Résultats numériques

L'état actuel de la dégradation est estimé en utilisant les mesures jusqu'au cycle de prédiction N p = 1200.

Les résultats obtenus avec l'EKF sont représentés sur la Figure B.2. L'étape de filtrage permet d'estimer la longueur de la fissure (en bleu) en utilisant les mesures et l'étape de propagation d'incertitudes permet de prédire la longueur de la fissure (en vert) dans le futur jusqu'au seuil de dégradation.

Abstract

In this manuscript, contributions to the development of methods for on-line model-based prognosis are presented. Model-based prognosis aims at predicting the time before the monitored system reaches a failure state, using a physicsbased model of the degradation. This time before failure is called the remaining useful life (RUL) of the system. Model-based prognosis is divided in two main steps: (i) current degradation state estimation and (ii) future degradation state prediction to predict the RUL. The first step, which consists in estimating the current degradation state using the measurements, is performed with filtering techniques. The second step is realized with uncertainty propagation methods. The main challenge in prognosis is to take the different uncertainty sources into account in order to obtain a measure of the RUL uncertainty. There are mainly model uncertainty, measurement uncertainty and future uncertainty (loading, operating conditions, etc.). Thus, probabilistic and set-membership methods for model-based prognosis are investigated in this thesis to tackle these uncertainties.

. The ability of an extended Kalman filter and a particle filter to perform RUL prognosis in presence of model and measurement uncertainty is first studied using a nonlinear fatigue crack growth model based on the Paris' law and synthetic data. Then, the particle filter combined to a detection algorithm (cumulative sum algorithm) is applied to a more realistic case study, which is fatigue crack growth prognosis in composite materials under variable amplitude loading. This time, model uncertainty, measurement uncertainty and future loading uncertainty are taken into account, and real data are used. Then, two set-membership model-based prognosis methods based on constraint satisfaction and unknown input interval observer for linear discrete-time systems are presented. Finally, an extension of a reliability analysis method to modelbased prognosis, namely the inverse first-order reliability method (Inverse FORM), is presented. In each case study, performance evaluation metrics (accuracy, precision and timeliness) are calculated in order to make a comparison between the proposed methods.