Je souhaite remercier en premier lieu mon directeur de thèse, Pr Christian Soize. J'apprécie la chance qui a été la mienne de pouvoir travailler sous la direction d'un chercheur de son envergure. J'ai été très impressionné par la puissance et la portée de ses travaux. Je le remercie pour ses conseils précieux, ses retours sur mon travail, pour le temps qu'il a su m'accorder pour orienter mes recherches.

Je tiens ensuite à remercier Christine. Son expertise a largement contribué au succès de ce travail de thèse. Elle m'a apporté une meilleure compréhension des subtilités et des grands enjeux techniques du système ferroviaire. Elle m'a par ailleurs gratifié d'un soutien et d'une confiance dont je lui suis très reconnaissant.

Je remercie également Guillaume pour ses conseils avisés. Il a souvent apporté un regard différent et toujours très pertinent sur mon travail. Ses propositions ont eu un impact certains sur le résultat auquel j'ai finalement abouti.

Enfin, un grand merci à toute la direction Innovation & Recherche de la SNCF, en particulier le groupe PSF, à l'Association des doctorants du GPF, à tous les collègues et amis qui m'ont accompagnés pendant ces trois années de thèse.

Abstract

The work presented here deals with the development of a health-state monitoring method for high-speed train suspensions using in-service measurements of the train dynamical response by embedded acceleration sensors.

A rolling train is a dynamical system excited by the track-geometry irregularities. The suspension elements play a key role for the ride safety and comfort. The train dynamical response being dependent on the suspensions mechanical characteristics, information about the suspensions state can be inferred from acceleration measurements in the train by embedded sensors. This information about the actual suspensions state would allow for providing a more efficient train maintenance.

Mathematically, the proposed monitoring solution consists in solving a statistical inverse problem. It is based on a train-dynamics computational model, and takes into account the model uncertainty and the measurement errors. A Bayesian calibration approach is adopted to identify the probability distribution of the mechanical parameters of the suspension elements from joint measurements of the system input (the track-geometry irregularities) and output (the train dynamical response).

Classical Bayesian calibration implies the computation of the likelihood function using the stochastic model of the system output and experimental data. To cope with the fact that each run of the computational model is numerically expensive, and because of the functional nature of the system input and output, a novel Bayesian calibration method using a Gaussian-process surrogate model of the likelihood function is proposed. This thesis presents how such a random surrogate model can be used to estimate the probability distribution of the model parameters. The proposed method allows for taking into account the new type of uncertainty induced by the use of a surrogate model, which is necessary to correctly assess the calibration accuracy.

The novel Bayesian calibration method has been tested on the railway application and has achieved conclusive results. Numerical experiments were used for validation. The long-term evolution of the suspension mechanical parameters has been studied using actual measurements of the train dynamical response.

Résumé

Ce travail de thèse traite du développement d'une méthode de télédiagnostique de l'état de santé des suspensions des trains à grande vitesse à partir de mesures de la réponse dynamique du train en circulation par des accéléromètres embarqués.

Un train en circulation est un système dynamique dont l'excitation provient des irrégularités de la géométrie de la voie ferrée. Ses éléments de suspension jouent un rôle fondamental de sécurité et de confort. La réponse dynamique du train étant dépendante des caractéristiques mécaniques des éléments de suspension, il est possible d'obtenir en inverse des informations sur l'état de ces éléments à partir de mesures accélérométriques embarquées. Connaître l'état de santé réel des suspensions permettrait d'améliorer la maintenance des trains. D'un point de vue mathématique, la méthode de télédiagnostique proposée consiste à résoudre un problème statistique inverse. Elle s'appuie sur un modèle numérique de dynamique ferroviaire et prend en compte l'incertitude de modèle ainsi que les erreurs de mesures. Les paramètres mécaniques associés aux éléments de suspension sont identifiés par calibration Bayésienne à partir de mesures simultanées des entrées (les irrégularités de la géométrie de la voie) et sorties (la réponse dynamique du train) du système.

La calibration Bayésienne classique implique le calcul de la fonction de vraisemblance à partir du modèle stochastique de réponse et des données expérimentales. Le modèle numérique étant numériquement coûteux d'une part, ses entrées et sorties étant fonctionnelles d'autre part, une méthode de calibration Bayésienne originale est proposée. Elle utilise un métamodèle par processus Gaussien de la fonction de vraisemblance. Cette thèse présente comment un métamodèle aléatoire peut être utilisé pour estimer la loi de probabilité des paramètres du modèle. La méthode proposée permet la prise en compte du nouveau type d'incertitude induit par l'utilisation d'un métamodèle. Cette prise en compte est nécessaire pour une estimation correcte de la précision de la calibration.

La nouvelle méthode de calibration Bayésienne a été testée sur le cas applicatif ferroviaire, et a produit des résultats concluants. La validation a été faite par expériences numériques. Par ailleurs, l'évolution à long terme des paramètres mécaniques de suspensions a été étudiée à partir de mesures réelles de la réponse dynamique du train. High-speed trains dynamical behavior strongly relies on their suspension elements. By limiting the vertical and lateral oscillations induced by the track-geometry irregularities, they ensure the train stability and thus the ride safety. Moreover, they play a key role in passengers comfort and in structural protection by filtering the vibrations transmitted from one suspended mass to another. The suspension elements undergo damage through time. The damage can be of two types: sudden faults and gradual degradation. Sudden faults require a regular maintenance, mostly based on visual inspection. The gradual degradation leads to the replacement of certain suspension elements after a period of use. Nowadays, the replacement is triggered according to age or mileage criteria that, by definition, do not take into account the real health state of the suspension elements. These criteria are designed so that a given high percentage of these elements are still valid when replaced, as far as their mechanical properties are concerned. They are thus not optimal, since they lead to the replacement of functional suspension elements that could be kept on the train without negative impact on the dynamical behavior. A better knowledge of the real health state of the suspension would allow for replacing the suspension elements only when necessary and consequently reduce maintenance costs. This thesis is part of a project that develops a health-state monitoring solution for high-speed train suspensions, based on embedded accelerometers.

Glossary

French high-speed trains (TGV) circulating nowadays are only equipped with one type of sensors that monitor their dynamical behavior for safety purposes: the instability sensors. They trigger an alarm when they detect abnormal yaw oscillations directly threatening the ride safety. However, they are not involved in suspension monitoring. With the current development of sensors and of data-analysis technologies, more and more industrial systems and in particular transportation systems are equipped with sensors for various purposes (maintenance, system optimization, autonomous circulations...). This work is part of this trend. Embedded sensors inside a vehicle raise several issues: powering the sensors, transmitting data, maintaining the sensors... Using acceleration sensors for trainsuspension monitoring has already led to several research works and publications. The choice of acceleration sensors is guided by two main reasons. First, they are common and nonwired technologies exist, powered thanks to the train vibrations and transmitting data Frequency (Hz) Amplitude (dB) Figure 1.1: Influence of the track geometry. The three curves represent the mean value of the amplitude of the bogie vertical accelerations in the frequency domain, obtained from simulations using geometry records measured at three different years. through radio signals. Second, they allow for using a reduced number of sensors. The idea is not to equip every single suspension element, but rather to equip the main suspended masses of the train. The consequence is that the measurements do not directly provide the desired information. Several processing and analysis steps are required. The present thesis focuses on this aspect: the development of the mathematical method allowing for extracting information about the suspensions health state from the measured acceleration signals.

For train-suspension monitoring, relying solely on measurements of the train dynamical behavior is not sufficient, because of its strong sensitivity to the track-geometry irregularities. The latter consist of small displacements of the rails with respect to the theoretical track design. They constitute the main excitation source of a rolling train and, consequently, have a major influence on its dynamical behavior (see [START_REF] Perrin | Quantification of the influence of the track geometry variability on the train dynamics[END_REF], Lestoille et al., 2016a[START_REF] Lestoille | Stochastic prediction of high-speed train dynamics to long-term evolution of track irregularities[END_REF]). Track geometry is also subject to damage caused by railway traffic (see [START_REF] Bing | Development of railroad track degradation models[END_REF]) that wears and packs the track ballast. Track geometry requires regular maintenance, especially on high-speed lines, to ensure safe circulations at maximum speed. Consequently, the track-geometry irregularities evolve through time, subject to gradual degradation, regular maintenance, and replacements. It has strong implications on the suspension-monitoring problem. Indeed, a bad train dynamical behavior may be caused either by faulty suspension elements or by a degraded track geometry. For this work, a large quantity of data is used for the suspension monitoring, which implies that measurements performed over track geometries of various ages after maintenance are available. One could expect that with a sufficiently large amount of data, the averaging effect would remove the influence of the track-geometry degradation. We observed that it is not the case. We performed train-dynamics simulations along the whole line considered for this work, with the exact same conditions (identical train model, same constant speed), but using track-irregularities measurements performed at different years. As shown on Figure 1.1, the simulation results demonstrated significant variations in the average train dynamical responses from one year to another. This highlights the need for taking into account the actual track irregularities in order to perform a correct monitoring of the suspensions state.

Track-geometry irregularities measurements are included in the analysis thanks to train-dynamics simulation. The latter is performed by the commercial multibodydynamics code VAMPIRE R . Using simulation implies building a computational model, which has the advantage of explicitly accounting for the various suspension elements through the input parameters that describe their mechanical characteristics. The objective of the monitoring method is then to identify these parameters from measurements of the system response, the train dynamical behavior. We are thus confronted to a statistical inverse identification problem. The use of simulation, in conjunction with a large amount of measured data, may lead to unaffordable costs in computational resources. Such issues are also at the center of the present work and motivate the main novelties developed in the thesis.

The system is affected by various sources of uncertainty that must be taken into account to perform a robust suspension-parameter identification. They can be of different nature: the uncertainty on the model-parameters values, the modeling errors, the measurement uncertainties. Including uncertainties implies dealing with probabilistic quantities, which require adapted mathematical methods. In a nutshell, the industrial objective of this thesis is the development of a health-state monitoring method for high-speed suspensions for maintenance purposes, exploiting measurements of the train dynamical behavior and of the track-geometry irregularities, relying on a computational model and under uncertainties.

Health-state monitoring for train suspensions: short review

During the last decades, various methods have been proposed for the health-state monitoring of train suspensions. They rely on the measurements provided by embedded accelerometers. The literature focuses on the mathematical processing of the dynamicalresponse signals, which allows for extracting information about the train suspension elements.

The preliminary tests presented in [START_REF] Sunder | Operational Experiences with Onboard Diagnosis System for High Speed Trains[END_REF] demonstrate the potential of accelerometers for this purpose. The suspension elements affecting the dynamical behavior of the rolling train, it is assumed that, reciprocally, observing this behavior can provide information about the state of these suspensions. In [START_REF] Sunder | Operational Experiences with Onboard Diagnosis System for High Speed Trains[END_REF], the authors underline the interest of the analysis of the acceleration signals in both time and frequency domains.

Some important differences can be made between the methods listed in the following sections. The first difference concerns the monitoring objective: the detection of a sudden fault in one of the suspension elements or the study of their gradual degradation or wear, if they are subject to the latter. Fault detection usually considers bigger changes in the mechanical characteristics of suspension elements. The fault-detection methods work on short time scales to allow for a quick action after the fault happens. On the contrary, the monitoring of wear is performed on a longer time-scale but may look for smaller evolution of the suspensions health-state. The literature (see the following references) focuses mostly on sudden-fault detection, even though it can be argued that a fault-detection method may be used for the study of suspensions wear.

The second difference to be made is between fault detection, fault isolation, and faultmagnitude quantification. Fault detection is a first step that consists in observing a change in the train dynamical behavior, which is considered to be the consequence of a fault happening. Fault isolation is the subsequent step that consists in identifying the faulty suspension element. Fault-magnitude quantification consists in estimating the value of the mechanical characteristics of the faulty suspension to determine the damage level. Most solutions dealing with sudden suspension fault focus on fault detection and isolation, often with no clear separation between these two steps.

The third difference concerns the need of a train-dynamics model. It can provide a precious insight about the relationship between the suspensions state and the train dynamical behavior. It requires the construction of an accurate vehicle model, in which the suspension elements are explicitly represented and associated with mechanical parameters. Most health-state monitoring methods developed lately rely on simulation at some point. Therefore, the quality and representativity of train-dynamics models is of high importance for monitoring purposes. Three different categories of methods can be considered:

• the methods that rely only on the measured data and do not require any model. Because they are appropriate for fault detection only, such methods are not considered in the present thesis. Refer for instance to [START_REF] Wei | A comparative study on fault detection methods of rail vehicle suspension systems based on acceleration measurements[END_REF];

• the fault-identification methods that do not need a train model for the monitoring procedure itself, but require a prior training. The latter is usually achieved thanks to simulation results, but ideally measured data should be used. These methods are called data-driven methods and developed in Section 1.2.1.

• the fault-identification methods where the train-dynamics model is at the core of the monitoring procedure, often through the use of Kalman filters. These methods are called model-based methods and developed in Section 1.2.2.

Data-driven methods in the context of train-suspension monitoring

In [START_REF] Mei | Condition monitoring of rail vehicle suspensions based on changes in system dynamic interactions[END_REF], the authors propose a fault-detection method relying on the analysis of the cross-correlation function between different body motions. For vehicles with a symmetrical design, they show that the coupling between the different motions is small. A faulty element alters the symmetry, which results in a coupling between motions that can be observed in the cross-correlation function. Focusing on the vertical primary suspension elements, they highlighted the impact of a faulty damper on the cross-correlations between the bounce, pitch, and roll acceleration signals. In the studied case, fault isolation seems possible but requires that each type of fault be associated with its effect on the cross-correlation signals.

In [START_REF] Martinod | Identification of the technical state of suspension elements in railway systems[END_REF], the authors perform an experimental modal identification of the train dynamical behavior from the acceleration measurements. The frequency and damping ratio of several modes are identified. Thanks to simulation results using the multibody-dynamics code VAMPIRE R , regression models are identified between the frequency, the damping ratio and the value of the parameter describing the suspension state, for several modes and various degradation configurations. The regression can then be used in inverse to determine the suspension state from frequencies and damping ratio estimated from experimental modal identification.

The method proposed in [START_REF] Gasparetto | Data-driven condition-based monitoring of high-speed railway bogies[END_REF] is also based on a frequency analysis of the acceleration signal. The Random Decrement Technique (RDT) provides an approximation of the signal auto-correlation. The decomposition in Prony series then allows for estimating the natural frequencies and damping ratios of the vehicle free oscillations. The mean and standard deviation of the first natural frequency and damping ratio on various acceleration samples are then used to determine the fault type and a level of magnitude thanks to a k-NN (k Nearest Neighbors) categorization method. The k-NN initial training is achieved thanks to simulation results generated using a multibody-dynamics software.

In [START_REF] Sakellariou | On board fault detection and identification in railway vehicle suspensions via a functional model based method[END_REF], a baseline autoregressive with exogenous excitation (ARX) model is identified from experimental data. It represents the relationship between the system input (the track excitation) and the response (the accelerations) as an infinite impulse response (IIR) filter with a white noise residual. On this baseline-model structure, a functional ARX (FARX) model is identified. It represents a faulty train, with ARX coefficients depending on the fault magnitude. Only one type of suspension fault can be represented by a given FARX model. The FARX identification is performed thanks to simulation results. The fault-detection and identification procedure works as follows: the measured accelerations are injected in the FARX model; the magnitude of the fault is determined by minimizing the FARX residual variance. When this fault magnitude is above a threshold, the fault type is validated if the residual is sufficiently uncorrelated.

Model-based methods in the context of train-suspension monitoring

Most model-based methods rely on Kalman filters to perform a health-state monitoring of train suspensions. The purpose of linear Kalman filtering is to follow the time-evolution of a dynamical system using a linear model describing the system behavior while taking into account the information provided by regular measurements of some observable quantities.

A basic linear Kalman filter can be described thanks to the equations:

X k+1 = [A k+1 ]X k + U k+1 , (1.1) Y k = [H k ]X k + V k , (1.2)
where X k and Y k are the state-space vector and the observation vector at step k, [A k ] the state-transition matrix, [H k ] the observation matrix, U k and V k are independent centered Gaussian vectors representing the process noise and the observation noise. At each time step k, the state-space vector is updated using Eq. (1.1) and then conditioned by the available experimental observation using Eq. (1.2). In the case of train dynamics, state-space vector X k usually gathers the degrees of freedom of the various masses constituting the vehicle and their derivatives, while observation vector Y k consists of the accelerations measured in different locations. The mechanical parameters of suspension elements are used to compute matrices

[A k ] and [H k ].
A Kalman filter is designed to estimate the time-evolution of the state-space variables. On the contrary, the parameters are supposed to be fixed and given. Adaptations are thus required in order to perform Kalman-based parameter identification.

A first approach described in [START_REF] Alfi | Condition monitoring of suspension components in railway bogies[END_REF] is to consider the parameters of interest as state-space variables. An evolution equation must be defined for the parameters. Moreover, using an augmented state-space vector yields nonlinear equations that are solved by using an Extended Kalman Filter (EKF). The latter uses the Jacobian matrix of the system to perform a linearization around the current value of the state-space vector. However, the simultaneous estimation of several parameters using EKF appears to achieve bad results. Consequently, in [START_REF] Alfi | Condition monitoring of suspension components in railway bogies[END_REF], they propose to use multiple EKF in parallel, each of them performing the identification of a single suspension parameter. A modification of the train response can lead to changes in the estimated value of several parameters. The most probable faulty suspension is then identified from a likelihood estimation of the measurements with respect to the various faulty models.

Several other papers propose to use multiple Kalman filters in parallel. In [START_REF] Jesussek | Fault detection and isolation for a full-scale railway vehicle suspension with multiple Kalman filters[END_REF], multiple Kalman filters are run in parallel, each one modeling a fault type, with fixed magnitude. It means that the state-transition matrix and the observation matrix are modified for each filter to take into account the degraded suspensions parameters. The fault identification is once again performed from a likelihood estimation (averaged over several time-steps) of the measurements with respect to the different models.

The IMM (Interacting Multiple Models) algorithm presented in [START_REF] Hayashi | Fault Detection of Railway Vehicle Suspension Systems Using Multiple-Model Approach[END_REF]] is similar but additionally includes mode mixing. The parallel Kalman filters are no longer isolated but interact with one another. For a given filter and at each time step, the input state-space vector is a combination of the output state-space vector of all filters at the previous time step. This combination is based on the modes likelihood and on given transition probabilities between modes. In a subsequent paper [START_REF] Mori | Condition Monitoring of Railway Vehicle Suspension Using Multiple Model Approach[END_REF], the authors propose a model-updating procedure to adapt the baseline model as a fault is detected and to allow for identifying simultaneous faults. The idea is to work with several groups of Kalman filters in parallel, each one focusing on one type of suspension. After a modification in the suspensions parameters detected by one group, the algorithm allows for updating the models of the other groups.

The method proposed in [Li et al., 2004, Li et al., 2007], although relying on multiple Kalman filters, differs from the previous ones because the associated models are not chosen in advance to represent fixed fault types and magnitudes. In a similar fashion to EKF-based parameter estimation, the aim is to estimate the probability distribution of an augmented state-space vector gathering the initial state-space variables (the degrees of freedom of the multibody-dynamics model) and the suspension parameters. The authors then propose to separate the estimation of the probability distributions of the state-space variables and of the parameters thanks to marginalization. A Rao-Blackwellized particle filter is used to represent the probability distribution of the parameters by a weighted sample, while linear Kalman filters are used to estimate the state-space variables.

The fault-identification method proposed in [START_REF] Jesussek | Fault detection and isolation for a nonlinear railway vehicle suspension with a hybrid extended kalman filter[END_REF] relies on a single EKF, for the case when the various degrees of freedom are decoupled. More precisely, the degrees of freedom associated with one suspension element of interest must not be linked to other suspension elements. The use of EKF allows for considering dampers with nonlinear properties. The fault-identification procedure focuses on the Kalman filter innovation. Its norm is minimized by adding an error to the degrees of freedom associated with a suspension element of interest. The magnitude of this error determines whether the corresponding suspension element must be considered as faulty.

In [Li and Goodall, 2004], the authors also focus on the Kalman filter innovation. A fault is detected when the weighted-sum squared residual (averaged over several time steps), equivalent to the log-likelihood of the observations, crosses a given threshold. A possible fault-identification procedure is then suggested, which consists in analyzing the norm of the power spectral density of the different components of the Kalman innovation.

In [START_REF] Liu | An efficient recursive least square-based condition monitoring approach for a rail vehicle suspension system[END_REF], instead of a Kalman filter, a closely related time-domain filter known as Recursive Least Squares (RLS) is used. RLS is an algorithm able to identify the parameters of an input-output linear system by filtering the error signal between the measured and modeled outputs. It requires measurements of both the input and the output signals. In this paper, the acceleration of the wheelsets is used as input, so the track irregularities do not need to be measured nor to be included in the model.

Most of the cited articles rely on Kalman filters to propose a fault-detection and isolation method for abrupt changes in the suspensions characteristics. This is underlined by the magnitudes of the fault that are considered (often 50% or 100%). Considering that this type of failures represents a risk for the ride safety, these methods mostly focus on a quick detection and exploit acceleration signals of a few seconds.

In the present thesis, we are interested in the long-time degradation of the suspension elements. We exploit large measurements data-bases, corresponding to rides along several hundreds of kilometers of track, to ensure the robustness of the monitoring results.

In the cited articles, the train models usually consist of a single car, and consider only certain degrees of freedom and the associated suspension elements, for which, in general, the nonlinearities are not taken into account. On the contrary, the model-based approach used in this thesis includes the full representation of an entire train, including nonlinearities for wheel-rail contact and in suspension mechanical characteristics. It allows for the simultaneous monitoring of several suspensions elements, with or without strong couplings. This requires simulations using a commercial nonlinear multibodydynamics code, VAMPIRE R . The modeling is thus more physical and allows for finer identifications.

In the cited articles, the inclusion of the various uncertainties of the system is often lacking. Even though Kalman filters naturally take measurement noise into account, no model error is considered in general. The necessity to include the various uncertainties of the system motivates our choice to rely on a Bayesian approach. Many papers model the track-geometry irregularities as a simple colored random noise. As explained previously, we are convinced of the importance of using the real geometry whenever possible. For this work, actual measurements of the track-geometry irregularities are used instead of considering them as a basic random input.

Finally, certain of the proposed method are tested on numerical experiments but rarely on actual acceleration measurements. In this thesis, we presents the results of the application of our method on several sets of measurements.

Scientific problematic

In general, a physical system can be represented by a mathematical model providing an output response depending on an input excitation and input parameters. The goal of this thesis is to identify the system input parameters given the system output response, using a computational model. Such problem belongs to the class of inverse problems, and is more precisely equivalent to a calibration problem. In a deterministic framework, problems of this type are generally solved using optimization approaches, the solution being the parameters values minimizing a given distance between the computational-model output response and the available measurements of this response. A central point of this thesis is however to propose a calibration approach in a statistical framework, which includes various uncertainties in the system. For this statistical inverse problem, the Bayesian framework provides a rigorous and robust mathematical approach to combine stochastic models with experimental data. Fundamental aspects of Bayesian calibration of computational models are detailed in [Kennedy and O'Hagan, 2001]. The authors propose an exhaustive description of the various uncertainties one may encounter when exploiting a computational model of a real physical system. Their calibration method relies on Gaussian-process (GP) modeling: the unknown functions are represented by Gaussian processes. In particular, they efficiently replace the expensive computer code that is only known thanks to a limited set of observations.

The approach developed in this thesis has many similarities: the objective is to perform the calibration of a nonlinear expensive computational model, in a Bayesian framework. However, the procedure developed in [Kennedy and O'Hagan, 2001] is hardly applicable because of the functional nature of the considered input excitation (the trackgeometry irregularities) and output response (the acceleration signal representing the train dynamical behavior). About the functional excitation input, it has been shown in [START_REF] Perrin | Track irregularities stochastic modeling[END_REF] that the track irregularities can be modeled as a nonstationary stochastic process. Because of their statistically rich content, even with dimensionality reduction techniques such as truncated Karhunen-Loeve expansion, an acceptable model requires hundreds of coefficients. Representing the computational model by a GP model would then require to index the latter on a very high dimensional definition set. This is not a feasible solution, because of the sparsity of the observations available to identify the GP model. About the functional output response, the authors in [Kennedy and O'Hagan, 2001] solve the issue by augmenting the GP argument. In the applicative case of the present thesis, the output response consists of acceleration signals in the frequency domain. Their solution would then consist in making a frequencydependent GP. The main difference with the spatial example treated in the paper lies in the fact that each of the available acceleration measurements covers the entire frequency band of interest. Considering the frequency as a GP argument would then multiply the number of available observations by the discretization size of the measurements, which is not without consequences. Indeed, the covariance matrix of the observation needs to be computed and inverted. In the present case, the solution proposed in [Kennedy and O'Hagan, 2001] would lead to a very large and probably ill-conditioned matrix. Their method cannot be applied in our case without significant adaptations.

Consequently, we propose a novel method of Bayesian calibration of an expensive computational model with functional input and output. Compared to the method in [Kennedy and O'Hagan, 2001], we do not rely on a GP surrogate model of the computational model. Instead, a GP surrogate model of the likelihood function is built. The latter is at the core of Bayesian approaches, and has the main advantage of being scalar and not depending on the system excitation input, but only on the model parameters. The main consequence is that the uncertainty induced by the approximation inherent in the use of a surrogate model is not included in the likelihood function anymore and must be handled differently. Moreover, dealing with a random likelihood function involves a random probability density function that is not conventional. These issues are addressed in this thesis.

Overview of the adopted strategy

The industrial goal of this thesis is the development of a health-state monitoring method for high-speed train suspensions. The key choices made to address this problem are:

• exploiting measurements of the train dynamical behavior performed by embedded accelerometers. The objective is to avoid reducing the information provided by the sensors, considering only a few scalar indicators. It implies working with functional quantities of interest.

• performing a simultaneous identification of multiple parameters associated with different suspension elements, in order to avoid possible compensating effects.

• including the measured track-geometry irregularities in the analysis.

• taking uncertainties into account to ensure the identification robustness. In the same perspective, a large measurement database is exploited.

In order to consider the whole variety of suspension elements and the associated degrees of freedom, a multibody-dynamics computational model is necessary. Simulation is required to study the influence of the track-geometry irregularities on the train dynamical behavior. We rely on the commercial multibody-dynamics code VAMPIRE R . The entire train is represented by a multibody model in which the various suspension elements are accounted for and associated with one or several mechanical parameters. Identifying the parameters of a computational model from measurements of the system output is equivalent to a calibration problem. The Bayesian framework is well-suited to take into account the system uncertainties. Hence, we deal with a Bayesian calibration problem.

Using a computational model in addition to a large measurements database raises an issue of computational resources. The use of a surrogate model is a solution to address such an issue. We propose to combine a Bayesian approach with Gaussian-process surrogate models to solve the parameters identification problem. Such combination has already been achieved in the past when only scalar (or low dimensional) outputs are considered. We developed a novel approach allowing for functional outputs to be considered.

In order to perform the Bayesian calibration, a GP surrogate model of the likelihood function is built. This surrogate model is a stochastic process indexed by the computational-model parameters. The classical likelihood function is then replaced by a random algebraic representation. The uncertainties induced by the introduction of a surrogate model will be called surrogate-model uncertainty. They can no longer be included in the likelihood-function computation. They must be taken into account thanks to the introduction of an additional step. For a given realization of the random surrogate model, called a trajectory, Markov Chain Monte Carlo (MCMC) is used to construct realizations of the posterior probability density function (PDF) of the model parameters conditioned by this trajectory. The set of all realizations for all trajectories allows for estimating the posterior PDF. In addition to this novel Bayesian calibration method, an approximation method for the trajectories of the GP surrogate model is proposed.

Outline of the thesis

The second chapter focuses on the mathematical aspects of the method proposed. Section 2.1 introduces the Bayesian formalism and the likelihood function that lies at the core of Bayesian approaches. Gaussian-process surrogate modeling is then presented in Section 2.2. We underline the fact that this section also details the new method of goal-oriented approximation of a trajectory of a GP surrogate model in high dimension. Section 2.3 presents how GP surrogate models have been used for calibration according to the literature. The novel method we propose that allows for considering functional outputs is detailed in Section 2.4. Finally, Section 2.5 recalls the principle of Markov Chain Monte Carlo and details the TMCMC algorithm used for this work.

The third chapter deals with the industrial problem, starting by the definition of the problem in Section 3.1. The system characteristics and quantities of interest are then detailed in Sections 3.3 and 3.4. The identification of the output-predictive error accounting for the system uncertainties is detailed in Section 3.5. The two last sections present the results obtained for the application of the new identification method. Numerical experiments are used for validation in Section 3.6. Section 3.7 analyses the results obtained from actual measurements of the train dynamical behavior.

Chapter 2 Bayesian calibration using a random surrogate model of the likelihood function

This chapter focuses on the mathematical aspects of Bayesian calibration. They are presented independently from the industrial application, although the specificities of the latter drive the major aspects of the novel method we developed: a Bayesian calibration method of an expensive computational model with functional input and output. Section 2.1 introduces the Bayesian formalism and especially the likelihood function. Gaussian-process surrogate models and Markov Chain Monte Carlo are classical mathematical tools described in Sections 2.2 and 2.5. We propose a short review of two types of calibration method relying on GP surrogate models in Section 2.3. Their limits lead to the novel Bayesian calibration method with GP surrogate model of the likelihood function we describe in Section 2.4. A secondary new development introduced in this thesis consists of an approximation technique for the trajectories of a GP surrogate model, described in 2.2.4.

Classical Bayesian formalism adapted to the parameter calibration of a computational model

Bayesian calibration is a type of statistical inverse problem, consisting in estimating the probability distribution of the model input parameters from measurements of the output response. Section 2.1.1 describes the various quantities of interest and the associated uncertainties. Section 2.1.2 presents the Bayesian formalism that is used to solve the statistical inverse parameter identification. Figure 2.1 shows a diagram summarizing the calibration procedure.

System definition and uncertainties

The first input of the model is the random excitation X. For the application, X will represent the spatial discretization of the track-geometry irregularities. It is defined by a unknown a probability distribution for which a set of realizations {x mes,i } i is known. The second input of the model is the random system parameters W that are of interest. For the application, W will represent the mechanical parameters of the train suspensions. A prior probability density function (PDF) p prior W of random parameter W is assumed to be constructed using the maximum entropy principle under the constraints defined by the available information (see for instance [Soize, 2017]). The support of p prior W is denoted as C W . The objective of Bayesian calibration is to update the prior PDF in order to obtain a posterior PDF p post W that takes into account the information provided by the experimental data.

The computational model is represented by a deterministic mapping h depending on random quantities X and W. The model output response Y sim is written as:

Y sim = h(X; W) . ( 2

.1)

Mapping h being deterministic implies that if the excitation and the parameters are fixed, then the model output response is deterministic. The quantity Y sim is random but

{Y sim | X = x, W = w} = h(x; w) (2.2)
is deterministic.

The physical system is represented by a model affected by parametric uncertainty (mechanical properties, geometrical parameters, etc) and model uncertainty induced by modeling errors (simplification, linearization, etc), see for instance [Soize, 2017]. We suppose that they can globally be taken into account thanks to the introduction of a random error ε sim that represents the distance between the model response and the "real" physical response Y real :

Y sim = Y real + ε sim . (2.
3)

The measurements may also be subject to uncertainties. For instance, the exact setup of the sensors may not be accurately known. Formally, we assume that they can be represented in the following way:

Y mes = Y real + ε mes , (2.4)
where Y mes is the measured system response, and ε mes is a random variable representing the measurement uncertainties. In practice, for Bayesian calibration, the random measured response Y mes is only known through a set of realizations {y mes,i } i .

In practice, it is often difficult to separate the measurement uncertainties from the model uncertainties. They may however be taken into account globally thanks to a single random output-predictive error B representing the distance between the measured response and the model output response. It is supposed to be independent from excitation X and parameters W. This error has to be identified from the available experimental data (more details are given for the industrial application case in Section 3.5). The stochastic model of the system response denoted as Y is defined by

Y = Y sim -ε sim + ε mes = h(X; W) + B .
(2.5)

By construction, the stochastic model of the output response Y is a predictive model of the measured response Y mes .

Introducing the likelihood function

Since Y depends on parameters W, the principle of Bayesian calibration is to adapt the probability distribution of W so that the probability distribution of Y "fits" the available realizations of Y mes .

In this work, we consider the particular situation where experimental data do not only consist of a set of independent measurements {y mes,i } i of the output response, but of a set of ν ≥ 1 joint measurements of the input excitation and of the corresponding output response {(x mes,i , y mes,i )} 1≤i≤ν , which are simultaneously measured.

Mathematically, the Bayesian calibration procedure consists in estimating the probability distribution of W, represented by PDF p post W , conditioned by the input-output couple (X, Y), for which set {(x mes,i , y mes,i )} 1≤i≤ν is given. We first consider the case where a single measurement (x mes , y mes ) is available (ν = 1, the corresponding superscript is omitted). For w in the support C W of density p prior W , we have

p post W (w) = p W | X,Y (w | x mes , y mes ) (2.6) ∝ p Y | W,X (y mes | w, x mes ) p W | X (w | x mes ) (2.7) ∝ L(w) p prior W (w) , (2.8)
where

L : C W → R + w → p Y | W,X (y mes | w, x mes ) (2.9)
is the likelihood function, and where p prior

W (w) = p W | X (w | x mes ) because W is assumed independent from X.
The measurements are supposed to be fixed and omitted from the function L arguments. Equation (2.7) results from the application of the Bayes formula.

When ν > 1, independent copies

(X i , Y i ) 1≤i≤ν of couple (X, Y) must be intro- duced. Random vector W is now conditioned by X i = x mes,i , Y i = y mes,i , 1 ≤ i ≤ ν.
The posterior PDF is then equal to (PDF arguments are omitted for simplicity):

p post W = p W | X 1 ,...X ν ,Y 1 ,...Y ν (2.10) ∝ p Y 1 ,...Y ν | W,X 1 ,...X ν p W | X 1 ,...X ν (2.11) ∝ ν i=1 p Y i | W,X i p prior W .
(2.12) Equation ( 2.12) is deduced from the assumption that for any

(i, j) ∈ {1, . . . ν} 2 , if i = j, (X i , Y i ) is independent of (X j , Y j ).
The likelihood function is then defined by

L : w → ν i=1 p Y i | W,X i (y mes,i | w, x mes,i ) . (2.13)
The likelihood function may take values very close to zero, because of dimensionality aspects: the normalization constant of a probability density function associated with a high-dimensional random quantity goes to zero when dimension goes to infinity. The log-likelihood function is thus commonly used instead:

L :      C W → R w → ν i=1 log p Y i | W,X i (y mes,i |w, x mes,i ) . (2.14)
From Eq. ( 2.5), it can be deduced that:

{Y | X = x, W = w} = h(x; w) + B . (2.15)
The log-likelihood function can then be rewritten as

L : w → ν i=1
log p B y mes,ih(x mes,i , w) .

(2.16) Equation (2.8) is at the core of Bayesian calibration. The computation of the likelihood function is not self-evident and requires the use of adapted numerical methods. It depends on the complexity used for describing the stochastic model of B, and on the computational cost for evaluating h(x; w) at given x and w. In our specific case (a deterministic computational model and an input excitation introduced via measurements corresponding to the output response measurements), Eq. ( 2.16) shows that it requires running the computational model and estimating the PDF of the output-predictive error B as many times as the number of available measurements. We choose to identify B as a Gaussian vector, which simplifies the computation of the corresponding PDF.

Assuming that we are able to compute the log-likelihood function L(w) for any w in C W , the following step is to apply Markov Chain Monte Carlo (MCMC, see Section 2.5) with Eq. (2.8) to generate samples distributed as p post W . The moments of posterior PDF p post W are estimated using these samples. PDF p post W can be estimated using nonparametric statistics, for example the kernel-density estimation method.

Gaussian process surrogate modeling

This section presents the mains aspects of Gaussian-process (GP) surrogate modeling. The main features are given in Section 2.2.1. For the complete mathematical details, we refer the reader to Appendix A. The question of noisy observations and of the interest of adding a nugget effect is analyzed in Section 2.2.2. Section 2.2.3 presents two refining methods from the literature aiming at locally reducing the approximation error of the GP surrogate model. As explained in the subsequent Section 2.4, we need to draw trajectories (realizations) of the GP surrogate model. The classical method of simulation on a grid do not suit our problem because of the process dimension. Thus, we propose a new alternative method in Section 2.2.4, which consists in approximating a given trajectory by further conditioning the GP surrogate model.

General principles

Gaussian process surrogate modeling of a deterministic real target function y, defined on a subset X of R q , consists in constructing a Gaussian stochastic process Y indexed by X , whose mean function is a good approximation of y and for which statistical fluctuations are small. The latter is then conditioned by a set of n observations (x obs i , y(x obs i )) 1≤i≤n on the graph of y, for a given training set {x obs i } 1≤i≤n that is a subset of X , in order to obtain the surrogate model Y . The Bayesian approach [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF] for the construction of Y is detailed in Appendix A. A second approach based on the mean-square error minimization exists [START_REF] Sacks | Design and Experiments of Computer Experiments[END_REF]. Both lead to similar results. The second approach is worth mentioning because it shows that for any x in X , the mathematical expectation E{ Y (x)} constitutes the best approximation of y(x) in a mean-square sense, based on the available knowledge about target function y provided by the observations. The function x → E{ Y (x)} is called the Kriging predictor. The randomness of the surrogate model accounts for the uncertainty about the actual value of the target function. The variance Var{ Y (x)} quantifies the accuracy of the approximation of y(x) by E{ Y (x)}.

GP surrogate modeling is popular because of the closed-form expressions it provides. It also has the advantage to quantify the lack of knowledge about the target function through its variance, which allows for implementing adaptive learning algorithms (see for example the refining algorithm in Section 2.2.3). As shown on Figure 2.2 (left graph), a GP surrogate model is exact at the observation points and is interpolating in X outside these points.

The evaluation of the target function value at the points of the training set is usually a numerically expensive step when building a surrogate model. Consequently, the training has to be designed so that it maximizes the information provided about the target function. A common choice is to design a space-filling training set that will allow for a good representation of the global behavior of the target function in its whole definition set. Appendix A.6 gives an example of such space-filling training set: the optimized Latin Hypercube Sample. When particular features of the target function are studied, the initial surrogate model can be exploited to perform a subsequent refining step. Refining algorithm examples are presented in Section 2.2.3.

Noisy observations and nugget effect

With very little change in the identification procedure, we can consider that the observations are noisy in the construction of the GP surrogate model, as illustrated in Figure 2.2, right graph. Appendix A.5 details the mathematical implications of noisy observations. The noise vector is considered centered and Gaussian, with covariance matrix proportional to identity matrix. It is often referred to as "adding a nugget effect", because it supposes adding a constant value to the diagonal terms of the observations covariance matrix.

The first reason for adding a nugget effect is because it improves numerical stability [START_REF] Ababou | On the condition number of covariance matrices in kriging, estimation, and simulation of random fields[END_REF]. In the construction of the GP surrogate model, the observations covariance matrix must be inverted. Even though the latter is theoretically a positive-definite symmetric matrix, it may be ill-conditioned, which affects the quality of the results or even makes the matrix numerically impossible to invert. Adding a nugget effect automatically improves the matrix conditioning.

In [START_REF] Gramacy | Cases for the nugget in modeling computer experiments[END_REF], the authors argue and demonstrate on simple examples that adding a nugget effect can be beneficial for other reasons, especially when building a surrogate model of a computational model. Necessary assumptions such as stationarity of the centered Gaussian process may be incorrect and affect the quality of the surrogate model, especially when the observation points are sparse and the target function contains small scale variations. The authors point out that the interpolation provided by the zero nugget assumption does not always guaranty a good representation of the target function. They argue that the smoothing effect provided by the addition of a nugget effect offers overall a better representation, in the sense of coverage (the target function is contained in the confidence band of the surrogate model), and of fitting to the large scale behavior of the target function. The introduction of a nugget was tested in our multidimensional case, with good results, as illustrated on Figure 2.3. The GP model without nugget effect displays oscillations, but is smooth with a nugget effect, without a significant increase of the width of the confidence region. Both graphs represent a slice of the GP surrogate model along the first dimension of set X . The dotted curve corresponds to the mean function and the surrounding colored area to the 95% quantiles. The observations are supposed exact on the left graph and noisy on the right one. They are not represented because none belongs to this particular slice of multidimensional set X .

Local refining of the GP surrogate model around the maximum of the target function y

Building the GP surrogate model with a space-filling training set allows for learning about the global behavior of the target function on the whole set X . However, it usually results in a lack of local precision of the surrogate model, especially in high-dimension, when the training set tends to be sparse. If particular features of the target function need to be explored, the initial surrogate model can be used to propose new observation points.

After the evaluation of the target function at these new points, the surrogate model is updated using the new information about the target function. The new observation points are determined by maximizing a criterion φ : X → R corresponding to the feature of interest. Criterion φ depends on the GP surrogate model. The iterative refining algorithm can be defined as follows. The stopping condition is a threshold ρ on the maximum of criterion φ. If the latter is lower than threshold ρ, then the surrogate model cannot be significantly refined any further around the feature of interest. Usually, the number of iterations is also limited, depending on the affordable number of calls to target function y.

while the maximum number of calls to y is not reached do Determine

x new ∈ X that maximizes φ if φ(x new ) > ρ then Evaluate y(x new )
Update Y with the new observation (x new , y(x new )) else Break while loop end if end while

In this section, two criteria are presented, originally used with the exact same algorithm to solve an optimization problem using GP surrogate models [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF].

Refining the surrogate model may also be seen as building an adaptive design of experiment [START_REF] Picheny | Adaptive Designs of Experiments for Accurate Approximation of a Target Region[END_REF].

In [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF]], the objective is to maximize an expensive target function using a GP surrogate model. The authors introduce the EGO algorithm (Efficient Global Optimization). Because each call to the target function is expensive, each new observation should bring as much information as possible about the target function while refining the maximum location. Criterion φ takes into account the variance of the surrogate model in order to offer a balance between getting closer to the maximum and learning about the target function where information is lacking. Consequently, the proposed criterion is the "expected improvement" (EI):

φ EI : x → E max(0, Y (x) -y max ) , (2.17)
where y max is the current observed maximum of the target function and Y is the GP surrogate model, conditioned by the available observations, including the possible new observations performed at the previous iterations of the algorithm. For this approach, the observations of target function y are supposed to be exact. The point in which y is maximum is looked for among the points of the training set (including the new observation points provided by the algorithm), at which the value of y is exactly known.

When observations are noisy, EI is no longer adapted. A similar criterion adapted to noisy observations called Knowledge Gradient for Continuous Parameters (KGCP) is proposed in [START_REF] Scott | The Correlated Knowledge Gradient for Simulation Optimization of Continuous Parameters using Gaussian Process Regression[END_REF], derived from the Knowledge Gradient criterion introduced in [START_REF] Frazier | The knowledgegradient policy for correlated normal beliefs[END_REF]. It can be seen as a generalization of the expected improvement. The centered Gaussian noise affecting the observations is denoted as ε. In the general case, its variance depends on x. The KGCP criterion is written as

φ KG : x → E max ξ∈X E{ Y (ξ) | Y (x) + ε(x)} -max ξ∈X E Y (ξ) .
(2.18)

In Eq. ( 2.18), the expectation

E{ Y (ξ) | Y (x) + ε(x)} is a function of random variable Y (x) + ε(x).
Consequently, it is random itself, which justifies the calculation of the expectation of its maximum. The idea is to look for the point that is expected to maximize the maximum of the Kriging predictor when added to the observations. As such, the criterion is hardly computable because of the two maximizations on the whole set X . That is why an approximation of the knowledge gradient is used in [START_REF] Scott | The Correlated Knowledge Gradient for Simulation Optimization of Continuous Parameters using Gaussian Process Regression[END_REF]. The maximum on X is approximated by the maximum on the training set augmented with the point in which φ is evaluated. In Eq. ( 2.18), the expression of φ KG in is kept, in which set X replaced by the set C * (x) = {x} ∪ {x obs i } 1≤i≤n . For each iteration, a new observation is added to the training set, so the size of C * (x) increases.

Goal-oriented approximation of GP surrogate model trajectories using a conditioning set

In this section, we propose a goal-oriented approximation method of the GP surrogate model trajectories. Here goal-oriented means that the trajectories are accurately approximated in a given subset X * of interest of X but not in all the set. Subset X * is chosen so as to cover certain features of GP surrogate model Y (.). An example designed to approximate the GP trajectories in the vicinity the maximum of the GP mean function is detailed at the end of Section 2.4. Let X c be a finite and nonordered subset of X * made up of N c points {x c j } 1≤j≤Nc . Subset X c is chosen in order to approximate a given trajectory y(.) by the mean function of stochastic process Y (.) conditioned by the value of y(.) at the points of X c . The latter is then called the conditioning set. Its construction depends on the definition of X * . For a good representation of X * , we suggest drawing conditioning set X c from a uniform distribution on X * . This can be achieved using TMCMC (see Section 2.5.2) with the indicator function on X * , as long as this indicator function can easily be computed.

We denote as m(.) and C(., .) the mean function and the covariance function of stochastic process Y (.). It is assumed that the covariance operator, whose kernel is covariance function C(., .), is positive definite. We denote as Y c and y c the vector gathering the values of process Y (.) and trajectory y(.) at the points of X c :

Y c =    Y (x c 1 ) . . . Y (x c Nc )    , y c =    y(x c 1 ) . . . y(x c Nc )    . (2.19)
The approximate y(.; y c ) of trajectory y(.) is defined as

y(x; y c ) = E Y (x) | Y c = y c (2.20) = m(x) + r c (x) t [C c ] -1 (y c -m c ) , x ∈ X , (2.21)
where

m c = E{Y c } =    m(x c 1 ) . . . m(x c Nc )    , (2.22) r c (x) = E {(Y (x) -m(x))(Y c -m c )} =    C(x, x c 1 ) . . . C(x, x c Nc )    , (2.23) [C c ] = E (Y c -m c )(Y c -m c ) t , [C c ] ij = C(x c i , x c j ) . (2.24)
As a matter of fact, trajectory y(.) is only known on X c . It is even unknown before the value of y c is randomly drawn, as a realization of Y c , using the classical generator for Gaussian vector. The approximate y(.; y c ) is completely dependent on y c , which is why we explicitly indicated this dependency. Note that if two different trajectories of Y (.) had the same value on X c , then their approximates with the proposed method would be the same. Conditioning set X c must be sufficiently space-filling in X * so that the variance of the random variable { Y (x) | Y c = y c } is small for any x in X * .

In Appendix A.4, the covariance function will be parameterized introducing hyperparameters. They can be identified using the Bayesian approach, in which they are assumed to be random, and defined by a prior probability distribution. They must be conditioned by the available data. The probability distributions of the system quantities of interest are then expressed conditioned by the hyperparameters and should be integrated with respect to them. In practice, this Bayesian approach cannot be applied analytically for the hyperparameters of the covariance function.

Instead of the Bayesian approach, these hyperparameters will be identified by optimization according to a criterion assessing how well the Gaussian process is fitting the data provided by the training set. The criterion used is the Maximum Likelihood Estimation (MLE). This approach guarantees that the mean function of the GP surrogate model is close to the target function and that the statistical fluctuation around this mean are small enough (small variance).

Using GP surrogate models for calibration

This section presents the existing solutions for the calibration of an expensive computational model thanks to GP surrogate models. Section 2.3.1 comments the Bayesian calibration method developed in [Kennedy and O'Hagan, 2001] and the reasons why it is unadapted to the problem we are dealing with in this thesis. Section 2.3.2 mentions two papers dealing with calibration problems that inspired the novel method proposed in this thesis.

Bayesian calibration with a surrogate model of a scalar output

Let us first consider the case of a simulation-based model with scalar output, represented for example by the equation

Y = h(X; W) + ε (2.25)
where Y is the model of the random output response, mapping h represents the computational model, X is the random vector of the input excitation, W is the random vector of the model parameters, and ε is an additive noise representing both the model uncertainties and the measurement errors. Random variable ε is modeled by a Gaussian centered random variable of variance σ 2 ε . Compared to [Kennedy and O'Hagan, 2001], we consider here that the model uncertainties and the measurement errors can be taken into account globally, and that they do not depend on the input. For a single measurement y mes corresponding to an input x mes , the likelihood function is written as

L(w) = p Y | X,W (y mes | x mes , w) (2.26) = p N y mes ; h(x mes ; w), σ 2 ε , (2.27)
where w ∈ C W and p N (.; µ, σ 2 ) stands for the normal density of mean µ and variance σ 2 :

p N (.; µ, σ 2 ) : u → 1 σ √ 2π exp - (u -µ) 2 2σ 2 .
(2.28)

Using MCMC requires numerous evaluations of the likelihood function L and consequently numerous calls the computational model h. An expensive computational model makes the procedure unaffordable. An efficient solution to adress this numerical issue is to rely on a GP surrogate model that consists in representing the computational model by a Gaussian process. Building the GP surrogate model consists in conditoning a prior Gaussian process by observations on a training set, as explained in Section 2.2. The subsequent Bayesian calibration step also consists in conditoning a probabilistic model, by the available measurements this time.

We suppose that the computational model is represented by the GP surrogate model H : (x; w) → ĥ(x; w) + Z(x; w) , (2.29) where ĥ is the deterministic Kriging predictor and Z a centered GP of variance σ 2 Z independent from noise ε. The previous model of the output response then becomes Y = ĥ(X; W) + Z(X; W) + ε , (2.30) and the associated likelihood function is

L : w → p N (y mes ; ĥ(x mes ; w), σ 2 Z (x mes ; w) + σ 2 ε ) , (2.31)
If several measurements are considered, the expression cannot be simply expressed as a product similar to Eq. ( 2.13), because the independence assumption does not hold anymore. More specifically, if one considers two different measurements (x mes,i , y mes,i ) and (x mes,j , y mes,j ), then for any w in C W , random variables Z(x mes,i ; w) and Z(x mes,j ; w) are not independent. The expression of the likelihood function then takes the form of a multivariate Gaussian density function.

The major interest of building a GP surrogate model of the computational model, apart from the computational efficiency, is that the surrogate model uncertainty can be readily introduced in the likelihood function. This is highlighted by Eq. ( 2.31) in which the variance of the GP surrogate model σ 2

Z is added to the variance of the noise σ 2 ε . The calibration procedure is not modified. The Bayesian formalism is kept, which allows for a correct evaluation of the posterior uncertainty on the parameters W.

The first limit of the approach introduced by [Kennedy and O'Hagan, 2001] concerns the dimension of X, which controls the dimension of the set the GP is indexed by. The higher the dimension, the sparser the observations, and the larger the number of hyperparameters that need to be identified. In this thesis, the computational model works with functional excitation inputs corresponding to the track-geometry irregularities. Their discretized form has a dimension of several thousands. The work of [START_REF] Perrin | Track irregularities stochastic modeling[END_REF] has shown that they are characterized by a rich statistical content. Consequently, dimensionreduction techniques such as truncated Karhunen-Loeve expansion still require hundreds of coefficients to provide a good approximation. The functional nature of the input of computational model tremendously increases the dimension of the GP surrogate model, which thus makes the representation given in Eq. ( 2.29) impossible to build. An alternative is the construction of a surrogate model of the computational model that takes only parameter w as argument, but that is specific to one realization of the input excitation. This solution implies that one surrogate model H i (.) is build for each input measurement x mes i . For 1 ≤ i ≤ ν, H i (.) then represents the mapping h(x mes i ; .). The increased number of surrogate models however counters the computational efficiency they are supposed to provide.

The second limit concerns the functional nature of the system output response Y . Indeed, computational models often output quantities depending on time, space or frequency. In the problem treated in this thesis, we study the train dynamical response as a function of the frequency ω. In [Kennedy and O'Hagan, 2001], the authors circumvent this difficulty by augmenting the input variable x with the variables on which the output response depends. It means considering the model

Y = h(ω; X; W) + ε .
(2.32)

Formally, this approach does not significantly complexify the overall method. It requires building a surrogate model H indexed by a set of slightly higher dimension. The dimension increase is equal to the dimension of ω, one in our case. In the paper, the authors present a geotechnique example with a dependance on a 2-D space variable. However, this augmentation approach may raise a significant issue regarding the number of observations used for the construction of the surrogate model. Indeed, the computational model generally provides a discretized response in m sampling points of ω. If ω becomes an argument of the GP surrogate model, then the number of observations used for its construction is multiplied by m. For fixed values of x and w, the value of Y provided by the computational model at each sampling point of ω is considered as a single observation.

Dealing with hundreds of thousands of observations is a major difficulty since the covariance matrix of observations has to be manipulated. Moreover, numerous observations that would be too close to each other may result in an ill-conditionning of the covariance matrix. An adaptation for functional outputs is proposed in [Perrin, 2018], which consists in identifying the statistical dependence structure of Y for different values of ω independently from the surrogate model covariance depending on x and w. They are then combined thanks to the Kronecker product.

Deterministic calibration with a functional output

For the calibration of a computational model with functional outputs, for which one evaluation is expensive, another approach is proposed in [START_REF] Ranjan | Inverse Problem for a Time-Series Valued Computer Simulator via Scalarization[END_REF] and [START_REF] Pratola | Fast sequential computer model calibration of large nonstationary spatialtemporal processes[END_REF]. In both papers, the calibration is solved as an optimization problem. Optimal parameters are obtained by minimizing a cost function ∆ : C W → R depending on the parameters and representing the distance between the measurements and the modeled output. To cope with the fact that one evaluation with the computational model is expensive, the authors rely in both cases on the optimization algorithm EGO [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF]] briefly presented in 2.2.3. This algorithm relies on a representation of function ∆ by a GP surrogate model. The use of GP surrogate models is made possible by the scalarization of the problem through the construction of function ∆.

In [START_REF] Ranjan | Inverse Problem for a Time-Series Valued Computer Simulator via Scalarization[END_REF], cost function ∆ is defined as the logarithm of the L 2 -norm of the difference between the measurement and the modeled output. In [START_REF] Pratola | Fast sequential computer model calibration of large nonstationary spatialtemporal processes[END_REF], the approach is based on the framework provided by [Kennedy and O'Hagan, 2001]. They also work on the difference between the measurement and the modeled output. Two stochastic models for this difference are defined: a simple model, equivalent to a centered Gaussian white noise, and colored Gaussian model with a non-zero mean function and a given time correlation. The basic assumption is that when the parameters are set to their "real" values, the distance between the measurements and the modeled output should fit equivalently on both models. When the parameters differ from their "real" values, the distance should fit better on the colored model. Function ∆ is then defined from the likelihood ratio of the white and colored models.

Although both methods achieve interesting results, they cannot be considered as Bayesian approaches. They identify optimal parameters values that provide the best fit between the experimental data and the model, according to the chosen distance. However, they do not take into account the uncertainties in the computational model and their propagation in the calibration procedure. In [START_REF] Pratola | Fast sequential computer model calibration of large nonstationary spatialtemporal processes[END_REF], the authors propose to estimate the distribution of the position of the minimum of the GP surrogate model. They achieve that by drawing and minimizing trajectories of the GP surrogate model. From this distribution, they deduce confidence intervals around the calibrated parameters values. This analysis allows for evaluating the uncertainty on the calibrated parameters stemming from the approximation by a surrogate model, but still ignores the other sources of uncertainty, such as model uncertainties or measurement errors.

Novel method of Bayesian calibration with GP surrogate model of the likelihood function

The objective of this thesis is the development of a Bayesian calibration method for which one evaluation of the computational model is expensive, in presence of a functional output. As shown in [START_REF] Ranjan | Inverse Problem for a Time-Series Valued Computer Simulator via Scalarization[END_REF] and [START_REF] Pratola | Fast sequential computer model calibration of large nonstationary spatialtemporal processes[END_REF], scalarization combined with GP surrogate modeling is an efficient solution. In the Bayesian formalism (presented in Section 2.1), a natural scalarization is provided by the likelihood function. Hence, we perform the calibration relying on a GP surrogate model of the likelihood function. The objective is to remain in the Bayesian framework while taking advantage of GP surrogate modeling. This approach raises various questions that will be addressed in this section:

• how should a GP model of the likelihood function be built?

• how to exploit this random surrogate model for calibration?

• how to take into account the uncertainty induced by the use of a surrogate model for the estimation of the posterior probability distribution of the parameters?

The idea of building a surrogate model of the likelihood function is not new. For instance, [START_REF] Orlande | Approximation of the likelihood function in the bayesian technique for the solution of inverse problems[END_REF] proposes an interpolation of the likelihood function using radial basis functions, while [START_REF] Dietzel | Bayesian inference of a lake water quality model by emulating its posterior density[END_REF]] also relies on GP surrogate modeling. However, the main consequence of these approaches, the introduction of a new type of uncertainty, the surrogate-model uncertainty, is never taken into account.

In addition, it should be noted that we focus on the likelihood function itself rather than on the product L×p prior W of the likelihood function by the prior PDF of the parameters. The reason is simply because we consider a prior PDF p prior W uniform on C W : p prior W ∝ 1 on C W . It may then be removed from the fundamental equation Eq. (2.8). If the prior PDF was different, the method detailed in the following section could be used by simply replacing the likelihood function L by the product L × p prior W , as long as the value of the prior PDF can be computed anywhere in admissible set C W .

GP surrogate modeling of the likelihood function

As explained before, we choose to work with the log-likelihood function L instead of the likelihood function L. The first reason is the dimensionality issue that makes the likelihood function take values very close to zero. It also circumvents the fact that the likelihood function has to respect positivity. In addition, the variations of the log-likelihood function tend to be slower (especially when the likelihood function corresponds to a very peaked density, as it is the case here), which make this function easier to represent by a GP surrogate model.

We choose the Matérn-5 2 function for the covariance shape, because it provides smooth results. A smooth surrogate model makes the calibration easier, especially if MCMC is used. For the same reason, we introduce a nugget effect, even though the observations of the log-likelihood function are a priori not noisy. Justifications for this choice are detailed in Section 2.2.2.

The GP surrogate model of log-likelihood L(.) is denoted as L(.; Θ). Random variable Θ need not to be explicitly defined but is introduced to indicate the randomness of the surrogate model. A realization θ of Θ yields the realization L(.; θ) of stochastic process L(.; Θ), which is a trajectory of the GP surrogate model. At any point w in C W , the Kriging predictor is denoted as the mathematical expectation E Θ {L(w; Θ)} and constitutes the best approximation of L(w) in a mean-square sense. The variance Var Θ {L(w; Θ)} quantifies the approximation error of the Kriging predictor. It should be noted that L(.; Θ) corresponds to the stochastic process Y (.) introduced in Section 2.2.1.

Calibration strategies

The GP surrogate model of the log-likelihood function may be used in various different ways to perform the calibration. First, it can simply be maximized, using EGO [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF]] if the observations are considered exact, KGCP [START_REF] Scott | The Correlated Knowledge Gradient for Simulation Optimization of Continuous Parameters using Gaussian Process Regression[END_REF] if a nugget effect is added. This approach would provide the most probable parameters values. The likelihood is equivalent to a distance between the measurements and the modeled output. This approach is thus equivalent to [START_REF] Ranjan | Inverse Problem for a Time-Series Valued Computer Simulator via Scalarization[END_REF], Pratola et al., 2013] and suffers the same lack: it does not follow the Bayesian formalism and only provides optimal parameters instead of a complete posterior PDF. To remain in the Bayesian framework, the posterior PDF must be estimated via MCMC.

MCMC must be applied on a deterministic likelihood function. The GP surrogate model being random, it cannot be used directly as such. To perform the posterior PDF estimation with MCMC, the most straightforward solution is to use the Kriging predictor E Θ {L(.; Θ)} instead of L. In the rest of this manuscript, we refer to this approach as the KP (for "Kriging predictor") method. Once the GP surrogate model is built, this solution is easy to implement. It can provide useful results, especially if only the mean value of the model parameters are needed. However, this solution does completely ignore the surrogate model uncertainty, which represents the approximation error intrinsic in such modeling. The uncertainty on the calibrated parameters would only stem from the system uncertainties represented by random vector B.

In the following section, we propose a method that takes into account the uncertainty of the GP surrogate model in order to estimate more correctly the calibration error. In the rest of this manuscript, we refer to this new approach as the MCT (for "Monte Carlo on the trajectories") method.

Novel MCT method: Monte Carlo on the trajectories of the GP surrogate model

If the deterministic log-likelihood function L(.) is replaced by the random surrogate model L(.; Θ), the corresponding posterior PDF becomes random as well. The random posterior PDF corresponding to L(.; Θ) can be expressed using the conditioning by Θ as

w → p W | X,Y,Θ (w | x mes , y mes , Θ) . (2.33)
For simplicity, we only present the case of a single measurement (x mes , y mes ). The explanations are completely equivalent in the case of multiple measurements.

The quantity we are looking for remains unchanged: the conditional PDF of parameters W knowing (X, Y). The rule of conditional probabilities states that it is equal to the expectation with respect to Θ of the previous PDF:

p post W (w) = p W | X,Y (w | x mes , y mes ) (2.34) = E Θ p W | X,Y,Θ (w | x mes , y mes , Θ) (2.35) ≈ 1 N N j=1 p W | X,Y,Θ (w | x mes , y mes , θ j ) , (2.36)
where {θ j } 1≤j≤N correspond to N independent realizations of Θ, and so {L(.; θ j )} 1≤j≤N correspond to N independent trajectories of the GP surrogate model. The estimation of the N PDF p W | X,Y,Θ (w | x mes , y mes , θ j ) in the right-hand side of Eq. ( 2.36) requires the use of the MCMC method. In practice, for each j, MCMC provides samples in C W , distributed following p W | X,Y,Θ (w | x mes , y mes , θ j ). If, for all j, these samples have the same number of points, their simple concatenation is distributed following p post W . To apply MCMC on a given trajectory L(.; θ j ) of the log-likelihood surrogate model, it is necessary to be able to compute this trajectory value anywhere in C W . Indeed, MCMC randomly picks new candidate points in C W and computes the likelihood function value in these points to build the resulting samples (see Section 2.5). The trajectory computation is performed using the method developed in Section 2.2.4. We denote as L(.; θ j ) the approximate of trajectory L(.; θ j ) using this method. Another approach could consist in iteratively conditioning the Gaussian process by the previously drawn samples at the previous steps of the MCMC. However, it implies the manipulation of a full covariance matrix whose size would increase at each step. Moreover, this covariance matrix is bound to become ill-conditioned because of the concentration of points in the vicinity of the maximum of the trajectory. This approach has been tested and has been judged too time consuming.

Section 2.2.4 develops a goal-oriented approximating method for trajectories of a GP surrogate model focused on a specific subset of the admissible set C W , denoted as X * in Section 2.2.4. The approximation method is designed to be accurate in X * but not elsewhere. Subset X * should then be chosen with care. In the present case, we define X * as the region of C W where "the relative value of the likelihood function is high". In the region of low relative likelihood, no samples will be drawn, an accurate representation of the trajectories is then useless. We define X * as the set where the surrogate model has a probability higher than a tolerance ρ ∈]0, 1[ to be greater than where its mean function is maximum:

X * = {w ∈ C W | Proba {L(w; Θ) > L(w max ; Θ)} ≥ ρ} (2.37)
with w max = arg max w∈C W E Θ {L(w; Θ)}. For any w fixed in C W , let Λ(w) be the Gaussian random vector

Λ(w) = L(w; Θ) L(w max ; Θ) . (2.38)
It is a Gaussian vector because the process L(.; Θ) is Gaussian. From the construction of the surrogate model, one can directly compute the mean vector m Λ (w) and the covariance matrix [C Λ (w)] of vector Λ(w) for any w in C W . The definition of set X * can be rewritten

X * = w ∈ C W | F N a t m Λ (w) a t [C Λ (w)]a ≥ ρ , (2.39)
where a = [1 , -1] t and F N : R → [0, 1] is the cumulative distribution function of the normal density

F N : x → 1 √ 2π x -∞ exp - u 2 2 du .
(2.40)

Posterior probability density function estimation with MCMC

This section presents the Markov Chain Monte Carlo (MCMC) class of algorithms, see for instance [START_REF] Givens | Computational statistics[END_REF]. MCMC is used at the end of the calibration procedure in order to estimate the posterior probability density function (PDF) of W.

We present the TMCMC algorithm that has been used in this thesis and the classical Metropolis-Hastings algorithm on which TMCMC is based.

The purpose of MCMC is to generate independent realizations of a random variable W, defined by its given PDF π. The principle is to simulate a Markov process that admits the probability distribution π(w)dw as its invariant probability distribution. PDF π is assumed to be known up to a normalizing constant. MCMC is generally used when random variable W cannot be defined as a function of a random variable for which a generator is known. From the drawn realizations of W, statistics of π can be estimated.

The classical Metropolis-Hastings algorithm

The algorithm presented in this section was proposed in [Hastings, 1970]. We sum up the main features of the algorithm in the continuous case. The different possible values of random variable W are considered as possible states of the Markov process. The transition probability density between state v and state w is written

p(v, w) = q(v, w)α(v, w) (2.41)
where q is the proposal distribution and α is defined as:

α(v, w) = 1 if π(w) π(v) ≥ 1 , π(w) π(v)
otherwise.

(2.42) At each step i, a proposal state v * is first randomly drawn using density w → q(v i-1 , w) centered on the previous state v i-1 of the Markov process. Function q depends only on the distance vw so that q(v, w) = q(w, v). A classical choice for q is the Gaussian density, which yields a proposal that can be written as

v * = v i-1 + σN , (2.43)
where σ ∈ R + * and N is a centered Gaussian random vector. Number σ controls the step size, the average distance between two subsequent states of the chain. Depending on the density ratio r = π(v * ) π(v i-1 ) the proposal v * is then kept or rejected. If r ≥ 1, the proposal is kept as the new state v i of the Markov process. Otherwise, it has a probability 1 -r to be rejected, in which case, state v i is set equal to the previous state v i-1 .

The Metropolis-Hastings algorithm can be summed up as follows:

Choose starting point v 0 for i = 1 : n do Draw a proposal v * as a realization of

v i-1 + σN Compute r = π(v * ) π(v i-1 ) if r ≥ 1 then Set v i = v * else Draw u from a uniform distribution on [0, 1] if u ≤ r then Set v i = v * else Set v i = v i-1 end if end if end for Return {v i } 0≤i≤n
Set {v i } 0≤i≤n is asymptotically distributed according to target PDF π. To allow for better exploring the support of π, several chains are often run in parallel using different starting points. We tested this approach, but ended up with a strong dependency between the final samples and the chosen distribution of the starting points of the parallel chains.

A specific algorithm: TMCMC

The Transitional MCMC (TMCMC) algorithm [START_REF] Ching | Transitional Markov Chain Monte Carlo Method for Bayesian Model Updating, Model Class Selection, and Model Averaging[END_REF]] is a MCMC algorithm derived from the Metropolis-Hastings algorithm. The principle is to start from samples distributed according to the prior PDF, and making this set of points gradually evolve toward samples distributed according to the posterior PDF. This algorithm has been developed to sample from multimodal or very peaked distribution, which is often challenging using the classical Metropolis-Hastings algorithm. In this Section, we adopted the same notations as in Section 2.1.

The algorithm works in m steps by drawing m successive sets of samples {w j,k } 1≤k≤Ns , with index j denoting the step number. Integer N s is the target number of samples, set by the user. For 1 ≤ j ≤ m, samples {w j,k } k are distributed according to PDF p j defined for w ∈ C W by p j (w) ∝ L(w) q j p prior W (w) (2.44) where 0 = q 0 < q 1 < • • • < q m-1 < q m = 1. One can immediately notice that p 0 = p prior W and p m = p post W . Exponent q j is determined as

q j = arg min q∈]q j-1 ,1] |CV j-1 (q) -1| (2.45)
where CV j-1 (q) is the coefficient of variation (equal to the standard deviation divided by the mean value) of set {L(w j-1,k ) q-q j-1 } 1≤k≤Ns . With this choice of q j , the idea is to control the dispersion of the updated likelihood values at the point of the current set of samples. The exponents q j are determined iteratively, which means that the number m of steps in the algorithm is unknown at the beginning. At the beginning of each step j, 1 ≤ j ≤ m, each point w j-1,k of the previous set of samples, distributed according to p j-1 , is affected a normalized weighting coefficient

ᾱj,k = α j,k Ns =1 α j, (2.46)
where α j,k = L(w j-1,k ) q j -q j-1 , 1 ≤ k ≤ N s . The next N s samples are drawn from {w j-1,k } k according to the probability distribution defined by these normalized weighting coefficients. These new samples are then distributed according to p j . In order to avoid the repetition of identical samples, MCMC steps are applied to disturb the samples while keeping the same distribution. The Metropolis-Hastings algorithm is used, with a Gaussian density for the proposal distribution. Its covariance matrix is estimated from the samples {w j-1,k } k as 2.47) where m j = Ns k=1 ᾱj,k w j-1,k . The factor β is introduced to control the step size. In [START_REF] Betz | Transitional Markov Chain Monte Carlo: Observations and Improvements[END_REF], some improvements are proposed in the algorithm. Among those is the adaptation of coefficient β, depending on the rejection rate of the MCMC phase, in order to control this rejection rate.

[Σ j ] = β 2 Ns k=1 ᾱj,k (w j-1,k -m j )(w j-1,k -m j ) t (
The algorithm can be summed up as follows:

Log-likelihood: compute L(w Chapter 3

Parameter identification for high-speed train suspensions from in-line acceleration measurements

This chapter presents the industrial problem that motivated the development of the novel Bayesian calibration method, and the results provided by the application of this method. Section 3.3 describes the studied train, focusing on the suspension elements of interest, how they are modeled and how the model is parameterized. Section 3.4 presents the computational model and describes the input and output quantities of interest. Section 3.5 details the identification procedure for the output-predictive error, whose purpose is to represent the various uncertainties of the system, essential to perform a robust monitoring.

The last two sections deal with the application of the monitoring method and present the results. The methods, choices, and results are discussed throughout these sections. Consequently, no particular section will be devoted to the discussion. Section 3.6 focuses on the validation of the monitoring method thanks to numerical experiments, while Section 3.7 presents its exploitation on real measurements.

Setting the problem

A rolling train constitutes a nonlinear dynamical system. It oscillates under the excitation of the track-geometry irregularities. Acceleration signals measured at different points in the train are the observable quantities of the system, and constitute the train dynamical response. The system depends on a large variety of parameters. Among them, we focus on the parameters that describe the mechanical characteristics of certain suspension elements. Our objective is to determine the functional health state of these suspension elements. We thus consider these mechanical parameters as variables describing the suspension health state.

The train dynamical response is influenced by the suspension parameters. Reciprocally, the latter can be identified in inverse from measurements of the train dynamical response. To solve this parameter identification problem, we rely on a computational model able to simulate the behavior of a rolling train on a given track geometry. The computational model allows for including the track-irregularities measurements in the analysis, and has the great advantage to explicitly model all suspension elements. To sum Computational model

Real system

Simulated train dynamical response

Measured train dynamical response

Suspension mechanical parameters

Track-geometry irregularities Parameter identification Figure 3.1: Diagram of the train-dynamics system up, the objective is to identify the suspension mechanical parameters of the computational model that simulates the train dynamics from joint measurements of the track-geometry irregularities and of the train dynamical response (see Figure 3.1). To this end, we adopt the Bayesian calibration approach presented in Chapter 2 and we have implemented the MCT method described in Section 2.4.

The available data to perform the parameter identification are:

• numerous measurements of the train dynamical response, on several track stretches, at different dates.

• the corresponding measurements of the track-geometry irregularities.

• a complete multibody-dynamics model of the train.

• the engineering design of the train suspension elements and their mechanical characteristics.

The track curvature has a strong impact on both the train dynamical behavior and the track-geometry irregularities. Consequently, the track has been divided into several stretches depending on the curvature in order to gather more statistically coherent datasets. The different types of track stretches are: straight lines, full curves, curve entrances, and curve exits.

Definition of the train-dynamics system

The railway system consists of a train in interaction with a ballasted track through the rails. A railway ballasted track is a complex mechanical structure: the rails, the sleepers, the ballast, the various layers of substructures, and the ground on which it is seated. It may be described and modeled in different ways and on different scales depending on the studied phenomenon.

The propagation of vibrations induced by railway traffic in the track and the surrounding environment is a widely studied mechanical problem [START_REF] Degrande | A numerical model for groundborne vibrations from underground railway traffic based on a periodic finite elementboundary element formulation[END_REF], Clouteau et al., 2001a, Kuo et al., 2016, Galvín et al., 2018]. The latter requires a fine description of the track to account for the propagation medium from the source (the wheelrail contact forces) to the structures of interest affected by the vibrations (surrounding buildings or the track itself, for instance). All the components of the track are subject to variability. The use of stochastic models and adapted methods allow for performing uncertainty quantification. The uncertainties may be considered at every level: the rail unevenness [START_REF] Lombaert | Quantification of uncertainty in the prediction of railway induced ground vibration due to the use of statistical track unevenness data[END_REF], the ballast [de Abreu Corrêa et al., 2017], the ground [START_REF] Clouteau | Stochastic simulations in dynamic soil-structure interaction[END_REF], Ta et al., 2010], etc.

In this thesis, we are focusing on train dynamics. A full modeling of the various components of the track is not used. In the train-dynamics model, the excitation at the wheels level is driven by:

• the track design.

• the track stiffness that accounts for the deflection of the different track components under the train static and dynamical loadings.

• the track-geometry irregularities, measured simultaneously with the train dynamical response, under loading (for more details, see Section 3.4.2).

The train-dynamics model consists of:

• the rail and wheel profiles.

• the wheel-rail coefficient of friction driving the wheel-rail contact force.

• the multibody-dynamics modeling, including the suspension elements (for more details, see Section 3.4).

Description of the train and of the suspensions of interest

This study focuses on the TGV Réseau, a french one-floor high-speed train. This type of train consists of eight passenger cars and two motor cars, one at each end. Two car junctions are studied: the first between the two first passenger cars, the second between the two last passenger cars. Depending on the orientation of the train, a given equipped junction may be located at the head or the rear of the train. To perform the calibration, for a given junction, we always consider measurements of the train dynamical response in both situations. Indeed, the dynamical behavior of the train significantly varies between the head and the rear. The most obvious illustration is the yaw oscillations that have a greater amplitude at the rear of the train. Certain suspension elements (typically the yaw dampers) may then undergo different solicitations depending on the orientation of the train. Considering both locations provides complementary data about a given junction.

For each junction, one sensor is located on a carbody, and one on the bogie. They measure both vertical and lateral accelerations. Since the sensors locations do not coincide with the centers of mass, the measured accelerations correspond to combinations of the different body motions. In addition to these eight acceleration signals, the location of the train along the track as well as its speed are recorded.

The train is equipped with many types of suspension elements, and often several elements of the same type. They can be divided into three groups according to their location: • the primary suspension gathering the elements linking the wheelsets to the bogie frame.

• the secondary suspension gathering the elements linking the bogie frame to the carbody.

• the inter-carbody suspension gathering the elements linking one carboby to another.

For this study, the suspension elements of interest are (see Figure 3.2):

• the axial joints integrated in the trailing arms, part of the primary suspension.

• the airsprings, part of the secondary suspension.

• the vertical primary dampers;

• the yaw dampers, part of the secondary suspension.

• the upper inter-carbody dampers.

More details about the behavior and models for these elements are provided in the following sections. This selection has been determined based on the needs expressed by the maintenance workshops: we focused on the elements involved in the train safety, requiring frequent inspection or inducing high maintenance costs. Moreover, we ruled out the elements whose influence on the accelerations measured by the sensors is too small (based on simulation results). It should be noted that other sensors configurations may be more sensitive to the degradation of these elements.

The mechanical characteristics of these elements are defined according to various parameters. The mechanical parameters are used to described the health state, considered here as equivalent to the functional state of the suspension elements. From the initial sensitivity analysis, we have figured out that with the given sensors configuration, we have not been able to locate a fault on a specific element among the other elements of the same type in the same suspension. For each junction, in order to minimize the dimension of the parameterization, it has been assumed that the set of parameters describing one type of suspension element remains the same for all the elements of this type (in the considered junction). The admissible set for each parameter is defined from the specifications detailed in the vehicle engineering design, as an interval centered around the nominal value. Since the objective is monitoring possible degradations of suspension elements, the interval width is set to stretched specification margins.

Axial joints

On the studied bogie, the motion of the wheelset with respect to the bogie frame is constrained by trailing arms. The connection between these arms and the bogie frame is ensured by an axial joint that is related to the vertical motion of the wheelset and contains an elastomeric stiffness. This stiffness works in torsion. It comes in addition to the primary vertical coil spring (which is not considered in the parameter identification).

The axial joint is modeled by a 1-D viscoelastic without memory rheological component (similar to a Zener model) consisting of a spring in parallel with a branch made of a damper and a spring in series, as represented in Figure 3.3. This model allows for representing, in a simplified way, the softening of the material at low frequencies. Stiffness K 0 is considered as the parameter to be identified, while stiffness K 1 and damping C 1 are determined as functions of K 0 [START_REF] Deltarail | Vampire Help Manual[END_REF]. Experience shows that 2C 1 is equal to at least 5% of critical damping with a parallel stiffness of 2K 0 (this corresponds to a wheelset of mass M , linked to two axial joints in parallel), and that K 1 is approximately equal to 2K 0 .

For a translation motion, C 1 is written as

C 1 = 0.05 × 2 √ 2K 0 M 2 . (3.1)
If we now consider a rotational stiffness, with a trailing arm of length , the previous equation becomes

C 1 = 0.05 × 2 √ 2K 0 M 2 . (3.2)
with K 0 in N×m/rad and C 1 in N×m×s/rad.

Airsprings

Airsprings are nowadays a very common secondary suspension type for passenger trains. Their ability to filter vibrations plays an important role in the ride comfort. Basically, an airspring consists of an elastomeric membrane under pressure supporting the vertical load of the carbody, linked by a surge pipe to a reservoir. They are mainly involved in the vertical motion (and consequently in the roll and pitch motions) of the carbody through Various rheological models of airsprings can be found in the literature (see for instance [Kraft, 2012]). They usually require tuning numerous parameters in order to reproduce the complete behavior of this type of suspension. We have chosen to rely on a much simpler model that represents the airspring as constant stiffnesses and dampers in parallel for the six degrees of freedom between the carbody and the bogie. Only the vertical damping coefficient and stiffness are considered for the parameter identification, all the others being fixed.

Dampers

Train-suspension systems are equipped with numerous dampers. In this study, we focus on the vertical primary dampers, the yaw dampers, and the upper inter-carbody dampers. For railway systems, dampers usually have a nonlinear mechanical behavior, with a damping rate decreasing with velocity. Figure 3.4 provides an example of a damper nonlinear characteristic. The latter is often approximated by a bilinear curve. In order to limit the number of parameters to identify, the parameter for each type of damper is defined as a multiplicative factor applied to the given nominal nonlinear characteristics.

The vertical primary dampers are designed to damp the pump, pitch, and roll motion of the bogie frame with respect to the wheelsets. The yaw dampers are designed to damp the yaw motion of the bogie. The yaw instability can be a cause of derailment, the yaw dampers thus have a safety function. The inter-carbody dampers, despite being oriented longitudinally, mostly damp the pitch and yaw motions between carbodies.

Multi-body simulation for train dynamics

This section briefly introduces the computational model on which the proposed suspension monitoring solution is based: the multibody-dynamics simulation software VAM-PIRE R . One should note that this computational model is used as a black-box, meaning that we do not have access to the system mechanical equations. The following sections present the principle of multibody modeling, the computational model input (the trackgeometry irregularities) and output (the train dynamical response).

In general, in the figures given in this section as illustrations for various physical quantities involved in the analysis, no axes scales are indicated for confidentiality reasons.

Description of the train multibody model

A multibody model consists of rigid bodies linked together by mechanical joints. In the case of the train model, the rigid bodies represent the wheelsets, the bogies, and the carbodies, while the mechanical joints represent the various suspension elements. The flexible modes of the rigid bodies are not taken into account. For train-dynamics simulation, a wheel-rail contact law must also be provided. It is used to determine the force applied to a wheelset as a function of its displacement and velocity. This nonlinear contact law depends on the wheel and rail profiles. In this thesis, the Kalker model of rolling contact [Kalker, 1990] is adopted. More details for the construction of a TGV multibody model with VAMPIRE R can be found in [Kraft, 2012].

Such multibody model contains numerous parameters : body masses and inertia, mechanical properties of the joints, relative positions of the bodies and the joints, contact parameters, etc. For the suspension-monitoring application, we focus on a small number q of parameters associated with the suspension elements of interest listed in Section 3.3. Following the model introduced in Section 2.1, they are gathered in the random vector W, with values in the admissible set C W , subset of R q . The admissible set C W simply consists of the set product of the parameters admissible intervals.

Description of the input: the track-geometry irregularities

The track-geometry irregularities consist of the small-scale displacements of the rails with respect to the track design. They are of four types: lateral, vertical, gauge and cross-level irregularities. For a track stretch of length S, they can be denoted as {x(s) ∈ R 4 , s ∈ [0 , S]}, argument s representing the curvilinear abscissa following the track design. An illustration of the four irregularities is provided in figure 3.5. Combined with the track design and the train speed, they are used to compute the time-varying displacement condition imposed at the wheel-rail interface for each wheelset.

It has been shown in [START_REF] Perrin | Track irregularities stochastic modeling[END_REF], and reused for other measurements in [START_REF] Panunzio | Construction of a stochastic model of track geometry irregularities and validation through experimental measurements of dynamic loading[END_REF], that the track-geometry irregularities can be modeled as a nonstationary R 4 -valued random field {X(s), s ∈ [0 , S]} indexed by distance interval [0 , S]. Though the irregularities are represented by random field X, the calibration procedure presented here need not use the model developed in [START_REF] Perrin | Track irregularities stochastic modeling[END_REF], because the irregularities are directly measured. The available measurements on various track stretches are considered as realizations of random field X. For the track-geometry irregularities, the measurement method used by SNCF shows a very good reproducibility. Consequently, the measurement noise affecting these irregularities is considered negligible compared to the other sources of uncertainty.

Description of the output: the train dynamical response

The train dynamical response consists of n = 4 acceleration signals in two specific points in the train, along the vertical and lateral axes (in the axis system attached to the train). These points correspond to the locations of the embedded sensors, on the carbody and the bogie at a given carbody junction.

The simulation is performed step by step in the time domain. The output time signals are transformed into the frequency domain. More precisely, in order to avoid systematic phase-shift between the measured and the simulated signals, the studied quantity is the amplitude of the accelerations in the frequency domain. This amplitude is taken in dB to characterize the resonances as well as the antiresonances of the system.

For 1 ≤ k ≤ n, {a k (t) , t ∈ [0 , T ]} denotes the k th -acceleration signal in the time domain for a duration T . Let Ω be the frequency band of interest. the corresponding response in the frequency domain, { a k (ω) , ω ∈ Ω}, is defined by

a k (ω) = 10 log 10 T 0 1 √ T a k (t)e -iωt dt . (3.3)
Figures 3.6 and 3.7 present an illustration of the four acceleration signals in the time domain. Since the spectral content of the accelerations in the carbody and in the bogie are significantly different, the signals are not represented on the same scales. Figure 3.8 presents the statistical characteristics of the carbody lateral acceleration in the frequency domain. The mean function and quantiles have been estimated from measurements performed on multiple track stretches, each one considered as an independent realization of the stochastic train dynamical response.

Using the relationship between the curvilinear abscissa s and time t (depending on the train speed), the train-dynamics computational model is represented by the deterministic (3.4)

It associates the response a(ω) = h(ω; x, w) in the frequency domain with the irregularity signal x and the vector w of suspension mechanical parameters. As detailed in Section 2.1.1, the various uncertainties of the system must be taken into account to perform a robust parameter identification. Because we rely on a blackbox computational model, the model uncertainties and measurement errors are globally taken into account thanks to the random output-predictive error B. The latter must be identified using the available information provided by the simultaneous measurements of {x(s) , s ∈ [0 , S]} and { a(ω) , ω ∈ Ω}. The adopted identification procedure is detailed in Section 3.5. Output-predictive error B is added to the simulation output in order to obtain the stochastic model of train dynamical response, denoted as {Y(ω), ω ∈ Ω}:

Y = h(.; X, W) + B . (3.5)
This model is a R n -valued stochastic process indexed by frequency band Ω. Each component corresponds to one of the four studied acceleration signals.

The definition given by Eq. (3.5) includes the random representation of the track irregularities and the suspension parameters. However, in practice, the latter must always be fixed to deterministic values x and w to run the computational model h. Consequently, once B has been identified, we only have access to

{Y | X = x, W = w} = h(.; x, w) + B .
(3.6)

Identification of the output-predictive error

The identification of the output-predictive error (OPE) is performed from a reference set of measurements for which the model parameters are known. This approach relies on the strong hypothesis that the OPE is independent from the train parameters, from the track-geometry irregularities, and does not evolve on a long-time scale. Nevertheless, the identification of the OPE is performed independently for the different locations of the studied bogie (at the head or the rear of the train) and for each class of track stretches, depending on the curvature types of the track design. The parameters values must be known in order to identify an OPE that allows for parameters identification. The OPE must not include the variability of the train response stemming from the uncertainty on the parameters of interest. If it did, the addition of the OPE would erase the influence of the parameters of interest and thus make the parameters identification impossible. In the present case, the reference set corresponds to measurements performed immediately after a major maintenance operation during which the suspension elements have been renewed. For these measurements, their mechanical characteristics are assumed to be nominal.

In this section, we detail the construction of stochastic process B. This stochastic process is defined as a Gaussian process that is then completely defined by its mean function and its covariance function. The estimates of these two functions will be performed by using a limited number of realizations. There is a difficulty related to the noninvertibility of the covariance matrix of the discretized process B because the methodology presented in Chapter 2 (see Eq. (2.16)) requires the explicitation of the Gaussian PDF of B.

In this section, we consider that all stochastic processes are discretized, yielding random vectors. The µ frequency sample points of frequency band Ω are denoted as (ω i ) 1≤i≤µ . For any vector-valued stochastic process, the frequency sampling of its components is concatenated to obtain a single column vector. For instance, the frequency sampling of stochastic process {Y(ω), ω ∈ Ω} is the vector of dimension nµ

[Y 1 (ω 1 ) . . . Y 1 (ω µ ) . . . Y n (ω 1 ) . . . Y n (ω µ )] t .
(3.7)

Available measurements defining the reference set

The reference set, denoted as U 0 , contains independent joint measurements of the trackgeometry irregularities and of the train dynamical response, denoted respectively x ref,i and y ref,i , for 1 ≤ i ≤ ν 0 , performed on ν 0 track stretches:

U 0 = {(x ref,i , y ref,i )} 1≤i≤ν 0 . (3.8)
For 1 ≤ i ≤ ν 0 , let then b ref,i be the difference between the measured train response and the simulated one (with nominal parameters) on the i th track stretch:

b ref,i = y ref,i -h(x ref,i , w 0 ) , (3.9)
with w 0 the nominal values of the train parameters W. The set {b ref,i } 1≤i≤ν 0 is considered as a set of independent realizations of output-predictive error B for which the stochastic model is constructed in the next sections.

Stochastic model of B and its identification

The stochastic model of the output-predictive error B in constructed by using the Maximum Entropy principle of Information Theory under the constraints defined by the available information as follows. Stochastic process {B(ω) , ω ∈ Ω} is with values in R n . There is no additional information concerning the support of the probability distribution of B, which is a priori all the set R n (since the amplitude of the train dynamical response is studied on a log-scale, it is not necessary that response Y respects positivity). As explained above, we now only consider the discretized version of B as a real-valued random Gaussian vector of dimension nµ.

The rest of the available information consists of the empirical estimates of the mean vector m and of the covariance matrix [ C] using the set of realizations {b ref,i } 1≤i≤ν 0 , which are computed as (3.11) where matrix

m = 1 ν 0 ν 0 i=1 b ref,i , (3.10) [ C] = 1 ν 0 -1 [ B][ B] t ,
[ B] = [b ref,1 -m . . . b ref,ν 0 -m]
gathers the centered realizations. We consider the case for which the number ν 0 of available realizations is less than the dimension nµ of the realizations. As a consequence, matrix [ C] is not invertible. Let r denote the rank of [ C] such that r ≤ ν 0 < nµ.

Validation of the method on numerical experiments

Ideally, the validation of the KP and MCT methods proposed in Section 2.4 should be performed from experimental results. The calibration method should be applied to real measurements of the train response in order to obtain the parameters posterior PDF. In parallel the mechanical characteristics of the suspension elements of the equipped train should be experimentally measured for comparison. Because of its cost and logistical complexity, such validation campaign could not be carried on. We relied on numerical experiments to validate the proposed methods. The principle of the validation is explained in Section 3.6.1. In Section 3.6.4, we analyze the results of the application of the KP and MCT methods and compare the two methods. The EGO algorithm mentioned in Section 2.3.2 has also been tested, as briefly described in Section 3.6.3. In addition, we give the main features on the the GP surrogate model of the loglikelihood function on which all these methods are based in Section 3.6.2.

Numerical experiment principle

In the present case, a numerical experiment consists of simulated train responses that are used as experimental data. They are generated using actual measurements of the trackgeometry irregularities on several track stretches and known degraded suspension parameters. Moreover, an independent realization of the output-predictive error B is added to the response signal on each track stretch in order to generate a quantity as close as possible to an actual measurement. The numerical experiments allows for validating the calibration procedure: the procedure is applied on the virtual train response, the calibration results can then be compared to the reference parameters used to generate the response. The "reference" parameters injected in the computational model in the context of a numerical experiment have nothing to do with the "reference" date and set mentioned in Section 3.5, used for the identification of the output-predictive error.

We suppose that a set of ν 1 track irregularities measurements {x mes,i } 1≤i≤ν 1 is available. The validation procedure from a numerical experiment can then be summed up as follows:

1. Choose vector w 1 of reference parameters;

2. Run the simulation on the ν 1 track stretches with w 1 ; 3. Generate ν 1 independent realizations {b i } 1≤i≤ν 1 of B from the model given by Eq. (3.13);

4. Add these realizations to the simulated response to obtain realizations of the train dynamical response y num,i = h(.; x mes,i , w 1 ) + b i , 1 ≤ i ≤ ν 1 ;

(3.17 For the validation of the methodology using numerical experiments, several tests have been performed on several reference values w 1 . The identification is performed for seven mechanical parameters:

(a) the vertical stiffness of the airsprings.

(b) the vertical stiffness of the airsprings located not at the same car junction as the sensors but at the other end of the car.

(c) the vertical damping of the airsprings.

(d) the torsion stiffness of the primary axial joints.

and the multiplicative coefficient applied on the velocity-effort characteristic of three types of dampers:

(e) the vertical primary dampers.

(f) the yaw dampers.

(g) the upper inter-carbody dampers.

All these parameters correspond to suspension elements located at the car junction where the sensors are located, except the second one (b). This parameter has been added to highlight the consequences of the sensitivity of the train dynamical response on the accuracy of the parameter identification.

GP surrogate model characteristics

The parameter identification is performed on seven parameters. Consequently, random vector W is of dimension q = 7, as well as the admissible set C W on which random surrogate model L(.; Θ) is indexed. A training set of 500 points is used to build the random surrogate model. The correlation function is Matérn-5 2 . It was chosen for the smoothness of the resulting surrogate model. For the reasons detailed in Section 2.2.2, the observations are considered noisy, except for the test of the EGO algorithm, which is adapted to exact observations. The nugget value is determined by optimization along with the other covariance parameters.

In the noisy case, we tried to refine the surrogate model using the KGCP approach as suggested in Section 2.2.3. However, it did not significantly improve the calibration results. Eventually, only the surrogate model build from the initial training set is used. The refinement could not increase the accuracy of the surrogate model probably because the accuracy improvement is small compared to the added nugget effect. It may be an indication that the nugget is too big, and that a smaller value should be used. The additional observation points proposed by the refinement procedure are mainly gathered in the vicinity of the likelihood maximum. An alternative solution could then be the introduction of a nugget with an variance varying from one observation to another. This variance could then be smaller in the vicinity of the likelihood maximum. The higher density of points in this region would then have an actual impact on the estimated accuracy of the surrogate model. This would however make the construction more complex.

Test of EGO algorithm

The optimization approach for the parameter identification presented in 2.3.2 is tested on the numerical experiments. The EGO algorithm [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF]] is used to maximize the likelihood function in order to obtain the most probable posterior values of the parameters of interest. The EGO algorithm relies on a GP surrogate model of the log-likelihood function, with observations that are supposed to be exact. The optimal parameters are searched among the parameters values for which the computational model has actually been run. Figure 3.9 presents the evolution of the known maximum of the log-likelihood and of the Euclidian distance between the values of the parameters maximizing the loglikelihood and the reference parameters values with iterations. More iterations have actually been run than the 120 shown on the graph, but with no improvement. One can observe that after a few iterations, the algorithm tends to be stuck on the same maximum for numerous iterations before finding a new one. This can be explained by the exploratory feature of the EI criterion. The first iterations provide the greatest improvements of the maximum because of the sparsity of the initial training set. There is no reason for any point of the initial training set to be close to the reference parameters values. Consequently, adding a single point in the vicinity of the maximum thanks to the EI criterion improves significantly the known maximum value.

Figure 3.10 presents, for each parameter, the evolution of the value maximizing the likelihood function, compared to the reference value. One can observe the convergence towards the reference for all parameters. The optimal parameters that maximize the likelihood function, resulting from the EGO algorithm, constitute a good approximation of the reference parameters. Figure 3.11: Results of the numerical experiment with nominal parameters. The parameters values are also normalized so that every parameter varies between 0 and 1; 0.5 then corresponds to the nominal value. The curves represent the normalized marginal posterior PDF for each parameter, obtained using the KP (red or dark gray) and MCT (orange or light gray) methods.

Validation of KP and MCT methods

Numerical experiments are used as a way of validating the novel KP and MCT methods for Bayesian calibration, as well as a way of highlighting the differences in the results they provide. The validation procedure has been performed with different values of the reference parameter w 1 . Figure 3.11 presents calibration results with w 1 fixed to the nominal values (w 1 = w 0 ). On this graph are displayed the marginal densities of the posterior PDF of W post 1

. The figure compares the results obtained with the KP method that solely uses the Kriging predictor provided by the surrogate model of the likelihood function and with the MCT method that also includes the surrogate-model uncertainty.

One can first observe that the distributions are close to the nominal values (corresponding to 0.5 on the graph). Except for the second parameter (b), the difference between the maximum of the marginal PDFs and the nominal value is always lower than 5% of the size of the admissible interval. The dispersion varies from one parameter to another. This can be explained by the initial choice of the admissible intervals, since the results are scaled according to their size, but more importantly by the sensitivity of the train dynamical response to the different mechanical parameters.

The results for parameter (b) are a good illustration of this sensitivity question. Contrary to the other parameters, it corresponds to suspension elements that are not located at the same junction as the sensors, but at the other end of the car. Consequently, we expect this parameter to have less influence on the train dynamical response measured by the sensors. This is coherent with the results: the distribution center is further away from the nominal value, and its dispersion is greater than for the other parameters. These considerations also highlight the interest of Bayesian calibration: because we measure the uncertainty on the calibrated parameters, we have a way of assessing the accuracy of the calibration. A parameter with a lower dispersion of its posterior marginal density is assumed to be more precisely identified.

The fact that the marginal PDFs are rather peaked comes from the large database that is available for the calibration. The comparison of the results obtained using the KP and MCT methods shows that the marginal PDF have a larger dispersion with the MCT method. With the KP method, a source of uncertainty, the error introduced by the approximation of the likelihood function by a surrogate model, is ignored. As a consequence, the uncertainty on the calibrated parameters is reduced. Using the KP method thus leads to an overestimation of the calibration accuracy. Nevertheless, the marginals maxima appear to be located at very similar parameter values with the two methods. The KP method seems to provide satisfying results if only the most probable parameter values is to be determined. One can also assume that the dispersion of the posterior PDF obtained with the KP method is similar to what one would obtain using the actual likelihood function, and not the GP surrogate model. From there, the KP method may be seen as a way to estimate the uncertainty on the calibrated parameters stemming only from the model uncertainties of the computational model and the measurements errors.

Figure 3.12 presents calibration results for arbitrary values of the reference parameter w 1 , representing degraded suspension elements. Only results obtained using the MCT method are displayed in this figure. One can note in the bottom graph that we have tested the case for which a parameter is fixed in the boundary of the admissible interval.

In these graphs, instead of the marginal posterior PDF, only the mean values of the calibrated parameters and the 98% confidence intervals are represented. These intervals are } as an estimate of w 1 can be discussed. Instead, one could choose the value maximizing the joint posterior PDF, or the vector gathering values maximizing the marginals corresponding to each parameter. Indeed, these latter better correspond to the "most probable values" of the parameters. In the present case, we observed that there is no significant difference between all these estimates, in particular because the asymmetry of the posterior PDF is small.

The various examples presented in this section show that the numerical experiments are conclusive for the validation of both the KP and the MCT methods. Limitations are expressed for the KP method, as far as the the estimation of the uncertainty on the calibrated parameters is concerned. The interest of the Bayesian approach is highlighted as well as the impact of the sensitivity of train dynamical response to the mechanical parameters. In the following section, the KP and MCT are applied to real measurements of the train dynamical response.

Parameter identification from actual measurements

Measurements of the train dynamical response are available at several different dates. These latter are denoted (T k ) k with k = -5, . . . , 4 and chronologically ordered. The time steps between dates are of several months. Date T 0 corresponds to the reference date, used to identify the output-predictive error B. As explained in Section 3.5, a major maintenance operation of the train took place between T -1 and T 0 , just before T 0 , during which the suspension elements have been renewed. Consequently, the parameters are supposed to be nominal at date T 0 . For all dates, the available data correspond to the same line. They gather measurements of the train dynamical response on the same track stretches.

This section presents the results obtained by applying the calibration method on real measurements of the train dynamical response. Section 3.7.1 analyzes the difference between the results obtained on real measurements and those obtained on numerical experiments, with the KP and MCT methods. Section 3.7.2 is devoted to the long-term evolution of the mechanical parameters of the suspension elements obtained by applying the MCT method successively at all date (T k ) k .

Differences with the numerical experiments

Figure 3.13 presents the calibration results obtained using the KP and MCT methods using actual acceleration measurements (at date T 1 ). The same seven parameters, as listed in Section 3.6.1, are used. The comparison between the two methods is discussed in Section 3.6.4. No additional comment on the matter is made from the results presented here.

However, one significant difference should be noted compared to Figure 3.11. The posterior PDF dispersion is larger for the calibration performed on real measurements than on numerical experiments. This larger dispersion is probably caused by the choice made for the construction of output-predictive error B. In the case of the numerical experiments, B exactly represents the distance between the measured and simulated train response, since the simulated experimental responses {y num,i } 1≤i≤ν 1 have been generated by adding independent realizations of B to simulated responses. When actual measurements are used, the model of B is not able to exactly represent the distance between the measured and simulated train response anymore, which results in a lower accuracy of the calibration.

Long-term evolution of the suspension parameters

In this section, the MCT method of Bayesian calibration is used to study the long-term evolution of the suspension mechanical parameters. To this end, the data available at all ten dates T k , -5 ≤ k ≤ 4 are used. They represent almost five years of measurements. The maintenance operation that happened between T -1 and T 0 constitutes however a discontinuity in the suspension elements degradation. The data before maintenance (-5 ≤ k ≤ -1) is exploited to compare the calibration results resulting from old, potentially degraded elements, and new ones.

The parameter (b) has been removed from the present study. Indeed, Figure 3.13 shows that the support of the corresponding marginal PDF covers more than half of the admissible interval of the parameter. Consequently, we conclude that trying to identify this parameter appears too inaccurate to be relevant. Removing this parameter is not an issue, since its purpose was to highlight the impact of the sensitivity of the train dynamical response on the parameter identification.

Figures 3.14 and 3.15 present the long-time evolution of the six suspension parameters of interest as identified using the MCT method from the measurements performed at the ten dates (T k ) k . It has to be noted that for the dates before maintenance, (T -5 , . . . T -1 ), The left graphs present the time evolution of the mean value of the calibrated parameter, surrounded by the 98% confidence interval. The red diamonds indicate the dates before maintenance (represented by the vertical line), the blue dots the dates after. The scale is normalized, meaning that the (initial) admissible interval has been shrunk to [0, 1] for each parameter (0.5 is then the nominal value). The right graphs present the time evolution of the standard deviation of the calibrated parameter, on the same normalized scale. it was necessary to enlarge the admissible intervals of certain parameters. If not, then the identified parameter values would remain stuck in the boundary. However, in the graphs, the normalization remains unchanged: the interval [0, 1] still corresponds to the initial admissible intervals.

We can first observe that the mean values of the calibrated parameters at reference date T 0 (first blue dot) are close to the nominal values. This is in accordance to the assumption made for identifying the output-predictive error B, and constitutes a second type of validation of the method.

Important differences can be observed between the results before and after maintenance. The calibrated parameters are further from the nominal values before maintenance, except for the yaw and inter-carbody dampers, for which the values are similar before and after maintenance. It corresponds to the fact that the suspension elements are more degraded before the maintenance operation. Each confidence interval (whose size is proportional to the standard deviation of the calibrated parameter) is much wider before maintenance. An explanation is the fact that B represents better the distance between the measured response and the simulated response when the parameters are close to their nominal values. This highlights the limit of the hypothesis that B does not depend on the parameters values.

The values identified for the primary vertical dampers, on graph (e), seem unrealistic before maintenance. When degraded, the dampers usually witness a decrease of their mechanical characteristics, while the identification shows a large increase with respect to the nominal value. The analysis of the measured train dynamical response suggests that nonlinearities affect the primary suspension, but are not accounted for in the train model. Figure 3.16 presents the average amplitude of the vertical bogie acceleration in the frequency domain before and after maintenance. We can observe that before maintenance, so in the degraded state, a second mode seems to appear, very flat, which could come from nonlinear stiffnesses in the primary suspension, not taken into account. Recently, additional simulations have been performed, which advocate for this hypothesis. The identified parameters values for the primary vertical dampers may result from an attempt of the calibration procedure to compensate an evolution of these nonlinearities.

We can observe that the long-time evolution of the parameters as identified by the MCT method is generally nonmonotonic, which is unexpected. As previously stated, more complex models for the suspension elements would allow for more realistic mechanical behaviors. The possible degradation patterns can also be improved. For instance, experimental tests on degraded dampers have shown that they may display asymmetric or hysteretic behaviors, as well as characteristics depending on the displacement and not only the velocity. A more complex train model is expected to improve calibration, but would increase the problem dimension and would require more intensive computations.

Conclusion on the industrial application

In this chapter, we have introduced the computational model on which the proposed suspension-monitoring solution is based, with a particular focus on the suspension elements of interests and on the input and output physical quantities of the model. We detailed the identification procedure adopted for the output-predictive error, through which the uncertainties of the system are accounted for.

Both the KP and MCT Bayesian calibration methods have been validated on numerical experiments, which in addition have highlighted the added value of the MCT method. The latter allows for taking into account the surrogate-model uncertainty. The good results also validated the choices made concerning the construction of the GP surrogate model. The application on actual measurements of the train dynamical response provided encouraging results, as they showed a clear gap between the suspension calibrated parameters before and after maintenance. However, no clear evolution of the mechanical parameters due to gradual degradation could be pointed out. The degradation predictions could be improved thanks a greater number of sensors for the measurements of the train dynamical response, and, probably, thanks to an improvement of the computational model.

Chapter 4 Conclusion and perspectives

In this thesis, we have developed a Bayesian calibration method for a computational model with functional input and output, for which one evaluation is expensive. The purpose is the statistical inverse identification of the parameters of a mechanical system, taking into account the various uncertainties it may be affected by. The motivation for this work is the development of a monitoring solution for high-speed train suspensions for maintenance purposes, using embedded accelerometers. This industrial problem has been used to validate the new calibration method proposed. Nevertheless, we are convinced of the relevancy of the method for applications to other systems and in other industrial fields.

The Bayesian calibration of a computational model consists in estimating the probability distribution of the parameters of the model from experimental data bringing information about the model outputs. In the framework of Bayesian calibration, the methods existing in the literature were not sufficient for the specific problem we are dealing with, due to the functional nature of the system input and the system output, and due to the numerical cost of the computational model. This is the reason why we have proposed a novel Bayesian calibration method. In particular, we have not used the classical representation of the computational model by a GP surrogate model. Instead, we have chosen an innovative approach consisting in using a random surrogate model to represent the loglikelihood function that is at the core of the Bayesian approach. A surrogate model of the likelihood function circumvents the issues related to the functional input and output and decreases the computational cost of the calibration procedure. The novel method fully exploits the GP surrogate model as it allows for taking into account the surrogate-model uncertainty. The latter quantifies the error induced by the approximation of any target function by a surrogate model. Because uncertainties have a significant influence on the results of the calibration, it is important to consider the surrogate-model uncertainty in the calibration procedure. Implementing the proposed method is not self-evident and requires the use of advanced computational statistics tools such that: the identification of the stochastic model of the output-predictive error, the construction of the GP surrogate model, the calculation of trajectories of the GP process, the estimation of the posterior PDF thanks to the TMCMC algorithm.

The validation of this novel Bayesian calibration method has been performed using several numerical experiments carried out in the framework of an industrial application. These tests have demonstrated the efficiency of the proposed procedure for calibrating parameters of a complex stochastic dynamical system.

One perspective for this work would be to replace the Gaussian-process surrogate model of the log-likelihood function by another non-Gaussian surrogate model, using either a direct construction of a parameterized non-Gaussian measure for which a generator has to be known or an indirect construction such as Polynomial Chaos representation for which the generator is known. With such constructions, all the procedure that has been proposed would be preserved and could be used.

A second perspective of this work would be the application and the extension of the proposed Bayesian calibration procedure to the case for which the vector of parameters is itself a non-Gaussian stochastic process indexed by time, for which there are significant time correlations/dependencies.

The industrial motivation of the present work was the development of a health-state monitoring method for high-speed train suspensions from acceleration measurements by embedded sensors. Such monitoring solution would allow for optimizing the trainmaintenance cycles by replacing the suspension elements at the most appropriate time. It could then allow for decreasing maintenance costs by reducing the number of manual inspections and by limiting the replacement of healthy and functional suspension elements.

Embedded accelerometers provide only an indirect information about suspension health state. Mathematical processing are necessary to extract the relevant information from the measured data. This is the aspect of the solution we focused on. The literature provides various solutions to this monitoring problem in a railway context. Compared to the already existing methods, the major contribution of this thesis lies in the fact that the proposed method is suited for the use of a complex nonlinear computational model, while including the uncertainties of the system. We have developed a solution based on a multibody-dynamics simulation software. It allows for modeling an entire train rather than single car. Moreover, it allows for injecting real measured track irregularities, which have a strong influence on the train dynamical behavior. Taking uncertainties into account is necessary to perform a robust monitoring. Moreover, a large database of measurements of the train dynamical response and of the corresponding measurements of the track-geometry irregularities was available, on which the method has been tested.

The proposed Bayesian calibration method has been validated on numerical experiments. The results have highlighted the relevancy of the Bayesian approach. Indeed, it allows for estimating confidence intervals around the calibrated parameter values. They are a way of assessing the accuracy of the calibration, but more importantly provide precious information for decision making. A fully experimental validation campaign is expected. It would involve experimental testing of the isolated suspensions elements in order to determine their real mechanical characteristics and consequently health state. Measurements of the train dynamic response of a train equipped with the exact same elements should then be performed. The calibrated values of the mechanical parameters could then be compared to their "real" experimental value in order to assess the actual efficiency of the full method for industrial application.

The proposed Bayesian calibration method has also been tested on actual measurements of the train dynamical response, in order to study the time evolution of the suspension mechanical parameters. The results obtained are encouraging, but also underline the crucial importance of the suspension-elements modeling and of the choice of parameterization. They control the complexity of the nominal mechanical behavior of the suspension elements, but also the possible degradation they may undergo. As a matter of fact, some elements such as dampers tend to display much more complex behavior when degraded. The parameterization has been kept rather simple, in order to limit the dimensionality of the problem. Some results appear to be affected by a lack in the modeling of certain nonlinearities of the mechanical system. More generally, the quality of the monitoring method strongly depends on the quality of the computational model, which must account accurately for the influence of the suspension parameters on the train dynamical behavior.

A next important step is the development of a criterion for the trigger of maintenance operations, based on the posterior PDF of the suspension parameters provided by the calibration method.

Publications and external communications

The work achieved in this thesis has given rise to a series of communications, which are listed hereinafter.

• GP(m; Σ) denotes the Gaussian process whose mean function is x → m(x) from X into R and covariance function is (x, x ) → Σ(x, x ) from X 2 into R;

• f : X → R p gathers p deterministic regression functions on X .

• β ∈ R p is the vector of regression coefficients parameterizing the GP mean function.

• It is assumed that the correlation function R(x, x | ψ) depends only on xx and the variance σ 2 of Y (x) is independent of x. Consequently, the covariance function is written as Σ(x, x ) = σ 2 R(x -x | ψ), and R(0 | ψ) = 1. It depends on a vector ψ of parameters (for example, correlation lengths).

The mean function (that depends on x) is decomposed as a linear combination of regression functions chosen by the user, typically polynomials. The first parameter is thus the vector of regression coefficients β. The second parameter is the GP variance σ. The correlation function R(. | ψ) is also chosen by the user according to the expected regularity of y, for instance an exponential or a Matérn. The third parameter is the vector of correlation parameters ψ. They depend on the correlation function. Parameters (β, σ, ψ) are a priori unknown. They are not chosen by the user but have to be identified using the information about y provided by the observations on the training set. The following sections detail how to deal with the fact that parameters (β, σ, ψ) are actually unknown. They are identified using the information provided by the training set. relevancy for the representation of vector B (taking them into account had a strong negative impact in the estimation of the PDF of B). In addition, we set the truncation size so that covariance matrix [C γ ] is invertible.

Coefficient α is determined by minimizing the Leave-One-Out error of projection (see [START_REF] Perrin | A Posteriori error and optimal reduced basis for stochastic processes defined by a finite set of realizations[END_REF]) of the realizations {b ref,i } i on the basis [Φ ( ) ] obtained from the spectral decomposition of [ C] that depends on α.

Figure 2

 2 Figure 2.1: Diagram of the Bayesian calibration method

Figure 2

 2 Figure 2.2: 1-D examples of GP surrogate models. The black solid curve represents the target function, the dotted green curve the GP mean function (or Kriging predictor), the green area the 95% quantiles, the red diamonds the observations. Observations are exact on the left graph and noisy on the right one.

Component x 1

 1 Figure 2.3: Smoothing of the surrogate model thanks to the addition of a nugget effect. Both graphs represent a slice of the GP surrogate model along the first dimension of set X . The dotted curve corresponds to the mean function and the surrounding colored area to the 95% quantiles. The observations are supposed exact on the left graph and noisy on the right one. They are not represented because none belongs to this particular slice of multidimensional set X .

Figure 2 . 4 :

 24 Figure 2.4: Diagram of the MCT method

Figure 3 . 2 :

 32 Figure 3.2: TGV bogie with the suspension elements of interest: (a) vertical primary damper, (b) elastomeric stiffness, (c) yaw damper, (d) airspring

  Figure 3.3: Rheological model for the axial joint

  Figure 3.5: Track irregularities example. From bottom to top: lateral, gauge, vertical, cross-level irregularities.

  Figure 3.6: Vertical (top curve) and lateral (bottom curve) accelerations in a carbody, in the time domain.

FrequencyFigure 3

 3 Figure3.8: For several stretches, statistical representation of the experimental lateral acceleration in a carbody in the frequency domain. The solid line represents the mean function, while the filled area is bounded by the 5% and 95% quantiles.

  mapping h : ({x(s) , s ∈ [0 , S]}; w) → { a(ω) , ω ∈ Ω} .

) 5 .; 6 .

 56 Perform the calibration using input data{(x mes,i , y num,i )} 1≤i≤ν 1 (3.18)to obtain the calibrated random vector W post 1 Compare statistics of W post 1 to w 1 .

  Figure 3.9: Test of EGO algorithm: Evolution with the algorithm iterations of the known maximum of the likelihood function on the training set (solid blue curve) and of the Euclidian distance between the values of the parameters maximizing the likelihood function and the reference parameters values (dotted red curve).

  Figure 3.12: Results of the numerical experiments with arbitrary reference parameters, using the MCT method. The reference parameters values w 1 (black triangles) are compared to the mean of the marginal posterior PDF (blue dots). The blue lines represents the 98% confidence intervals around these calibrated values. The parameter scale normalization is identical to figure 3.11.

  Figure 3.13: Comparison of calibration results obtained at time T 1 using the KP (red or dark gray curves) and the MCT (orange or light gray curves) methods. The graph layout is identical to figure 3.11.

  Figure 3.14: Parameter identification at ten dates T k for three mechanical parameters: (a) Airspring stiffness, (c) Airspring damping coefficient, (d) Primary elastomeric stiffness.The left graphs present the time evolution of the mean value of the calibrated parameter, surrounded by the 98% confidence interval. The red diamonds indicate the dates before maintenance (represented by the vertical line), the blue dots the dates after. The scale is normalized, meaning that the (initial) admissible interval has been shrunk to [0, 1] for each parameter (0.5 is then the nominal value). The right graphs present the time evolution of the standard deviation of the calibrated parameter, on the same normalized scale.

Frequency

  Figure 3.15: Parameter identification at ten dates T k for three mechanical parameters: (e) Primary vertical damper, (f) Yaw damper, (g) Upper inter-carbody damper. The axes layout is identical to Figure 3.14.

A. 2

 2 Conditioning by the observationsLet y n and Y n be respectively the values of function y and of process Y at points {x obs i } 1≤i≤n : any x fixed in X , the vector concatenatingY (x) | {β, σ, ψ} and Y n | {β, σ, ψ} is Gaussian: Y (x) Y n | {β, σ, ψ} ∼ N f(x) t [F ] β ; σ 2 1 r(x) t r(x) ij = R(x obs i -x obs j | ψ), 1 ≤ i, j ≤ n. Matrix [R]is assumed to be invertible. The surrogate model is then obtained by conditioning process Y by Y n = y n . When doing so, Y | {Y n , β, σ, ψ} remains Gaussian. Its mean function and covariance function are immediately deduced from the formula of the conditioned Gaussian random variables:Y | {Y n = y n , β, σ, ψ} ∼ GP x → f(x) t β + r(x) t [R] -1 (y n -[F ]β) ; (x, x ) → σ 2 (R(x -x | ψ) -r(x) t [R] -1 r(x )) . (A.5)

  obs,k )

		Training set {w obs,k } k
	Experimental data
	Inputs {x mes,i } i Outputs {y mes,i } i	Computational model: run h(x mes,i ; w obs,k )
			Error B
		GP model L(.; Θ)	Subset X *
		Refinement
			Conditioning
			set X c
		Trajectories L(.; θ j )
		on X Estimate p post W

c

Approximates L(.; θ j ) TMCMC on L(.; θ j )

Acknowledgments

B Alternative construction of the covariance matrix of B

Set j = 1 Draw samples {w 0,k } 1≤k≤Ns distributed according to p prior W while q j-1 < 1 do Determine exponent q j from Eq. (2.45) for k = 1 : N s do Compute ᾱj,k from Eq. (2.46) end for Compute the covariance matrix [Σ j ] from Eq. ( 2.47) Initialize {w c k } 1≤k≤Ns such that ∀k, w c k = w j-1,k for k = 1 : N s do Draw index ∈ {1 . . . N s } with a probability ᾱj, Draw proposal w * from a normal distribution N (w c ; [Σ j ]) Draw r from a uniform distribution on [0, 1]; if r ≤ p j (w * ) p j (w c ) then w c = w * end if Set w j,k = w c end for Set j = j + 1 end while

Summary of the MCT calibration method

The complete calibration method we proposed articulates several techniques. For clarity, we present a short summary of the method in this last section, together with the diagram in Figure 2.4. The MCT method may be divided into the following steps:

1. Define the training set {w obs,k } k , space filling in the admissible set C W of the parameters to be identified;

2. For each parameters values of the training set, run the computational model on the available measurements of the input excitation {x mes,i } i , that is to say run h(x mes,i ; w obs,k ) for all i and k;

3. Using the chosen model for error B and the available measurements of the output response {y mes,i } i , compute the log-likelihood function values on the training set according to Eq. (2.16);

4. From the observations (w obs,k , L(w obs,k )) k , build the GP surrogate model L(.; Θ) of the log-likelihood function, as detailed in Appendix A. If necessary, perform a refinement step as described in Section 2.2.3; 5. Draw the conditioning set X c in the vicinity of the maximum of the mean value of the GP surrogate model, using TMCMC with the indicator function on the subset X * , deduced from the definition of the latter given by Eq. ( 2.39);

6. Draw trajectories of the GP surrogate model on X c and build the approximates L(.; θ j ) of these trajectories by further conditioning the GP surrogate model, as explained in Section 2.2.4;

We first consider a partial definition of B as a second-order random vector, of dimension nµ, which admits m as mean vector and [ C] as covariance matrix (at this point, the probability distribution of B is not defined yet). The spectral decomposition of [ C] yields (3.12) where [Φ] is a rectangular orthonormal matrix (such that [Φ] t [Φ] = [I r ]), of dimension nµ × r, and contains, in columns, the r eigenvectors of [ C] associated with the strictly positive eigenvalues that constitute the diagonal of matrix [Λ].

The use of the Principal Component Analysis (PCA) allows random vector B to be written as

where ξ is a second-order centered random vector, of dimension r, which admits the identity matrix as covariance matrix.

From the Maximum Entropy principle of Information Theory, under the constraints defined by the above available information about ξ, it can be deduced that ξ is a Gaussian random vector (and consequently, not degenerated). Since any affine transformation of a Gaussian vector is a Gaussian vector, from Eq. (3.13), it can be concluded that B is Gaussian as well. However, because [ C] is not invertible, Gaussian vector B does not admit a probability density function. Equation (3.13) shows that Bm is with values in the vectorial subspace of R nµ spanned by the columns of [Φ]. For the calculation of the log-likelihood, as given by Eq. ( 2.16), we thus choose to replace (this is a modeling) the PDF p B by the pseudo-PDF p B defined as

where

• p ξ is the PDF of random vector ξ, so the canonical multivariate Gaussian density of dimension r, written as

• ϕ is the projection on the vectorial subspace to which the values of Bm belong, defined as

Function p B is called a pseudo-PDF because it is not integrable on R nµ with respect to the Lebesgue measure dy in R nµ . This pseudo-PDF has to be viewed as a formal writing that means that for any realization y j of B, the corresponding realization of ξ is given by

It should be noted that another model has been studied consisting in introducing a regularization of covariance matrix [ C] in order to make it invertible. Such a model has not given satisfactory results.

For the numerical applications the construction of the vector basis [Φ] has been improved using the methodology presented in [START_REF] Perrin | A Posteriori error and optimal reduced basis for stochastic processes defined by a finite set of realizations[END_REF] and detailed in Appendix B. Appendix A

Conferences without proceedings

Gaussian-process surrogate modeling

Gaussian-process surrogate modeling consists in representing a given deterministic scalar target function by a conditioned Gaussian process (GP). Information about the target function is provided by a set of observations of the function value in some points of the definition set. These points constitute the initial training set. Randomness is introduced as a way to quantify the approximation stemming from the use of a surrogate model. The objective of GP surrogate modeling is, for example, to build an efficient approximation of a computational model for which one evaluation is expensive.

The evaluation of the target function in the points of the training set is generally the computationally expensive step of GP surrogate modeling. Consequently, this set should contain a limited number of points while maximizing the information provided. Therefore, a common choice is to define a space-filling training set to get as much information as possible about the global behavior of the target function in its whole definition set. Appendix A.6 provides an example on how to build a space-filling training set. After this initial training phase, a second refining step can be performed in order to explore particular features of the target function (see Section 2.2.3).

In this appendix, we present the principle of GP surrogate modeling based on the Bayesian approach (see for instance [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]). Note that another equivalent formulation, based on the minimization of the mean-square error, could have been presented (see for instance [START_REF] Sacks | Design and Experiments of Computer Experiments[END_REF]).

A.1 Problem statement

The goal is to build a surrogate model of a given deterministic real-valued function y defined on a subset X of R q , whose value is known only in the n training points {x obs i } 1≤i≤n in X .

Function y is supposed to be a particular sample path of an underlying Gaussian second-order real-valued stochastic process {Y (x), x ∈ X }, indexed by X . The surrogate model consists of this GP conditioned by the observations of y in the points of the training set. The underlying process is defined according to a parametric formulation, and denoted

where:

A.3 Mean function parameter β

In this section, σ and ψ are supposed to be known. Only the regression coefficients β are supposed to be unknown. The fact that no information is a priori available about β is taken into account by following a hierarchical approach. Parameter β is represented by a random vector with non-informative prior: p β ∝ 1. The training set is then used to learn about the distribution of β. The principle is to determine the distribution of β | Y n = y n for a fixed value of σ and ψ. In this section, the conditioning on σ, ψ is not systematically repeated for simplicity. Using the Bayes formula, one can write:

Knowing the probability distribution of β | Y n for a fixed value of σ and ψ, process Y | {Y n , β, σ, ψ} can be statistically averaged with respect to random vector β. It means that Y | {Y n , σ, ψ} can be used instead of Y | {Y n , β, σ, ψ}. This is achieved by relying on the rule of conditional expectation. For x in X , the mean value and variance of Y (x) | {Y n = y n , σ, ψ} are:

. This can be written as

A.4 Variance and correlation parameters

Parameters σ and ψ could be estimated using the same hierarchical approach, by putting prior distributions on these variables. However, in this case, no closed form can be determined in general for Y . Instead, parameters σ and ψ are determined according to a criterion assessing how well the Gaussian process is fitting the data provided by the training set. The criterion usually used is the Maximum Likelihood Estimation (MLE). The principle of MLE criterion is to maximize the density

The value of σ maximizing this density can be determined explicitly:

With this value of σ, the previous maximization is equivalent to minimizing the quantity

1 n in order to determine the optimal value of ψ. In general, no closed form exists for the covariance parameter; this step has to be performed numerically.

A.5 Case with noisy observations

Building a surrogate model based on a conditioned Gaussian process is also possible when the values of y on the training set are not computed exactly but affected by a random noise. A real-valued Gaussian random variable (noise) of zero mean and variance σ 2 is independently added to each component of Y n such that Y n = Y n + ε. In terms of components, for 1 ≤ i ≤ n, we have Y n i = Y n i + ε i in which ε 1 . . . , ε n are n independent copies of . Random vector ε is a Gaussian centered vector with covariance matrix σ 2 [I n ], with [I n ] the identity matrix of dimension n. Random vector ε is assumed to be independent of Y n .

In such a case, the available data are not y n as previously defined, but rather y n that gathers the observations of y in each point x obs i of the training set plus a realization e i of the random noise :

Stochastic process Y and random vector β must now be conditioned by Y n = y n and not by Y n = y n . For x ∈ X , the joint probability distribution of Y (x) and Y n can be expressed as follows:

The results given in Sections A.2 and A.3 hold, with correlation matrix

No closed form can be found anymore for the optimal value of σ using the MLE criterion. Consequently, it has to be optimized numerically along with parameter ψ. The variance σ 2 of the noise can be set by the user or optimized along with parameters σ and ψ.

A.6 Space-filling training set

This appendix presents one method for the generation of a space filling training set in a set X consisting of the set product of bounded intervals of R. The chosen method is the optimization of a Latin Hypercube Sample (LHS) according to a maximin distance [START_REF] Chen | A review of design and modeling in computer experiments[END_REF].

To build a N -points LHS, admissible set X must be divided into N cells along each dimension. The points are then scattered in X so that in every dimension, each cell only contains one sample point. The interest of a LHS is that if you consider a particular dimension, the sample points are regularly spaced, with no redundancy.

However, a LHS is not necessarily space-filling. We choose to measure the spacefilling property with the distance δ that consists of the smaller Euclidian distance between two points of the training set: for a finite discrete subset W of X , distance δ is defined as

The greater δ(W), the most space-filling W is. In practice, numerous LHS candidates are drawn. The one for which the value of criterion δ is the greatest is kept as the best training set in X .

Appendix B

Alternative construction of the covariance matrix of B

The PCA representation given by Eq. (3.13) using the empirical estimate [ C] of the covariance matrix can be improved by using the following methodology proposed by [START_REF] Perrin | A Posteriori error and optimal reduced basis for stochastic processes defined by a finite set of realizations[END_REF]. When the number of realizations is smaller than the dimension of the random vector to be identified, the projection basis The discretization of the restriction of a mean-square stationary process, indexed by R, to a bounded interval of R, admits a covariance matrix with constant diagonals. We denote as J µ the stationarization function defined on the set of all µ × µ real matrices:

where the entries of matrix [M stat ] are defined as

Process B has n components, concatenated in the discretized version of the process as shown in Eq. (3.7). Consequently, covariance matrix [ C] is a block matrix:

The µ × µ submatrices [C (kk ) ], (k, k ) ∈ {1, . . . , n} 2 , correspond to the covariance matrices of each component (when k = k ) and to the correlation matrices between two components (when k = k ). One should be careful to apply the stationarization independently on each submatrix [C (kk ) ]. Matrix [ C stat ] can then be written as where p ξ is the canonical multivariate Gaussian density of dimension .

The size of truncation of the basis is usually defined according to a convergence tolerance expressed as a function the eigenvalues, typically 10 -6 . In the present case, we had to set a greater tolerance. The eigenvectors of small eigenvalues have no statistical