E Margo Molhoek

Albert Van Dijk

Edwin J A Veldhuizen

Henk P Haagsman
email: h.p.haagsman@uu.nl

Floris J Bikker

A cathelicidin-2-derived peptide effectively impairs Staphylococcus epidermidis biofilms

Keywords:

Staphylococcus epidermidis is a major cause of nosocomial infections owing to its ability to form biofilms on the surface of medical devices. Biofilms are surface-adhered bacterial communities. In mature biofilms these communities are encased in an extracellular matrix composed of bacterial polysaccharides, proteins and DNA. The antibiotic resistance of bacteria present in biofilms can be up to 1000-fold higher compared with the planktonic phenotype. Host defence peptides (HDPs) are considered to be excellent candidates for the development of novel antibiotics. Recently, we demonstrated that a short variant of the HDP chicken cathelicidin-2, peptide F 2,5,12 W, has potent antibacterial and lipopolysaccharide-neutralising activities. This study reports on the antibiofilm activity of peptide F 2,5,12 W against two strains of S. epidermidis, including a multiresistant strain. Peptide F 2,5,12 W potently inhibited the formation of bacterial biofilms in vitro at a low concentration of 2.5 M, which is below the concentration required to kill or inhibit growth (minimal inhibitory concentration = 10 M).

Moreover, peptide F 2,5,12 W also impaired existing S. epidermidis biofilms. A 4-h challenge of pre-grown biofilms with 40 M F 2,5,12 W reduced the metabolic activity of the wild-type strain biofilm completely and reduced that of the multiresistant strain biofilm by >50%. It is concluded that F 2,5,12 W prevents biofilm formation and impairs mature S. epidermidis biofilms.

A c c e p t e d M a n u s c r i p t

Introduction

Staphylococcus epidermidis is a major cause of infections related to indwelling medical devices such as intravascular catheters, cerebrospinal fluid shunts, peritoneal dialysis catheters, intraocular lenses, cardiac pacemakers and prosthetic joints [START_REF] Raad | Infectious complications of indwelling vascular catheters[END_REF]. It is an important member of the human skin and mucous membrane microflora and can easily be transmitted to the surfaces of these devices when they are implanted or manipulated.

Staphylococcus epidermidis adheres to and grows on medical devices, which may lead to a mature biofilm consisting of cells embedded in a sticky extracellular slime composed of bacterial polysaccharides, proteins and DNA. This matrix renders the S. epidermidis biofilm resistant to antibiotics and host defence systems and makes it difficult or almost impossible to eradicate. Bacteria in biofilms can be up to 1000-fold more resistant to antibiotic treatment than the planktonic phenotype [START_REF] Gilbert | Biofilm susceptibility to antimicrobials[END_REF].

The rising number of infections caused by bacterial isolates resistant to conventional antibiotics has led to an intense search for novel antimicrobials and chemotherapeutics, including the use of host defence peptides (HDPs). HDPs comprise a class of antibiotic compounds, i.e. naturally occurring peptides, that exhibit broad antibacterial and/or immunomodulatory properties [START_REF] Zasloff | Antimicrobial peptides of multicellular organisms[END_REF].

Recently, the chicken HDP cathelicidin-2 (CATH-2) has been identified [START_REF] Van Dijk | CMAP27, a novel chicken cathelicidin-like antimicrobial protein[END_REF]. We demonstrated that CATH-2 displays potent antibacterial activity against a variety of bacteria, including biowarfare agents, and is able to inhibit lipopolysaccharide (LPS)induced cytokine release by peripheral blood mononuclear cells [START_REF] Van Dijk | Chicken heterophils are recruited to the site of Salmonella infection and release antibacterial mature cathelicidin-2 upon stimulation with LPS[END_REF]. Furthermore, we is the core element of CATH-2 required for antibacterial activity. The antibacterial and LPS-neutralising activities of C1-15 were subsequently improved by tryptophan substitution of all three phenylalanine residues, resulting in peptide F 2,5,12 W [START_REF] Van Dijk | Identification of chicken cathelicidin-2 core elements involved in antibacterial and immunomodulatory activities[END_REF][START_REF] Molhoek | Chicken cathelicidin-2-derived peptides with enhanced immunomodulatory and antibacterial activities against biological warfare agents[END_REF].

Here we report the inhibitory effect of F 2,5,12 W on S. epidermidis biofilm formation as well as on existing biofilms.

Material and methods

Peptide

F 2,5,12 W (RWGRWLRKIRRWRPK) was purchased from Pepscan Therapeutics (Lelystad, The Netherlands). Purity was determined using high-performance liquid chromatography (HPLC) and was >90%.

Micro-organisms

The strains used were the S. epidermidis clinical isolate BM185 and the S. epidermidis multiresistant strain BM492.

Bacterial susceptibility assay

Antimicrobial susceptibility of S. epidermidis was determined using the broth microdilution method. Briefly, two-fold serial dilutions of F 2,5,12 W were prepared in

Biofilm formation and quantification of activity against biofilms

Biofilm formation was induced in 96-well flat-bottom tissue culture plates (Costar ® ). First, overnight cultures of S. epidermidis strains were diluted 1:100 in fresh MHB containing 5% glucose and a 100 L aliquot of the diluted culture was added to each well. Following incubation at 35 C for 24 h without shaking, each well was washed twice with phosphate-buffered saline (PBS) to remove planktonic bacteria.

The biomass of the biofilm was quantified by semiquantitative crystal violet (CV) staining. After removal of planktonic bacteria, wells were stained with 0.4% CV for 15 min. Excess CV was removed by washing each well three times with PBS. Bound CV was solubilised in 100 L of 33% acetic acid and the optical density at 630 nm was measured.

A c c e p t e d M a n u s c r i p t 7 Metabolic activity of the biofilm was determined with the alamarBlue ® assay (Invitrogen, Carlsbad, CA). Following removal of planktonic bacteria, 100 L of MHB containing 5% glucose and 5% alamarBlue was added to each well. alamarBlue is a redox indicator that both shows fluorescence and changes colour in response to chemical reduction.

The extent of reduction is a reflection of bacterial viability [START_REF] Pettit | Microplate Alamar blue assay for Staphylococcus epidermidis biofilm susceptibility testing[END_REF]. Following 1 h incubation at 37 C, fluorescence was measured at excitation and emission wavelengths of 550 nm and 590 nm, respectively. The effects of F 2,5,12 W on biofilm formation were determined by addition of different concentrations of peptide (0-40 M) during the 24-h biofilm formation period. Cell growth was determined by measuring the absorbance at 595 nm.

Similarly, the effects of F 2,5,12 W on mature biofilms were determined by incubating washed 24-h-old biofilms with different concentrations of peptides (0-40 M) for 4 h.

Positive controls were untreated biofilms with only MHB plus 0.5% glucose. All experiments were performed in triplicate.

Statistical analysis

Differences between mean values of groups were evaluated using one-way analysis of variance (ANOVA) and Dunnett's post hoc test and were considered significant at P < 0.05.

A c c e p t e d M a n u s c r i p t 8

Results

Minimal inhibitory concentration

The concentration of F 2,5,12 W required to inhibit the growth of planktonic bacteria (MIC) was 10 M (data not shown). No differences in MIC values were found between the antibiotic-sensitive strain and the multiresistant strain. Both the clinical isolate and the multiresistant strain showed no visible growth at a concentration of 10 M.

Biofilm formation of two Staphylococcus epidermidis strains

Both strains formed biofilms in vitro. Culturing for 24 h resulted in a significantly larger biomass and greater metabolic activity for the multiresistant strain (BM492) compared with the clinical wild-type strain (BM185). Furthermore, a strong correlation was observed between biomass quantified by CV staining (Fig. 1A) and biofilm metabolic activity quantified by alamarBlue reduction (Fig. 1B) in the 24-h-old biofilm.

Peptide F 2,5,12 W impairs Staphylococcus epidermidis biofilm formation

An in vitro assay was performed to investigate the effects of F 2,5,12 W on biofilm formation of S. epidermidis. As shown in Fig. 2A,B, all tested concentrations of the peptide inhibited biofilm formation. No differences were found between the two strains tested.

The amount of biofilm biomass as well as the metabolic activity of the biofilm was significantly inhibited at 2.5 M peptide, a concentration four times below the MIC. This A c c e p t e d M a n u s c r i p t 9 inhibitory effect was increased to control levels at 10 M, a concentration equal to the MIC (Fig. 2).

F 2,5,12 W impairs the metabolic activity of existing biofilms

Mature biofilms are known to be more resistant to eradication compared with inhibition of de novo biofilm formation. Therefore, the activity of peptide F 2,5,12 W against mature biofilms was evaluated. Fig. 3 show the amount of biomass quantified by CV staining as well as metabolic activity quantified by the alamarBlue assay in treated and untreated biofilms.

Metabolic activity of S. epidermidis strain BM185 was significantly suppressed at 5 M peptide, a concentration two times below the MIC (Fig. 3C). Higher concentrations further suppressed the metabolic activity; at 40 M the metabolic activity was similar to the negative control levels. The effect of F 2,5,12 W treatment on the biomass of strain BM185 was also dose-dependent (Fig. 3A). The CV staining method is less reproducible than the alamarBlue assay and therefore was only significant at a concentration of 20

M.

Mature biofilms of the multiresistant strain of S. epidermidis BM492 were less susceptible to F 2,5,12 W treatment. However, a significant suppression of metabolic activity of the biofilm was found following treatment with 40 M F 2,5,12 W (Fig. 3D).

A c c e p t e d

M a n u s c r i p t 10

Discussion

The involvement of S. epidermidis in infections associated with indwelling medical devices as well as the rapid development of multiple antibiotic resistance have continued to receive significant attention [START_REF] Donlan | Biofilms: microbial life on surfaces[END_REF]. These types of infections are mediated by the propensity of S. epidermidis to adhere to surfaces and form biofilms. In general, bacteria grown as a biofilm are protected from host defences and often exhibit reduced susceptibility to antibiotics, contributing to the persistence of the biofilm [START_REF] Gilbert | Biofilm susceptibility to antimicrobials[END_REF]. F 2,5,12 W, a short peptide based on chicken CATH-2, has been reported to have broad antibacterial activity [START_REF] Molhoek | Chicken cathelicidin-2-derived peptides with enhanced immunomodulatory and antibacterial activities against biological warfare agents[END_REF], but so far its effects on bacterial biofilms had not been evaluated.

The primary aim of this study was to investigate the effects of F 2,5,12 W on S. epidermidis biofilms using a simple screening method to quantify biofilms in vitro. The viability of biofilm cells was considered as the most important aspect when evaluating the effect of antimicrobial agents. CV stains polysaccharides in the extracellular matrix but does not distinguish between living and dead cells. Therefore, a quantification model based on reduction of alamarBlue by metabolically active cells was also included.

First, the antibacterial activity of F 2,5,12 W against two S. epidermidis strains, including a multiresistant strain, was determined. The results demonstrated that F 2,5,12 W exhibits activity with a MIC value of 10 M for both strains. Therefore, the inhibitory effect of F 2,5,12 W on biofilm formation likely reflects decreased initial adhesion of bacteria to the surface and subsequent development of the biofilm.

Mature biofilms are known to be more refractory to eradication relative to inhibition of de novo biofilm formation. Therefore, the effect of peptide F 2,5,12 W on mature biofilms was evaluated. Low concentrations of F 2,5,12 W had a better antimicrobial effect on the weak biofilm-producing S. epidermidis strain BM185 than on the strong biofilm-forming strain BM492. Both the architecture of a thinner biofilm and the total number of bacteria in such a biofilm may explain this finding.

Several other groups have investigated the effect of natural HDPs or synthetic analogues on different types of biofilms [START_REF] Hell | Human cathelicidin peptide LL37 inhibits both attachment capability and biofilm formation of Staphylococcus epidermidis[END_REF][START_REF] Overhage | Human host defense peptide LL-37 prevents bacterial biofilm formation[END_REF][START_REF] Wei | Effect of MUC7 peptides on the growth of bacteria and on Streptococcus mutans biofilm[END_REF][START_REF] Singh | A component of innate immunity prevents bacterial biofilm development[END_REF][START_REF] Flemming | High in vitro antimicrobial activity of synthetic antimicrobial peptidomimetics against staphylococcal biofilms[END_REF][START_REF] Beckloff | Activity of an antimicrobial peptide mimetic against planktonic and biofilm cultures of oral pathogens[END_REF]. It was reported that various peptides could be used to prevent biofilm formation. However, only three groups reported that peptides were effective in killing bacteria in a 24-h-old biofilm [START_REF] Wei | Effect of MUC7 peptides on the growth of bacteria and on Streptococcus mutans biofilm[END_REF][START_REF] Flemming | High in vitro antimicrobial activity of synthetic antimicrobial peptidomimetics against staphylococcal biofilms[END_REF][START_REF] Beckloff | Activity of an antimicrobial peptide mimetic against planktonic and biofilm cultures of oral pathogens[END_REF].

In summary, we have demonstrated that at low concentrations, below those that kill or inhibit growth, F 2,5,12 W strongly prevents bacterial biofilm formation in vitro. Moreover,

  peptide C1-15, comprising the first 15 N-terminal amino acids of CATH-2,

  broth (MHB) containing 0.5% glucose at a volume of 50 L/well in 96-well flat-bottom tissue culture plates (Costar ® ; Corning Inc., Corning, NY). Overnight cultures of S. epidermidis were diluted 1:50 in MHB containing 0.5% glucose and 50 L of this diluted culture was added to each well and plates were incubated for 24 h at 35 C. The final concentration of F 2,5,12 W ranged between 0.6 M and 40 M. The minimal inhibitory concentration (MIC) was defined as the lowest peptide concentration that yielded no visible growth.

A

  c c e p t e d M a n u s c r i p t 11 In addition, it was demonstrated that F 2,5,12 W suppressed biofilm formation of both strains at a concentration four times below the MIC. Inhibition by F 2,5,12 W of biofilm formation may be the result of the inhibition of bacterial growth. However, inhibition of biofilm formation was observed at concentrations that did not affect bacterial growth.

F 2 , 5 ,

 25 [START_REF] Wei | Effect of MUC7 peptides on the growth of bacteria and on Streptococcus mutans biofilm[END_REF] W also impairs existing S. epidermidis biofilms.

15 Fig. 1 .

 151 Fig. 1. Biomass of a 24-h-old biofilm of two different Staphylococcus epidermidis strains,

Fig. 2 .

 2 Fig. 2. Effect of F 2,5,12 W on Staphylococcus epidermidis biofilm formation by (A,C) strain

Fig. 3 .

 3 Fig. 3. Effect of F 2,5,12 W on Staphylococcus epidermidis pre-formed biofilm of (A,C)

  ) (A,B) and metabolic activity was determined by the alamarBlue ® reduction assay (grey bars) (C,D). OD 630 , optical density at 630 nm. A c c e p t e d M a n u s c r i p

Funding This work was supported financially by the Dutch Ministry of Defense (V502 program) and Immuno Valley (ALTANT program of the Dutch Ministry of Agriculture, Nature and Food Quality).

Competing interests

None declared.

Ethical approval

Not required.

A c c e p t e d M a n u s c r i p

Edited Figure 1 A c c e p t e d M a n u s c r i p Edited Figure 3


Nicolas Brisebarre

Mioara Joldeş

Érik Martin-Dorel

Micaela Mayero

Jean-Michel Muller

Ioana Paşca

Laurence Rideau

Laurent Théry

Rigorous Polynomial Approximation using Taylor Models in Coq

Keywords: certified error bounds, Taylor models, Coq proof assistant, rigorous polynomial approximation

One of the most common and practical ways of representing a real function on machines is by using a polynomial approximation. It is then important to properly handle the error introduced by such an approximation. The purpose of this work is to offer guaranteed error bounds for a specific kind of rigorous polynomial approximation called Taylor model. We carry out this work in the Coq proof assistant, with a special focus on genericity and efficiency for our implementation. We give an abstract interface for rigorous polynomial approximations, parameterized by the type of coefficients and the implementation of polynomials, and we instantiate this interface to the case of Taylor models with interval coefficients, while providing all the machinery for computing them. We compare the performances of our implementation in Coq with those of the Sollya tool, which contains an implementation of Taylor models written in C. This is a milestone in our long-term goal of providing fully formally proved and efficient Taylor models.

Rigorous Approximation of Functions by Polynomials

It is frequently useful to replace a given function of a real variable by a simpler function, such as a polynomial, chosen to have values close to those of the given function, since such an approximation may be more compact to represent and store but also more efficient to evaluate and manipulate. As long as evaluation is concerned, polynomial approximations are especially important. In general the basic functions that are implemented in hardware on a processor are limited to addition, subtraction, multiplication, and sometimes division. Moreover, division is significantly slower than multiplication. The only functions of one variable that one may evaluate using a bounded number of additions/subtractions, multiplications and comparisons are piecewise polynomials: hence, on such systems, polynomial approximations are not only a good choice for implementing more complex functions, they are frequently the only one that makes sense.

Polynomial approximations for widely used functions used to be tabulated in handbooks [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF]. Nowadays, most computer algebra systems provide routines for obtaining polynomial approximations of commonly used functions. However, when bounds for the approximation errors are available, they are not guaranteed to be accurate and are sometimes unreliable.

Our goal is to provide efficient and quickly computable rigorous polynomial approximations, i.e., polynomial approximations for which (i) the provided error bound is tight and not underestimated, (ii) the framework is suitable for formal proof (indeed, the computations are done in a formal proof checker), while requiring computation times similar to those of a conventional C implementation.

Motivations

Most numerical systems depend on standard functions like exp, sin, etc., which are implemented in libraries called libms. These libms must offer guarantees regarding the provided accuracy: they are heavily tested before being published, but for precisions higher that single precision, an exhaustive test is impossible [START_REF] De Dinechin | Assisted verification of elementary functions using Gappa[END_REF]. Hence a proof of the behavior of the program that implements a standard function should come with it, whenever possible. One of the key elements of such a proof would be the guarantee that the used polynomial approximation is within some threshold from the function. This requirement is even more important when correct rounding is at stake. Most libms do not provide correctly rounded functions, although the IEEE 754-2008 Standard for Floating-Point (FP) Arithmetic [START_REF][END_REF] recommends it for a set of basic functions. Implementing a correctly rounded function requires rigorous polynomial approximations at two steps: when actually implementing the function in a given precision, and-before that-when trying to solve the table maker's dilemma for that precision.

The 1985 version of the IEEE Standard for FP Arithmetic requires that the basic operations (+, -, ×, ÷, and √ •) should produce correctly rounded results, as if the operations were first carried out in infinite precision and these intermediate results were then rounded. This contributed to a certain level of portability and provability of FP algorithms. Until 2008, there was no such analogous requirement for standard functions. The main impediment for this was the table maker's dilemma, which can be stated as follows: consider a function f and a FP number x. In most cases, y = f (x) cannot be represented exactly. The correctly rounded result is the FP number that is closest to y. Using a finite precision environment, only an approximation ŷ to y can be computed. If that approximation is not accurate enough, one cannot decide the correct rounding of y from ŷ. Ziv [START_REF] Ziv | Fast evaluation of elementary mathematical functions with correctly rounded last bit[END_REF] suggested to improve the accuracy of the approximation until the correctly rounded value can be decided. A first improvement over that approach derives from the availability of tight bounds on the worst-case accuracy required to compute some functions [START_REF] Lefèvre | Worst cases for correct rounding of the elementary functions in double precision[END_REF], which made it possible to write a libm with correctly rounded functions, where correct rounding is obtained at modest additional costs [START_REF] Project | CRlibm, Correctly Rounded mathematical library[END_REF]. The TaMaDi project [START_REF] Muller | Projet ANR TaMaDi -Dilemme du fabricant de tables -Table Maker's Dilemma[END_REF] aims at computing the worst-case accuracy for the most common functions and formats. Doing this requires very accurate polynomial approximations that are formally verified.

Beside the Table Maker's Dilemma, the implementation of correctly rounded elementary functions is a complex process, which includes finding polynomial approximations for the considered function that are accurate enough. Obtaining good polynomial approximations is detailed in [START_REF] Brisebarre | Computing Machine-efficient Polynomial Approximations[END_REF][START_REF] Brisebarre | Efficient polynomial L ∞ -approximations[END_REF][START_REF] Chevillard | Évaluation efficace de fonctions numériques[END_REF]. In the same time, the approximation error between the function and the polynomial is very important since one must make sure that the approximation is good enough. The description of a fast, automatic and verifiable process was given in [START_REF] Joldeş | Rigourous Polynomial Approximations and Applications[END_REF].

In the context of implementing a standard function, we are interested in finding polynomial approximations for which, given a degree n, the maximum error between the function and the polynomial is minimum: this "minimax approximation" has been broadly developed in the literature and its application to function implementation is discussed in detail in [START_REF] Chevillard | Évaluation efficace de fonctions numériques[END_REF][START_REF] Muller | Elementary Functions, Algorithms and Implementation[END_REF]. Usually this approximation is computed numerically [START_REF] Remez | Sur un procédé convergent d'approximations successives pour déterminer les polynômes d'approximation[END_REF], so an a posteriori error bound is needed. Obtaining a tight bound for the approximation error reduces to computing a tight bound for the supremum norm of the error function over the considered interval. Absolute error as well as relative errors can be considered. For the sake of simplicity, in this paper, we consider absolute errors only (relative errors would be handled similarly). Our problem can be seen as a univariate rigorous global optimization problem, however, obtaining a tight and formally verified interval bound for the supremum norm of the error function presents issues unsuspected at a first sight [START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF], so that techniques like interval arithmetic and Taylor models are needed. An introduction to these concepts is given below.

Interval arithmetic and Taylor models. The usual arithmetic operations and functions are straightforwardly extended to handle intervals. One use of interval arithmetic is bounding the image of a function over an interval. Interval calculations frequently overestimate the image of a function. This phenomenon is in general proportional to the width of the input interval. We are therefore interested in using thin input intervals in order to get a tight bound on the image of the function. While subdivision methods are successfully used in general, when trying to solve this problem, one is faced with what is known as a "dependency phenomenon": since function f and its approximating polynomial p are highly correlated, branch and bound methods based on using intervals of smaller width to obtain less overestimation, end up with an unreasonably high number of small intervals. To reduce the dependency, Taylor models are used. They are a basic tool for replacing functions with a polynomial and an interval remainder bound, on which basic arithmetic operations or bounding methods are easier.

Related work

Taylor models [START_REF] Makino | Rigorous Analysis of Nonlinear Motion in Particle Accelerators[END_REF][START_REF] Neumaier | Taylor forms -use and limits[END_REF][START_REF] Makino | Taylor models and other validated functional inclusion methods[END_REF] are used for solving rigorous global optimization problems [START_REF] Makino | Rigorous Analysis of Nonlinear Motion in Particle Accelerators[END_REF][START_REF] Berz | Rigorous global search using Taylor models[END_REF][START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF][START_REF] Berz | Long-term stability of the tevatron by verified global optimization[END_REF] and obtaining validated solutions of ODEs [START_REF] Neher | On Taylor model based integration of ODEs[END_REF] with applica-tions to critical systems like particle accelerators [START_REF] Berz | Long-term stability of the tevatron by verified global optimization[END_REF] or robust space mission design [START_REF] Lizia | Robust Space Trajectory and Space System Design using Differential Algebra[END_REF]. Freely available implementations are scarce. One such implementation is available in Sollya [START_REF] Chevillard | Sollya: An environment for the development of numerical codes[END_REF]. It handles univariate functions only, but provides multiple-precision support for the coefficients. It was used for proving the correctness of supremum norms of approximation errors in [START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF], and so far it is the only freely available tool that provides such routines. However, this remains a C implementation that does not provide formally proved Taylor models, although this would be necessary for having a completely formally verified algorithm.

There have been several attempts to formalize Taylor models (TMs) in proof assistants. An implementation of multivariate TMs is presented in [START_REF] Zumkeller | Formal Global Optimization with Taylor Models[END_REF]. They are implemented on top of a library of exact real arithmetic, which is more costly than FP arithmetic. The purpose of that work is different than ours. It is appropriate for multivariate polynomials with small degrees, while we want univariate polynomials and high degrees. There are no formal proofs for that implementation. An implementation of univariate TMs in PVS is presented in [START_REF] Cháves | Utilisation et certification de l'arithmétique d'intervalles dans un assistant de preuves. These[END_REF]. Though formally proved, it contains ad-hoc models for a few functions only, and it is not efficient enough for our needs, as it is unable to produce Taylor models of degree higher than 6. Another formalization of Taylor models in Coq is presented in [START_REF] Collins | A Taylor Function Calculus for Hybrid System Analysis: Validation in Coq[END_REF]. It uses polynomials with FP coefficients. However, the coefficients are axiomatized, so we cannot compute the actual Taylor model in that implementation. We can only talk about the properties of the involved algorithms.

Our purpose is to provide a modular implementation of univariate Taylor models in Coq, which is efficient enough to produce very accurate approximations of elementary real functions. We start by presenting in Section 2 the mathematical definitions of Taylor models as well as efficient algorithms used in their implementation. We then present in Section 3 the Coq implementation. Finally we evaluate in Section 4 the quality of our implementation, both from the point of view of efficient computation and of numerical accuracy of the results.

Presentation of the Taylor Models

Definition, Arithmetic

A Taylor model (TM) of order n for a function f which is supposed to be n + 1 times differentiable over an interval [a, b], is a pair (T, ∆) formed by a polynomial T of degree n, and an interval part ∆, such that f

(x) -T (x) ∈ ∆, ∀x ∈ [a, b].
The polynomial can be seen as a Taylor expansion of the function at a given point. The interval remainder ∆ provides an enclosure of the approximation errors encountered (truncation, roundings).

For usual functions, the polynomial coefficients and the error bounds are computed using the Taylor-Lagrange formula and recurrence relations satisfied by successive derivatives of the functions. When using the same approach for composite functions, the error we get for the remainder is too pessimistic [START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF]. Hence, an arithmetic for TMs was introduced: simple algebraic rules like addition, multiplication and composition with TMs are applied recursively on the structure of function f , so that the final model obtained is a TM for f over [a, b]. Usually, the use of these operations with TMs offers a much tighter error bound than the one directly computed for the whole function [START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF]. For example, addition is defined as follows: let two TMs of order n for f 1 and f 2 , over [a, b]: (P 1 , ∆ 1 ) and (P 2 , ∆ 2 ). Their sum is an order n TM for f 1 + f 2 over [a, b] and is obtained by adding the two polynomials and the remainder bounds:

(P 1 , ∆ 1 ) + (P 2 , ∆ 2 ) = (P 1 + P 2 , ∆ 1 + ∆ 2 )
. For multiplication and composition, similar rules are defined.

We follow the definitions in [START_REF] Joldeş | Rigourous Polynomial Approximations and Applications[END_REF][START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF], and represent the polynomial T with tight interval coefficients. This choice is motivated by the ease of programming (rounding errors are directly handled by the interval arithmetic) and also by the fact that we want to ensure that the true coefficients of the Taylor polynomial lie inside the corresponding intervals. This is essential for applications that need to handle removable discontinuities [START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF]. For our formalization purpose, we recall and explain briefly in what follows the definition of valid Taylor models [23, Def. 2.1.3], and refer to [START_REF] Joldeş | Rigourous Polynomial Approximations and Applications[END_REF]Chap. 2] for detailed algorithms regarding operations with Taylor models for univariate functions.

Valid Taylor Models

A Taylor model for a function f is a pair (T, ∆). The relation between f and (T, ∆) can be rigorously formalized as follows.

Definition 1. Let f : I → R be a function, x 0 be a small interval around an expansion point x 0 . Let T be a polynomial with interval coefficients a 0 , . . . , a n and ∆ an interval. We say that (T, ∆) is a Taylor model of f at x 0 on I when

   x 0 ⊆ I and 0 ∈ ∆, ∀ξ 0 ∈ x 0 , ∃α 0 ∈ a 0 , . . . , α n ∈ a n , ∀x ∈ I, ∃δ ∈ ∆, f (x) - n i=0 α i (x -ξ 0 ) i = δ.
Informally, this definition says that there is always a way to pick some values α i in the intervals a i so that the difference between the resulting polynomial and f around x 0 is contained in ∆. This validity is the invariant that is preserved when performing operations on Taylor models. Obviously, once a Taylor model (T, ∆) is computed, if needed, one can get rid of the interval coefficients a i in T by picking arbitrary α i and accumulating in ∆ the resulting errors.

Computing the Coefficients and the Remainder

We are now interested in an automatic way of providing the terms a 0 , . . . , a n and ∆ of Definition 1 for basic functions. It is classical to use the following Lemma 1 (Taylor-Lagrange Formula). If f is n + 1 times differentiable on a domain I, then we can expand f in its Taylor series around any point x 0 ∈ I and we have: ∀x ∈ I, ∃ξ between x 0 and x such that

f (x) = n i=0 f (i) (x 0 ) i! (x -x 0 ) i T (x) + f (n+1) (ξ) (n + 1)! (x -x 0 ) n+1 ∆(x,ξ)
.

Computing interval enclosures a 0 , . . . , a n , for the coefficients of T , reduces to finding enclosures of the first n derivatives of f at x 0 in an efficient way. The same applies for computing ∆ based on an interval enclosure of the n+1 derivative of f over I. However, the expressions for successive derivatives of practical functions typically become very involved with increasing n. Fortunately, it is not necessary to generate these expressions for obtaining values of {f (i) (x 0 ), i = 0, . . . , n}. For basic functions, formulas are available since Moore [START_REF] Moore | Methods and Applications of Interval Analysis[END_REF] (see also [START_REF] Griewank | Evaluating Derivatives -Principles and Techniques of Algorithmic Differentiation[END_REF]). There one finds either recurrence relations between successive derivatives of f , or a simple closed formula for them. And yet, this is a case-by-case approach, and we would like to use a more generic process, which would allow us to deal with a broader class of functions in a more uniform way suitable to formalization.

Recurrence Relations for D-finite Functions. An algorithmic approach exists for finding recurrence relations between the Taylor coefficients for a class of functions that are solutions of linear ordinary differential equations (LODE) with polynomial coefficients, called D-finite functions. The Taylor coefficients of these functions satisfy a linear recurrence with polynomial coefficients [START_REF] Stanley | Differentiably finite power series[END_REF]. Most common functions are D-finite, while a simple counter-example is tan. For any D-finite function one can generate the recurrence relation directly from the differential equation that defines the function, see, e.g., the Gfun module in Maple [START_REF] Salvy | Gfun: a Maple package for the manipulation of generating and holonomic functions in one variable[END_REF]. From the recurrence relation, the computation of the first n coefficients is done in linear time. Let us take a simple example and consider f = exp. It satisfies the LODE f = f , f (0) = 1, which gives the following recurrence for the Taylor coefficients (c n ) n∈N : (n + 1)c n+1 -c n = 0, c 0 = 1, whose solution is c n = 1/n!.

This property lets us include in the class of basic functions all the D-finite functions. We will see in Section 3.2 that this allows us to provide a uniform and efficient approach for computing Taylor coefficients, suitable for formalization. We note that our data structure for that is recurrence relation + initial conditions and that the formalization of the isomorphic transformation from the LODE + initial conditions, used as input in Gfun is subject of future research.

Formalization of Taylor Models in Coq

We provide an implementation 5 of TMs that is efficient enough to produce very accurate approximating polynomials in a reasonable amount of time. The work is carried out in the Coq proof assistant, which provides a formal setting where we will be able to formally verify our implementation. We wish to be as generic as possible. A TM is just an instance of a more general object called rigorous polynomial approximation (RPA). For a function f , a RPA is a pair (T, ∆) where T is a polynomial and ∆ an interval containing the approximation error between f and T . We can choose Taylor polynomials for T and get TMs but other types of approximation are also available like Chebyshev models, based on Chebyshev polynomials. This generic RPA structure will look like:

Structure rpa := { approx: polynomial; error: interval } In this structure, we want genericity not only for polynomial with respect to the type of its coefficients and to its physical implementation but also for the type for intervals. Users can then experiment with different combinations of datatypes. Also, this genericity lets us factorize our implementation and will hopefully facilitate the proofs of correctness. We implement Taylor models as an instance of a generic RPA following what is presented in Section 2. Before describing our modular implementation, we present the Coq proof assistant, the libraries we have been using and how computation is handled.

The Coq proof assistant

Coq [START_REF] Bertot | Interactive Theorem Proving and Program Development. Coq'Art: The Calculus of Inductive Constructions[END_REF] is an interactive theorem prover that combines a higher-order logic and a richly-typed functional programming language. Thus, it provides an expressive language for defining not only mathematical objects but also datatypes and algorithms and for stating and proving their properties. The user builds proofs in Coq in an interactive manner. In our development, we use the SSReflect [START_REF] Gonthier | A Small Scale Reflection Extension for the Coq system[END_REF] extension that provides its own tactic language and libraries.

There are two main formalizations of real numbers in Coq: an axiomatic one [START_REF] Mayero | Formalisation et automatisation de preuves en analyses réelle et numérique[END_REF] and a constructive one [START_REF] Geuvers | Constructive Reals in Coq: Axioms and Categoricity[END_REF]. For effective computations, several implementations of computable real numbers exist. A library for multiple-precision FP arithmetic is described in [START_REF] Boldo | Flocq: A Unified Library for Proving Floating-point Algorithms in Coq[END_REF]. Based on this library, an interval arithmetic library is defined in [START_REF] Melquiond | Proving Bounds on Real-Valued Functions with Computations[END_REF]. It implements intervals with FP bounds. Also, the libraries [START_REF] O'connor | Certified exact transcendental real number computation in coq[END_REF] and [START_REF] Krebbers | Computer Certified Efficient Exact Reals in Coq[END_REF] provide an arbitrary precision real arithmetic. All these libraries are proved correct by deriving a formal link between the computational reals and one of the formalizations of real numbers. We follow the same idea: implement a computable TM for a given function and formally prove its correctness with respect to the abstract formalization of that function in Coq. This is done by using Definition 1 and the functions defined in the axiomatic formalization.

The logic of Coq is computational: it is possible to write programs in Coq that can be directly executed within the logic. This is why the result of a computation with a correct algorithm can always be trusted. Thanks to recent progress in the evaluation mechanism [START_REF] Boespflug | Full reduction at full throttle[END_REF], a program in Coq runs as fast as an equivalent version written directly and compiled in OCaml. There are some restrictions to the programs that can be executed in Coq: they must always terminate and be purely functional, i.e., no side-effects are allowed. This is the case for the above mentioned computable real libraries. Moreover, they are defined within Coq on top of the multiple-precision arithmetic library based on binary tree described in [START_REF] Grégoire | A purely functional library for modular arithmetic and its application to certifying large prime numbers[END_REF]. So only the machine modular arithmetic (32 or 64 bits depending on the machine) is used in the computations in Coq.

For our development of Taylor models we use polynomials with coefficients being some kind of computable reals. Following the description in Section 2, we use intervals with FP bounds given by [START_REF] Melquiond | Proving Bounds on Real-Valued Functions with Computations[END_REF] as coefficients. Since the interval and FP libraries are proved correct, so is the arithmetic on our coefficients. By choosing a functional implementation for polynomials (e.g., lists), we then obtain TMs that are directly executable within Coq. Now, we describe in detail this modular implementation.

A Modular Implementation of Taylor Models

Coq provides three mechanisms for modularization: type classes, structures, and modules. Modules are less generic than the other two (that are first-class citizens) but they have a better computational behavior: module applications are performed statically, so the code that is executed is often more compact. Since our generic implementation only requires simple parametricity, we have been using modules.

First, abstract interfaces called Module Types are defined. Then concrete "instances" of these abstract interfaces are created by providing an implementation for all the fields of the Module Type. The definition of Modules can be parameterized by other Modules. These parameterized modules are crucial to factorize code in our structures.

Abstract polynomials, coefficients and intervals. We describe abstract interfaces for polynomials and for their coefficients using Coq's Module Type. The interface for coefficients contains the common base of all the computable real numbers we may want to use. Usually coefficients of a polynomial are taken in a ring. We cannot do this here. For example, addition of two intervals is not associative. Therefore, the abstract interface for coefficients contains the required operations (addition, multiplication, etc.) only, where some basic properties (associativity, distributivity, etc.) are ruled out. The case of abstract polynomials is similar. They are also a Module Type but this time parameterized by the coefficients. The interface contains only the operations on polynomials (addition, evaluation, iterator, etc.) with the properties that are satisfied by all common instantiations of polynomials. For intervals, we directly use the abstract interface provided by the Coq.Interval library [START_REF] Melquiond | Proving Bounds on Real-Valued Functions with Computations[END_REF].

Rigorous polynomial approximations.

We are now able to give the definition of our rigorous polynomial approximation.

Module RigPolyApprox (C : BaseOps)(P : PolyOps C)(I : IntervalOps). Structure rpa : Type := RPA { approx : P.T; error : I.type }.

The module is parameterized by C (the coefficients), by P (the polynomials with coefficients in C), and by I (the intervals). Generic Taylor polynomials. Before implementing our Taylor models, we use the abstract coefficients and polynomials to implement generic Taylor polynomials. These polynomials are computed using an algorithm based on recurrence relations as described in Section 2.3. This algorithm can be implemented in a generic way. It takes as argument the relation between successive coefficients, the initial conditions and outputs the Taylor polynomial.

We detail the example of the exponential, which was also presented in Section 2.3. The Taylor coefficients (c n ) n∈N satisfy (n + 1)c n+1 -c n = 0. The corresponding Coq code is Definition exp_rec (n : nat) u := tdiv u (tnat n).

where tdiv is the division on our coefficients and tnat is an injection of integers to our type of coefficients. We then implement the generic Taylor polynomial for the exponential around a point x 0 with the following definition.

Definition T_exp n x0 := trec1 exp_rec (texp x0) n.

In this definition, trec1 is the function in the polynomial interface that is in charge of producing a polynomial of size n from a recurrence relation of order 1 (here, exp_rec) and an initial condition (here, texp x0, the value of the exponential at x0). The interface also contains trec2 and trecN for producing polynomials from recurrences of order 2 and order N with the appropriate number of initial conditions. Having specific functions for recurrences of order 1 and 2 makes it possible to have optimized implementations for these frequent recurrences. All the functions we currently dispose of in our library are in fact defined with trec1 and trec2. We provide generic Taylor polynomials for constant functions, identity, x → 1

Taylor models. We implement TMs on top of the RPA structure by using polynomials with coefficients that are intervals with FP bounds, according to Section 2. Yet we are still generic with respect to the effective implementation of polynomials. For the remainder, we also use intervals with FP bounds. This datatype is provided by the Coq.Interval library [START_REF] Melquiond | Proving Bounds on Real-Valued Functions with Computations[END_REF], whose design is also based on modules, in such a way that it is possible to plug all the machinery on the desired kind of Coq integers (i.e., Z or BigZ).

In a TM for a basic function (e.g., exp), polynomials are instances of the generic Taylor polynomials implemented with the help of the recurrence relations. The remainder is computed with the help of the Taylor-Lagrange formula in Lemma 1. For this computation, thanks to the parameterized module, we reuse the generic recurrence relations. The order-n Taylor model for the exponential on interval X expanded at the small interval X0 is as follows:

Definition TM_exp (n : nat) X X0 := RPA (T_exp n X0) (Trem T_exp n X X0).

We implement Taylor models for the addition, multiplication, and composition of two functions by arithmetic manipulations on the Taylor models of the two functions, as described in Section 2. Here is the example of addition: Definition TM_add (Mf Mg : rpa) := RPA (P.tadd (approx Mf) (approx Mg)) (I.add (error Mf) (error Mg)).

The polynomial approximation is just the sum of the two approximations and the interval error is the sum of the two errors. Multiplication is almost as intuitive. We consider the truncated multiplication of the two polynomials and we make sure that the error interval takes into account the remaining parts of the truncated multiplication. Composition is more complex. It uses addition and multiplication of Taylor polynomials. Division of Taylor models is implemented in term of multiplication and composition with the inverse function x → 1/x. The corresponding algorithms are fully described in [START_REF] Joldeş | Rigourous Polynomial Approximations and Applications[END_REF].

Discussion on the formal verification of Taylor models. The Taylor model Module also contains a version of Taylor polynomials defined with axiomatic real number coefficients. These polynomials are meant to be used only in the formal verification when linking the computable Taylor models to the corresponding functions on axiomatic real numbers. This link is given by Definition 1. The definition can be easily formalized in the form of a predicate validTM.

Definition validTM X X0 M f := I.subset X0 X /\ contains (error M) 0 /\ let N := tsize (approx M) in forall x0, contains X0 x0 -> exists P, tsize P = N /\ ( forall k, (k < N) -> contains (tnth (approx M) k) (tnth P k) ) /\ forall x, contains X x -> contains (error M) (f x -teval P (x -x0)).

The theorem of correctness for the Taylor model of the exponential TM_exp then establishes the link between the model and the exponential function Rexp that is defined in the real library.

Lemma TM_exp_correct : forall X X0 n, validTM X X0 (TM_exp n X X0) Rexp.

Our goal is to formally prove the correctness of our implementation of Taylor models. We want proofs that are generic, so a new instantiation of the polynomials would not require changing the proofs. In a previous version of our Coq development we had managed to prove correct Taylor models for some elementary functions and addition. No proofs are available yet for the version presented here but adapting the proofs to this new setting should be possible.

Benchmarks

We want to evaluate the performances of our Coq implementation of Taylor models. For this we compare them to those of Sollya [START_REF] Chevillard | Sollya: An environment for the development of numerical codes[END_REF], a tool specially designed to handle such numerical approximation problems.

The Coq Taylor models we use for our tests are implemented with polynomials represented as simple lists with a linear access to their coefficients. The coefficients of the approximating polynomial in our instantiation of Taylor models as well as the interval errors are implemented by intervals with multiple-precision FP bounds as available in the Coq.Interval library described in [START_REF] Melquiond | Proving Bounds on Real-Valued Functions with Computations[END_REF]. Since we need to evaluate the initial conditions for recurrences, only the basic functions already implemented in Coq.Interval can have their corresponding Taylor models.

In Sollya, polynomials have interval coefficients and are represented by a (coefficient) array of intervals with multiple-precision FP bounds. Sollya's autodiff() function computes interval enclosures of the successive derivatives of a function at a point or over an interval, relying on interval arithmetic computations and recurrence relations similar to the ones we use in our Coq development. Thus, we use it to compute the Taylor models we are interested in.

Timings, accuracy and comparisons

We compare the Coq and the Sollya implementations presented above on a selection of several benchmarks. Table 1 gives the timings as well as the tightness obtained for the remainders. These benchmarks have been computed on a 4-core computer, Intel(R) Xeon(R) CPU X5482 @ 3.20GHz.

Each cell of the first column of Table 1 contains a target function, the precision in bits used for the computations, the order of the TM, and the interval under consideration. When "split" is mentioned, the interval has been subdivided into a specified amount of intervals of equal length (1024 subintervals for instance in line 3) and a TM has been computed over each subinterval. Each TM is expanded at the middle of the interval. The symbols RD t (ln 4), resp. RU t (ln 2), denote ln(4) rounded toward -∞, resp. ln(2) rounded toward +∞, using precision t.

Columns 2 and 3 give the total duration of the computations (for instance, the total time for computing the 1024 TMs of the third line) in Coq and Sollya respectively. Columns 4 and 5 present an approximation error obtained using Coq and Sollya, while the last column gives, as a reference, the true approximation error, computed by ad-hoc means (symbolically for instance), of the function by its Taylor polynomial. Note that when "split" is mentioned, the error presented corresponds to the one computed over the last subinterval (for instance, [2 -1/256, 2] for the arctan example). For simplicity, the errors are given using three significant digits.

In terms of accuracy, the Coq and Sollya results are close. We have done other similar checks and obtained the same encouraging results (the error bounds returned by Coq and Sollya have the same orders of magnitude). This does not prove anything but is nevertheless very reassuring. Proving the correctness of an implementation that produces too large bounds would be meaningless.

Coq is 6 to 10 times slower than Sollya, which is reasonable. This factor gets larger when composition is used. One possible explanation is that composition implies lots of polynomial manipulations and the implementation of polynomials as simple lists in Coq maybe too naive. An interesting alternative could be to use persistent arrays [START_REF] Armand | Extending Coq with Imperative Features and Its Application to SAT Verification[END_REF] to have more efficient polynomials. Another possible improvement is at algorithmic level: while faster algorithms for polynomial multiplication exist [START_REF] Von Zur Gathen | Modern computer algebra[END_REF], currently in all TMs related works O(n 2 ) naive multiplication is used. We could improve that by using a Karatsuba-based approach, for instance.

Conclusion and Future Works

We have described an implementation of Taylor models in the Coq proof assistant. Two main issues have been addressed. The first one is genericity. We wanted our implementation to be applicable to a large class of problems. This motivates our use of modules in order to get this flexibility. The second issue is efficiency. Working in a formal setting has some impact in terms of efficiency. Before starting to prove anything, it was then crucial to evaluate if the compu-tational power provided by Coq was sufficient for our needs. The results given in Section 4 clearly indicate that what we have is worth proving formally. We are in the process of proving the correctness of our implementation. Our main goal is to prove the validity theorem given in Section 2 formally. This is tedious work but we believe it should be completed in a couple of months. As we aim at a complete formalization, a more subtle issue concerns the Taylor models for the basic functions and in particular how the model and its corresponding function can be formally related. This can be done in an ad-hoc way, deriving the recurrence relation from the formal definition. An interesting future work would be to investigate a more generic approach, trying to mimic what is provided by the Dynamic Dictionary of Mathematical Functions [START_REF] Benoit | The Dynamic Dictionary of Mathematical Functions (DDMF)[END_REF] in a formal setting.

Having Taylor models is an initial step in our overall goal of getting formally proved worst-case accuracy for common functions and formats. A natural next step is to couple our models with some positivity test for polynomials, for example some sums-of-squares technique. This would give us an automatic way of verifying polynomial approximations formally. It would also provide another way of evaluating the quality of our Taylor approximations. If they turned out to be not accurate enough for our needs, we could always switch to better kinds of approximations such as Chebyshev truncated series, thanks to our generic setting.

Table 1 .

 1 Benchmarks and timings for our implementation in Coq. .39s 6.23 × 10 -182 6.22 × 10 -182 6.22 × 10 -182 10 -29 2.93 × 10 -29 2.85 × 10 -29 .92 × 10 -50 6.10 × 10 -50 5.89 × 10 -50 .01 × 10 -43 9.33 × 10 -44 8.97 × 10 -44 .90 × 10 -47 4.92 × 10 -47 4.90 × 10 -47

		Execution time		Approximation error
		Coq Sollya	Coq	Sollya	Mathematical
	exp
	prec=120, deg=20 I=[1, RD53(ln 4)]	7.40s	0.01s	7.90 × 10 -35 7.90 × 10 -35 6.57 × 10 -35
	no split
	exp
	prec=120, deg=8 I=[1, RD53(ln 4)]	20.41s 3.77s	3.34 × 10 -39 3.34 × 10 -39 3.34 × 10 -39
	split in 1024
	exp
	prec=600, deg=40 I=[RU113(ln 2), 1] split in 256 38.10s 16arctan
	prec=120, I=[1, 2] split in 256 7.43 × exp × sin 11.45s 1.03s
	prec=200, deg=10 I=[1/2, 1] split in 2048 1m22s 12.05s 6exp/sin
	prec=200, deg=10 I=[1/2, 1] split in 2048 3m41s 13.29s 4exp • sin
	prec=200, deg=10 I=[1/2, 1]	3m24s 12.19s 4
	split in 2048

It is available at http://tamadi.gforge.inria.fr/CoqApprox/

This research was supported by the TaMaDi project of the French ANR (ref. ANR-2010-BLAN-0203-01).


Nabil Rachdi

Jean-Claude Fort

Thierry Klein

STOCHASTIC INVERSE PROBLEM WITH NOISY SIMULATOR -APPLICATION TO AERONAUTICAL MODELby

Keywords: Inverse Problem, Uncertainty Analysis, Robust characterization, M-estimation

), arising from some variables of the model. The phenomenon we study is the fuel mass needed to link two given countries with a commercial aircraft, where we only consider the Cruise phase . From a data base of fuel mass consumptions during the cruise phase, we aim at identifying the Specic Fuel Consumption (SF C) in a robust way, given the uncertainty of the cruise speed V and the lift-to-drag ratio F . In this paper, we present an estimation procedure based on Maximum-Likelihood estimation, taking into account this uncertainty.

Résumé. Le problème inverse est une pratique assez courante en ingénierie, où le but est de déterminer les causes d'un certain phénomène à partir d'observations de ce dernier. Le phénomène mis en jeu est représenté par un modèle numérique, dont certaines composantes peuvent comporter une part de variabilité (voir [4]). Le phénomène étudié est la masse de fuel nécessaire pour eectuer une liaison xée avec un avion commercial, en ne considérant que la phase de Croisière. Le but étant, à partir de données de masses de fuel consommées en croisière, d'identier de manière robuste la consommation spécique SF C de la motorisation en tenant compte de l'incertitude sur la vitesse de croisière V et sur la nesse F de l'avion. Dans cet article, nous proposons une procédure d'estimation basée sur une méthode de maximum de vraisemblance, prenant en compte cette incertitude.

Introduction

One of engineering activities is to model real phenomena. Once a model is built (physical principles, state equations, etc.), some parameters have to be identied and some variables of the model may present some intrinsic variability. Hence, the identication of parameters should implicitly take into account the uncertainty of variables of the model. In this paper, we present a likelihood-based method to estimate aeronautic parameters in a Fuel mass model. We use an analytical model that can be viewed as a black-box simulator. From a data base M * ,1 f uel , ..., M * ,n f uel giving the mass of fuel consumed for n lines between two given cities with a specic commercial aircraft, we aim at identifying the Specic Fuel Consumption (SF C) which corresponds to a characteristic value of engines. The model we use depends in particular on the cruise speed (V ) and on the lift-to-drag ratio (F ). These variables present intrinsic variability: the cruise speed may depend on atmospheric conditions and the lift-to-drag ratio is also subjected to variability potentially caused by turbulent phenomena. As a matter of fact, the identication of the parameter SF C should take into account the variability of the cruise speed and the lift-to-drag ratio.

In this paper, we propose an algorithm taken from the work of N. Rachdi et al. [START_REF] Rachdi | Risk bounds for new M-estimation problems[END_REF].

It allows a characterization of SF C from the observed data M * ,1 f uel , ..., M * ,n f uel and model simulations when the number of observations n is small. This article is organized as follows. In Section 2 we describe the setting of the problem. In Section 3 we build the algorithm for the inverse problem with a Maximum-Likelihood based method. In Section 4 we apply the algorithm given in Section 3. In Section 5 we illustrate the eect of modeling conditions, particularly the random modeling on SF C and the number of observed data. In Section 6 we establish the Theorem 6.2 providing an upper bound of the estimation error of the proposed algorithm. Section 7 is devoted to proving the Theorem 6.2.

2. General setting 2.1. Observations. In our study, the data M * ,1 f uel , ..., M * ,n f uel were generated from an aeronautic software which simulates gas turbine congurations used for power generation. In particular, it can simulate the consumed mass of fuel at some conguration of engines, altitude, speed, atmospheric conditions, etc. (See Figure 1). This software is very complex and very much time consuming. In fact only 200 outputs from this software are available by choosing various atmospheric conditions. This sample constitutes our data reference.

In a rst time, we pick up a small sample of size n = 32 from this reference sample. The data are given in Table 1.

Next, we will suppose that the observations M * ,1 f uel , ..., M * ,n f uel are drawn from an unknown probability distribution Q with associated Lebesgue density f with support The dierence M sup -M inf = 500 kg have to be thought as an overconsumption of about 7% approximatively.

I := [M inf = 7600, M sup = 8100] .

2.2.

A simplied aeronautical model. We recall that we are interested in identifying the specic fuel consumption SF C. It is a signicant factor determining the fuel eciency of a particular engine. To handle this problem we introduce a classical simplied Fuel mass model given by the Bréguet formula:

M f uel = (M empty + M pload ) e SF C•g•Ra V •F 10 -3 -1 . (1)
The xed variables are • M empty : Empty weight = basic weight of the aircraft (excluding fuel and passengers),

• M pload : Payload = maximal carrying capacity of the aircraft,

• g : Gravitational constant,

• Ra : Range = distance traveled by the aircraft.

The uncertain variables mentioned in the introduction are • V : Cruise speed = aircraft speed between ascent and descent phase, • F : Lift-to-drag ratio = aerodynamic coecient.

Noise modeling.

As said in the introduction, we have to take into account the uncertainty of the cruise speed V and the lift-to-drag ratio F . Given the nominal value of each variable (see Table 2), an expert judgment can derive the uncertainty bounds. (see Table 3).

The uncertainty on the cruise speed V represents a relative dierence of arrival time of 8 minutes. 4. The probability density function of a beta distribution on [a, b] with shape parameters (α, β) is

variable distribution parameter V Uniform (V min , V max ) F Beta (7, 2, F min , F max )
g (α,β,a,b) (x) = (x -a) (α-1) (b -x) β-1 (b -a) β-1 B(α, β) 1 [a,b](x) ,
where B(•, •) is the beta function.

Figure 2 shows the probability density functions of V and F . In order to emphasize the "noisy" feature of the variables V and F , we will use the writing

• V = V nom + V , • F = F nom + F , where V is a centered uniform random variable on the interval [ min V , max V ] with min V = V min -V nom and max V = V max -V nom .
Variable F , supposed to be independent of V , is a beta random variable on the interval [ min F , max F ] with shape parameters [START_REF] Rachdi | Risk bounds for new M-estimation problems[END_REF][START_REF] Billingsley | Convergence of probability measures[END_REF] where

min F = F min -F nom and max F = F max -F nom .
2.4. Robust identication of SF C. In our developments, we will not consider parameter SF C to be deterministic but it will be supposed random. We indeed do not only need to compute SF C. We also want to take into account its own variability in order to have a robust characterization of this parameter.

Assumption 2.1. Let us assume that the random variable SF C is compactly supported.

This assumption will allow to apply the theoretical framework (presented in annex) which needs the noise to be compactly supported. Such assumption does not impact the numerical results.

As a rst approach, let us assume that

SF C = µ SF C + σ SF C SF C , SF C ∼ N T (0, 1) . (2)
with unknown parameters µ SF C and σ SF C . We give the following ranges of variation Remark 2.2. All the numerical results we will present have been generated without truncation, which will not have a signicant eect. Indeed, the truncations are set to high quantiles (more than 99% for upper quantiles and less than 1% for lower quantiles).

Statistical modeling.

Let us denote by (E, P ) the probability space associated to the noise vector = ( SF C , V , F ) T , and denote the vector of parameters by

θ = (µ SF C , σ SF C ) T .
Then, we consider the analytical and simplied model of M f uel the mass of consumed fuel, as the function

M f uel = h( , θ) ,
where h : (E, P ) × Θ → I h is given by

h( , θ) = (M empty + M pload ) exp (µ SF C + σ SF C SF C ) • g • Ra (V nom + V ) • (F nom + F ) • 10 -3 -1 (3) with Θ = [15, 20] × [s, 1]
and I h is the interval

I h = h(E, Θ) = [M h inf , M h sup ] . (4)
We denote by |I h | its length.

Remark 2.3. We observe through simulations that

I ⊂ I h , where I = [M inf , M sup ] is the observation interval given above.
The purpose is now to estimate parameter θ ∈ Θ from the set of data M * ,1 f uel , ..., M * ,n f uel . In the next section, we propose an estimation procedure taken from [START_REF] Rachdi | Risk bounds for new M-estimation problems[END_REF].

Parameter estimation

In our previous framework it is possible to apply the procedures developed in [START_REF] Rachdi | Risk bounds for new M-estimation problems[END_REF].

In particular, we choose to work with the log -contrast which can be understood as a Maximum Likelihood based estimation.

The sample M * ,1 f uel , ..., M * ,n f uel is drawn from an unknown distribution Q. We will use the parametric family of distributions {Q θ , θ ∈ Θ} where Q θ is the pushforward measure of P by the (measurable) application u → h(u, θ) .

That means that we consider the models h( , θ) to be a reasonable rst approximation in order to obtain statistical information about Q.

Denoting ρ θ the Lebesgue density associated to the measure Q θ , the maximum likelihood procedure is given by

θ = Argmin θ∈Θ - 1 n n i=1 log(ρ θ (M * ,i f uel )) . (5)
However, the above procedure is unfeasible because the density ρ θ is not analytically tractable. As suggested in [START_REF] Rachdi | Risk bounds for new M-estimation problems[END_REF], we replace ρ θ by an estimator denoted ρ m θ . There are many ways to estimate a density. We choose the simplest one, which is commonly used in industrial modeling. Let 1 , ..., m be m random variables i.i.d from P , we consider the kernel estimate of our density

ρ m θ (•) = 1 m m j=1 K b m θ (• -h( j , θ)) , (6)
where K b m θ is the Gaussian kernel

K b m θ (x) = 1 √ 2 π b m θ e -x 2 2 (b m θ ) 2 ,
and b m θ is computed from the sample h( j , θ), j = 1, ..., m for θ ∈ Θ, by Silverman's rule-of-thumb :

b m θ = 4 3 1/5 m -1/5 σθ . (7)
The quantity σθ is the empirical standard deviation of the sample h( j , θ), j = 1, ..., m

σθ = 1 m m j=1 h( j , θ) - 1 m m j=1 h( j , θ) 2 .
Other popular estimates are truncated projections on a suitable basis of functions (sine, cosine, wavelets, etc.). Here we do not discuss the optimization of the density estimation, but all what follows could be applied in the same way. The numerical results may be slightly dierent, but qualitatively the same. Replacing ρ θ by ρ m θ in ( 5) and simplifying by the multiplying constant 1/n yields the estimation procedure

θ = Argmin θ∈Θ - n i=1 log 1 m m j=1 K b m θ h( j , θ) -M * ,i f uel . (8)
Hence, our problem is an inverse problem. More precisely, it is an inverse problem in presence of uncertainties, also called probabilistic inverse problem or stochastic inverse problem. This topic is often treated in the eld of uncertainty management: the goal may for instance be to identify the intrinsic uncertainty of a system, see for instance the PhD works [START_REF] Barbillon | méthodes d'interpolation à noyaux pour l'approximation de fonctions type boîte noire coûteuses[END_REF] and [START_REF] Kuhn | [END_REF]. Another reference is the paper of E. de Rocquigny and S. Cambier [START_REF] De Rocquigny | Inverse probabilistic modelling of the sources of uncertainty: a non-parametric simulated-likelihood method with application to an industrial turbine vibration assessment[END_REF], where the purpose is to identify a parameter of interest which controls the vibration amplication of stream turbines. Our framework is dierent. The main dierence lies in the absence of assumptions, in the present paper, on the distribution of the error between observation data and reference data. Thus, it dierentiates the estimation procedures we propose from the ones developed in [START_REF] De Rocquigny | Inverse probabilistic modelling of the sources of uncertainty: a non-parametric simulated-likelihood method with application to an industrial turbine vibration assessment[END_REF].

In the following section, we provide a numerical analysis using the algorithm given by [START_REF] Talagrand | Sharper bounds for Gaussian and empirical processes[END_REF]. Theoretical aspects will be addressed in Section 6.

4. Numerical study : rst approach 4.1. Estimation. Setting

J(θ) = - n i=1 log 1 m m j=1 K b m θ h( j , θ) -M * ,i f uel with θ = (µ SF C , σ SF C ) T ,
our problem is a minimization problem where we want to compute

θ = Argmin θ∈Θ J(θ) .
We recall that n = 32. The data (M * ,i f uel ) i=1,...,n are provided by Table 1. We choose m = 10000 (the number of calls to the model given by ( 3)), and for j = 1, ..., m, j ∼ P where

P (du, dv, dw) = 1 √ 2π L e -u 2 /2 g (7,2, min F , max F ) (v) 1 max V -min V 1 [ min V , max V ] (w) du dv dw ,
with L = Φ(3)-Φ(-3) (Φ is the cumulative distribution function of a standard Gaussian random variable). This optimization procedure can be solved by Quasi-Newton methods. We present the results in Table 5. In Figure 3 given by N T ( µ SF C , σ SF C ) . We notice that the minimum θ = ( µ SF C , σ SF C ) is correctly located.

Comparison with reference sample.

In order to analyse the results obtained in the previous subsection, we need some reference sample of SF C values at the same simulation conditions. We take for reference the sample of SF C values of size 200 described in the introduction, provided by the aeronautic software. The characteristics of this sample are given in Table 6 The data in Table 6 have to be compared with those in Table 5 where the mean and standard deviation are µ SF C = 17.397 and σ SF C = 0.201, respectively. Table 7 provides the associated relative errors. Figure 5 shows the histogram of the reference sample and the estimated distribution of SF C obtained in Figure 3. It appears that the location of the variable of interest SF C is well reached whereas the standard deviation estimation provides an error of 60%. The "error" has roughly two origins:

Reference sample Estimated

Statistical error : this error is mainly due to the limited number (n = 32) of data from the observed masses of fuel M * ,i f uel . It is also due to the error induced by the kernel approximation of ρ θ .Yet the choice of m = 10000 calls to the analytical model garantees that the error on ρ θ is small . Model error : this error is relative to the use of Fuel mass model (1) with uncertain variables V and F (Figure 2), and includes the Gaussian hypothesis for SF C (2).

Thus, the model error can be separated into 2 parts: physical model error and uncertainty modeling error. We observe on Figure 5 that the SF C parameter does not behave like a Gaussian variable. This can be qualied as model error.

However, if one just wants to estimate the mean value of SF C, the Gaussian hypothesis does not have a signicative impact (0.5% of error). On the other hand, if one wants more information about SF C, other modeling tools are needed to allow a robust In the next subsection, we will discuss the uncertainty modeling, more precisely, the Gaussian hypothesis for SF C given by (2).

5. On the probabilistic modeling of SF C 5.1. Considering Wiener-Hermite representation in the previous analysis.

The characterization of a random variable by the mean and the standard deviation only could be too approximative in order to study the whole behavior of the variable. In this study, we have made an a priori (a model) on the variable of interest SF C. In (2) we supposed that

SF C ∼ N T (µ SF C , σ 2 SF C
) , which we now rewrite

SF C = µ SF C + σ SF C ξ T , ξ T ∼ N T (0, 1) . ( 9
)
We will see that this (truncated) Gaussian hypothesis on SF C is a particular case of a more general representation.

The so called Wiener Chaos Expansion, developed in the 30's by Wiener [START_REF] Wiener | The homogeneous chaos[END_REF], gives a representation of any second-order random variable Z:

Z = ∞ l=0 z l Υ l ((ξ k ) k≥1
) , (with convergence in L 2 (P) ) [START_REF] Van Der Vaart | Weak Convergence and Empirical Processes[END_REF] where (ξ k ) k≥1 is a (innite) sequence of independent standard normal random variables and the Υ l 's are the multivariate Hermite polynomials. This expansion is also called Wiener-Hermite expansion.

In practice, we have to consider a nite sequence (ξ 1 , ..., ξ M ) where M is called the order of the expansion, and the sum in ( 10) is truncated at p which is the degree of the expansion. Hence, considering all M -dimensional Hermite polynomials of degree lower than p, the representation ( 10) is approximated by

Z Z p,M = P -1 l=0 z l Υ l (ξ) , ξ = (ξ 1 , ..., ξ M ) , (11)
where

P = (M + p)! M ! p! .
The integer P corresponds to the number of coecients to be estimated. For our purpose, by Assumption 2.1, we will consider the following approximation

Z Z p,M = P -1 l=0 z l Υ l (ξ T ) , ξ T = (ξ 1 T , ..., ξ M T ) , (12)
where ξ T is a vector of truncated standard Gaussian variables. Moreover, one can notice that by orthogonality arguments in (12), we have without truncation We illustrate this aspect concerning SFC in the next subsection.

E(Z p,M ) = z 0 = E(Z) (13)
5.2. Application to the Specic Fuel Consumption. In our purpose, if we suppose that E(SF C 2 ) < ∞ (it is implicitly supposed in the Gaussian hypothesis), we can set the following modeling

SF C p,M = P -1 l=0 z l Υ l (ξ T ) , ξ T = (ξ 1 T , ..., ξ M T ) , M, p ≥ 1
which we rewrite by (13)

SF C p,M = µ SF C + P -1 l=1 z l Υ l (ξ T ) , ξ T = (ξ 1 T , ..., ξ M T ) , M, p ≥ 1 . ( 15
)
It appears now that the Gaussian representation ( 9) is the particular case of (15) with p = 1 and M = 1. Moreover, in view of the Wiener representation [START_REF] Van Der Vaart | Weak Convergence and Empirical Processes[END_REF], the Gaussian one (9) may lead to a rough approximation (if SF C is not Gaussian) and thus contributes to a non negligible model error. It is clearly observed in Figure 5 where the reference data does not seem to be drawn from a Gaussian distribution. As a matter of fact, one can hope to reduce the model error (described in the previous subsection), at least the error corresponding to SF C modeling, by considering a less restrictive representation (15) with some appropriate order M ≥ 1 and degree p ≥ 1 . Let us consider the Wiener-Hermite expansion of order M = 2 and degree p = 2

SF C 2,2 = µ SF C + 5 l=1 θ l Υ l (ξ T ) , ξ T = (ξ 1 T , ξ 2 T ) or SF C 2,2 = µ SF C + θ 1 ξ 1 T + θ 2 ξ 2 T + θ 3 ξ 1 T ξ 2 T + θ 4 ((ξ 1 T ) 2 -1) + θ 5 ((ξ 2 T ) 2 -1) , (16)
that leads to estimate P = (2+2)! 2!2! = 6 coecients. Table 8 shows the result obtained by the algorithm developed in the previous section where we change the function h( , θ) in Let us compare the relative errors of the rst two statistical moments by considering SF C 1,1 (i.e the Gaussian hypothesis Table 7) and SF C 2,2 (Table 9). The Wiener-Hermite modeling seems to improve the mean estimation of SF C whereas the standard deviation is poorly estimated in the two cases with an error of about 60%. There is no signicative dierence between the two methods regarding the rst two moments. However, the behavior of density functions corresponding to SF C 1,1 (see Figure 5) and SF C 2,2 is clearly not the same. We present in Figure 6 the result obtained when SF C is modeled by a Wiener expansion of order M = 2 and degree p = 2.

(3) replacing σ SF C SF C by θ 1 ξ 1 T + θ 2 ξ 2 T + θ 3 ξ 1 T ξ 2 T + θ 4 ((ξ 1 T ) 2 -1) + θ 5 ((ξ 2 T ) 2 -1) , with = (ξ 1 T , ξ 2 T , V , F ) and θ = (µ SF C , θ 1 , ..., θ 5 ) . θ 0 θ 1 θ 2 θ 3 θ 4 θ 5 SF C 2,
The distribution of SF C given by the Wiener expansion in Figure 6 seems to have a behavior close to the reference sample one, despite the fact that there is a non negligible bias. As mentioned in the previous subsection, this is due to the statistical and model errors. Indeed, let us recall that we have at disposal n = 32 reference fuel masses from which we characterize the SF C parameter. It would be interesting to see what happens when adding reference fuel masses, i.e by reducing the statistical error.

Wiener-Hermite analysis with augmented reference fuel mass sample.

We present here numerical results obtained by adding 50 new samples from the data basis built with the complex software with the same initial conditions. Figure 7 shows the characterization of SF C obtained by a Wiener expansion of order M = 2 and degree p = 2 from the augmented reference fuel mass sample of size n = 82. The Table 10 gives the coecients corresponding to this simulation. Hence, by adding reference data we improve signicatively the characterization of SF C on the rst two statistical moments as well as on the whole probability density function of SF C. In the next section we introduce some knowledge on the SF C modeling through an expert judgment inducing a new and statistical modeling that will improve to be better. 5.4. Analysis with a "good" a priori knowledge. In the previous analyses, we only consider truncated Wiener-Hermite expansions. This is more of a mathematical hypothesis than a knowledge brought to the modeling. Suppose now that an expert judgment says that the distribution of the SF C is of exponential form. Mathematically, it is equivalent to supposing that the probability density of SF C belongs to the family

θ 0 θ 1 θ 2 θ 3 θ 4 θ 5 SF C
p(u; θ) = θ 2 e -θ 2 (u-θ 1 ) 1 [θ 1 ,+∞[ , θ = (θ 1 , θ 2 ) ∈ R + × R * + .

One can check that this suggestion induces the modeling

SF C exp = θ 1 - 1 θ 2 log(ξ) , ξ ∼ U([0, 1]) , ( 17
)
where U([0, 1]) is the uniform distribution on the interval [0, 1] . In order to satisfy Assumption 2.1, we consider the truncated version

SF C exp = θ 1 - 1 θ 2 log(ξ T ) , ξ T ∼ U([c, 1]) , (18)
for some small c > 0.

Remark 5.1. Representation (18) seems quite dierent from the one provided by the Wiener-Hermite expansions (see ( 9) and ( 16)). Yet, as the random variable SF C exp has nite variance, the modeling (18) could be seen as a practical alternative to a Wiener expansion [START_REF] Van Der Vaart | Weak Convergence and Empirical Processes[END_REF]. Such a Wiener expansion would be given by choosing a high order M exp and a high degree p exp in (12).

In what follows, we present the results of the numerical analysis corresponding to n = 32 and n = 82. We clearly see that the informative knowledge contributes to improving signicatively the characterization of the Specic Fuel Consumption. With n = 82 fuel mass data, the results are satisfying as shown in Figure 8 and Table 15. 5.5. Conclusion. Section 5 was dedicated to illustrating the eect of the modeling conditions for SF C characterization. In particular, we showed the impact of a "model error" through the modeling of the random variable SF C. We also illustrated how the statistical error, through the number of fuel mass data, appears in the performance of the estimation. In all cases, we computed a parameter θ. If we supposed that there is no model error, that is Q ∈ {Q θ , θ ∈ Θ} (where Q is the "true" distribution of M * f uel ), the error is only due to the limited number of data. So it makes sense to investigate the dierence θθ * , where θ * is dened as Q = Q θ * . If this is not the case it gives an insight on the statistical error part. It is the topic of the following section.

Theoretical result

In this paper, the study of the procedure performance (8) will be non-asymptotic, i.e for a xed number of observations M * ,i f uel (n) and a xed number of variables j (m). The asymptotic study is let to a forthcoming work.

The quality of such estimation procedure can be investigated by giving an upper bound of the distance between the reachable parameter θ and the best parameter θ * (unknown).

The latter can be seen as the parameter obtained if one has an innite number of observations M * ,i f uel and variables j . More precisely,

θ * = Argmin θ∈Θ E Q log ρ θ (M * f uel ) , (19)
where E Q log ρ θ (M * f uel ) can be seen as the "limit" of the quantity

1 n n i=1 log 1 m m j=1 K b m θ h( j , θ) -M * ,i f uel (20)
in [START_REF] Talagrand | Sharper bounds for Gaussian and empirical processes[END_REF] when n and m go to innity. The Maximum Likelihood equation ( 19) turns out to be the minimization of the Kullback-Leibler divergence between Q and the family {Q θ , θ ∈ Θ}, while equation ( 20) is a smoothed empirical counterpart.

We consider the model h( , θ) given in [START_REF] De Rocquigny | Inverse probabilistic modelling of the sources of uncertainty: a non-parametric simulated-likelihood method with application to an industrial turbine vibration assessment[END_REF], but what follows can be generalized to any other one. Then, denoting by • the Euclidian norm in R 2 , it makes sense to bound the quantity

θ -θ * 2 .
Let us denote by R(θ) := E Q log ρ θ (M * f uel ) , and by f the Lebesgue density associated to the measure Q . -A2 It exists η > 0 such that for all θ ∈ Θ, the density probability of h( , θ) we noted ρ θ , satises

ρ θ > η .
-A3 For all θ ∈ Θ, the second derivative of ρ θ , we note ρ θ , exists and

C := sup θ∈Θ ρ θ 2 < +∞ .
-A4 We suppose that

0 < δ < inf θ∈Θ σθ and sup θ∈Θ σθ < σ < +∞ ,
where σθ is dened in [START_REF] Rachdi | Risk bounds for new M-estimation problems[END_REF].

We prove the following consistency theorem: Theorem 6.2. Let us consider the estimator θ in (8) and the Assumptions (6.1).

Then, for all 0 < τ < 1/2, with probability at least

1 -2 τ θ -θ * 2 ≤ c 1 log(a 1 τ -1 ) n + c 2 log(a 2 τ -1 ) + c 3 m 1/10 √ m ,
for some constants c 1 , c 2 , c 3 , a 1 and a 2 .

The risk bound of this theorem seems surprising since we obtain a rate of n 1/4 , whereas one expects a rate close to √ n for the treated parametric problem. This theorem is a consistency result. Consequently, it does not give information about the rate of convergence. The obtained bound can be explained by the fact that, by Assumption A1, we have [START_REF] Rachdi | Risk bounds for new M-estimation problems[END_REF]) then obviously θθ * ≈ n -1/4 . However, the √ n-rate can be reached by considering the approach of Corollary 5.53 (pp. 77) in [START_REF] Van Der | Asymptotic statistics[END_REF] where an additional assumption is made on the risk function θ → R(θ). More precisely, this assumption relies on the function θ → log(ρ θ ) which is supposed to satisfy a Lipschitz condition. This work is let to a forthcoming paper which will deal with a central limit theorem for the parameter θ (the rate of convergence will therefore be reachable).

R( θ) -R(θ * ) ≈ θ -θ * 2 . Indeed, if R( θ) -R(θ * ) ≈ 1/ √ n (bound given in

Proof of Theorem 6.2

We give a general proof of Theorem 6.2.

By Assumption A1, we have the Taylor-Lagrange formula

R( θ) = R(θ * ) + 1 2 ( θ -θ * ) T ∇ 2 R(ξ) ( θ -θ * ) , (21)
for some ξ ∈ Θ . Then, we will use the following lemma Lemma 7.1 (Rayleigh's quotient). Let H be a real symmetric matrix p × p and denote by λ 1 < ... < λ p the ordered eigenvalues of H. It holds that for all x ∈ R p -{0}

λ 1 ≤ x T H x x T x ≤ λ p .
Now, applying this lemma with H = ∇ 2 R(ξ) and x = ( θθ * ) yields

λ min θ -θ * 2 ≤ ( θ -θ * ) T ∇ 2 R(ξ) ( θ -θ * ) ,
where λ min > 0 is the smallest eigenvalue of the set of matrices {∇ 2 R(θ), θ ∈ Θ} . Then, using this last inequality with equality (21) gives

θ -θ * 2 ≤ 2 λ min R( θ) -R(θ * ) . ( 22
)
The problem turns to bound the positive quantity R( θ) -R(θ * ), where θ is given by [START_REF] Talagrand | Sharper bounds for Gaussian and empirical processes[END_REF]. Such bound can be investigated by applying Theorem 4.1 in [START_REF] Rachdi | Risk bounds for new M-estimation problems[END_REF], which is a general result. We will aim at computing constants K τ 1 and K τ 2 such that, with high probability (23)

R( θ) -R(θ * ) ≤ 2 f 2 η 1 √ n η 2 δ 2 f 2 γ K τ 1 + 1 √ m 1 √ 2 π δ K τ 2 + 1 m 2/5 C (1.06 σ) 2 √ 3 .
In our framework, the main work is to compute the concentration constants K τ 1 and K τ 2 derived from [START_REF] Rachdi | Risk bounds for new M-estimation problems[END_REF] in the following particular case.

7.1. On concentration constants K τ 1 and K τ 2 . Since we apply Theorem 4.1 in [START_REF] Rachdi | Risk bounds for new M-estimation problems[END_REF], one has to guarantee that the constants K τ 1 and K τ 2 are nite. For this, let us rst recall some denitions and notations relative to empirical processes. Denition 7.2. Empirical process. Let W be some probability measure on some space T and let us suppose given a k i.i.d sample ξ 1 , ..., ξ k drawn from W . Let us denote by W k the empirical measure

W k := 1 k k i=1 δ ξ i
and G some class of real valued functions g : T → R . We call W -empirical process indexed by G the following application

G k : G -→ R g -→ G k := √ k T g(t) (W k -W ) (dt) ,
also written

G k g := 1 √ k k i=1 (g(ξ k ) -E W (g(ξ))) .
We denote the supremum of an empirical process by

G k G := sup g∈G |G k g| .
Following the proof lines of Theorem 2.1, Table 2 p.11 in [START_REF] Rachdi | Risk bounds for new M-estimation problems[END_REF] (giving classes of functions) and considering the inequality (23), it is easy to verify that K τ 1 is dened as for all n ≥ 1 ,

P( U n A ≤ K τ 1 ) ≥ 1 -τ (24) where U n is the Q-empirical process (Q n = 1 n n i=1 δ M * ,i f uel
) indexed by the class of functions

A = {y ∈ I -→ (y -λ) 2 , λ ∈ I h } (25)
where we recall

I = [M inf , M sup ] and I h = h(E, Θ) . Similarly, the constant K τ 2 is dened as follows for all m ≥ 1 , P( V m B ≤ K τ 2 ) ≥ 1 -τ (26)
where V m is the P -empirical process (P m = 1 m m j=1 δ j ) indexed by the class of functions

B = {x ∈ E -→ e -(h(x,θ)-λ) 2 /2 b 2 , (θ, λ, b) ∈ Θ × I h × [δ, σ]} . ( 27
)
By the writings (24) and (26), the constants K τ 1 and K τ 2 arise from the "concentration of the measure phenomenon" (see [START_REF] Ledoux | The concentration of measure phenomenon[END_REF], [START_REF] Billingsley | Convergence of probability measures[END_REF]). More precisely, these constants characterize the tightness of the sequences of random variables U n A (which is (M * ,i f uel ) i=1,...,n dependent) and V m B (which is ( j ) j=1,...,m dependent). Now, we aim at computing (upper bound) these constants using concentration inequalities where the classes of functions A and B will play a crucial role. In particular, we will apply the following theorem which is Theorem 2.14.9 in [START_REF] Van Der Vaart | Weak Convergence and Empirical Processes[END_REF].

Before, let us recall the denition of the bracketing numbers (taken from [START_REF] Van Der Vaart | Weak Convergence and Empirical Processes[END_REF] p. 83-85). Denition 7.3. Bracketing numbers. Let G be some class of functions on T and denote by W a probability measure on T . Given two functions l, u, the bracket [l, u] is the set of all functions g with l ≤ g ≤ u. An

ε-bracket is a bracket [l, u] with ||u-l|| 2,W < ε. The bracketing number N [ ] (ε, G, L 2 (W ))
is the minimum number of ε-brackets needed to cover the class of functions G.

The entropy with bracketing is the logarithm of the bracketing number.

Remark 7.4. The bracketing numbers measure the "size", the complexity of a class of functions. Theorem 7.5. Let G be a uniformly bounded class of (measurable) functions g : T → [0, 1] and denote by W a probability measure on T . If the class G satises, for some constants K and L

N [ ] (ε, G, L 2 (W )) ≤ K ε L for every 0 < ε < K . (28)
Then, for every t > 0,

P( G k G > t) ≤ D t √ L L e -2t 2 ,
for a constant D that only depends on K.

The proof of this theorem can be found in [START_REF] Talagrand | Sharper bounds for Gaussian and empirical processes[END_REF]. Now, let K τ be a constant (to determine) which satises

P( G k G ≤ K τ ) ≥ 1 -τ . This is equivalent to P( G k G > K τ ) ≤ τ . (29)
By Theorem 7.5, applied with t = K τ , we have

P( G k G > K τ ) ≤ D K τ √ L L e -2(K τ ) 2 . (30)
Hence, the constant K τ can be taken such that

D K τ √ L L e -2(K τ ) 2 ≤ τ ,
which is similar to

(K τ ) 2 - L 2 log(K τ ) ≥ log(a L,D τ -1 ) 2 , with a L,D = D √ L L . (31)
Then, for small enough τ > 0, let us consider the constant

K τ = log(a L,D τ -1 ) 2 (32)
which satises (31).

Finally, we see that the constant K τ can be characterized (only) by the class of functions G through the constants D and L provided by Theorem 7.5.

In our purpose, the classes of interest are A and B dened in (25) and (27), respectively. Next, one can easily check that these classes are uniformly bounded and it is suitable to work with normalized classes

Ā = α A + 1 β A A , (33)
B = α B + 1 β B B , (34)
such that all the functions take values in [0, 1] . Now, we have to prove that the classes Ā and B have polynomial bracketing numbers following (28). This will give the constants L Ā, D Ā and L B, D B needed to identify the key constants K τ 1 and K τ 2 dened in (24) and (26), respectively. where f is the density associated to the measure Q, and using the fact that I ⊂ I h , we obtain that

F 2,Q ≤ 4 β A M h sup f 2 . This last inequality yields N ε 2 F 2,Q , I h , | • | ≤ N β A ε 8 M h sup f 2 , I h , | • | .
Since I h = [M h inf , M h sup ], the quantity (covering number) in the right member is bounded by

8 |I h | M h sup f 2 β A ε , |I h | = M h sup -M h inf . We nally get N [ ] (ε, Ā, L 2 (Q)) ≤ 8 |I h | M h sup f 2 β A ε , that is N [ ] (ε, Ā, L 2 (Q)) ≤ K Ā ε L Ā ,
with L Ā = 1 and

K Ā = 8 |I h | M h sup f 2 β A
that determines D Ā by [START_REF] Talagrand | Sharper bounds for Gaussian and empirical processes[END_REF] .

A similar work gives the constant L B = 1 and a constant D B . But, the constants of interest K τ 1 and K τ 2 dened in (24) and ( 26) are relative to non normalized classes A and B. Let us remark that if Ḡ = α + 1 β G, then

G k Ḡ = 1 β G k G . (37)
Now, let us denote by K τ Ḡ the constant that satises P( G k Ḡ ≤ K τ Ḡ ) ≥ 1 -τ , and denote by K τ G the constant that satises P( G k G ≤ K τ G ) ≥ 1 -τ . By (37), it is easy to check that we can take

K τ G = β K τ Ḡ . We deduce that K τ 1 = β A K τ Ā and K τ 2 = β B K τ B .
Finally, by ( 22) and (23) we have with probability 1 -2τ

θ -θ * 2 ≤ 4 f 2 λ min η 1 √ n η 2 δ 2 f 2 γ K τ 1 + 1 √ m 1 √ 2 π δ K τ 2 + 1 m 2/5
C (1.06 σ) 2 √ 3 which we rewrite

θ -θ * 2 ≤ √ 2 c 1 √ n K τ Ā + √ 2 c 2 √ m K τ B +
c 3 m 1/5 with corresponding constants c 1 , c 2 and c 3 and K τ Ā, K τ B are given by ( 35) and (36).

That concludes the proof of Theorem 6.2.

Figure 1

 1 Figure 1. Aeronautic software Reference Fuel Masses [kg] 7918 7671 7719 7839 7912 7963 7693 7815 7872 7679 8013 7935 7794 8045 7671 7985 7755 7658 7684 7658 7690 7700 7876 7769 8058 7710 7746 7698 7666 7749 7764 7667

Figure 2 .

 2 Figure 2. (a) Uncertainty on F . (b) Uncertainty on V .

µ

  SF C ∈ [15, 20] and σ SF C ∈ [s, 1] , for a small s > 0 . The distribution N T (0, 1) of SF C is a symmetric truncated standard Gaussian on the interval [-3, 3]. Now, our problem amounts to estimating the location parameter µ SF C and the standard deviation σ SF C .

  we show the resulting probability density function of SF C estimator value of J( θ) estimated SF C location estimated SF C dispersion θ J( θ) = 199.465 µ SFC = 17.397 σ SFC = 0.201

Figure 3 .

 3 Figure 3. Estimated Specic Fuel Consumption distribution.

Figure 4

 4 Figure 4 provides prole views of the criterion function θ = (µ SF C , σ SF C ) → J(θ), rst at σ SF C = σ SF C (Figure 4(a), we show log(J(θ))) and then at µ SF C = µ SF C (Figure 4(b)).We notice that the minimum θ = ( µ SF C , σ SF C ) is correctly located.

Figure 4 .

 4 Figure 4. (a) Prole view of log(J) at σ SF C = σ SF C . (b) Prole view of J at µ SF C = µ SF C .

Figure 5 .

 5 Figure 5. Reference and estimated SF C distributions.

  that by the decompositionZ p,M = Z + Z p,M -Z ,each choice of p and M will induce a model error mod err := Z p,M -Z .

Figure 6 .

 6 Figure 6. Estimations of SF C probability density with a Wiener Expansion p = M = 2 .

Figure 7 .

 7 Figure 7. Characterization of SF C with an augmented sample of fuel mass

θ 1 θ 2 SFC

 2 exp 17.23 3.45

Figure 8 .

 8 Figure 8. Characterization of SFC with an exponential hypothesis

Assumptions 6 . 1 .

 61 Let us consider the following assumptions.-A1 The map θ → R(θ) is twice dierentiable with∇R(θ * ) = 0and has a symmetric positive denite Hessian matrix ∇ 2 R . Let us denote by λ min > 0 the smallest eigenvalue of the set of matrices {∇ 2 R(θ), θ ∈ Θ} .

7. 2 .

 2 Characterization of L Ā, D Ā, L B, D B. We consider the Theorem 2.7.11 in[START_REF] Van Der Vaart | Weak Convergence and Empirical Processes[END_REF] (p. 164) which deals with classes that are Lipschitz in a parameter. It reads: Theorem 7.6. Let G = {t ∈ T → g s (t) , s ∈ S} be a class of functions satisfying for all t ∈ T , s, s ∈ S , |g s (t) -g s (t)| ≤ d(s, s ) G(t) , for some metric d on S and some functionG : t → G(t).Then, for any normN [ ] (2ε G , G, • ) ≤ N (ε, S, d) ,where N (ε, S, d) is the minimal number of balls {r , d(r, s) < ε} of radius ε needed to cover the set S .In what follows, we detail the case of the class Ā. The case of the class B is exactly in the same spirit.Let us recall that Q is the probability measure considered on I (observation space) and that we haveĀ = {f λ : y ∈ I -→ α B + 1 β A (y -λ) 2 , λ ∈ I h } , where I = [M inf , M sup ] and I h = [M h inf , M h sup ] (with I ⊂ I h ). So |f λ (y) -f λ (y)| = 1 β A |(y -λ 1 ) 2 -(y -λ 2 ) 2 | ≤ |λ 1 -λ 2 | F (y) ,with F (y) = 2 β A (y + M h sup ), and by Theorem 7.6 applied with • = • 2,Q , it holds that N [ ] ( , Ā, L 2 (Q)) ≤ N 2 F 2,Q , I h , | • | . Moreover, since F 2,Q ≤ sup y∈I F (y) f 2

7. 3 .

 3 End of the proof. By the previous subsection, we get the constants K τ Ā and K τ B given by (32) with associated constants L and D:K τ Ā = log(a 1 τ -1 ) 2 , a 1 = a L Ā,D Ā = D Ā (35) K τ B = log(a 2 τ -1 ) 2 , a 2 = a L B ,D B = D B (36)where initially a L,D = D √ L L (by (31)).

Table 1 .

 1 Simulated mass of fuel consumptions from aeronautic software

Table 2

 2 gives the xed variables values and the nominal values considered for uncertain variables.

Table 2 .

 2 Values of Fuel mass model inputs

	input value or nominal value unit
	M empty	42600	kg
	M pload	19900	kg
	g	9.8	m/s 2
	Ra	3000	km
	V nom	231	m/s
	F nom	19
	variable nominal value min max V 231 226 234 F 19 18.7 19.05

Table 3 .

 3 Minimal and maximal values of uncertain variables

Moreover, specialists in turbine engineering propose to model the uncertainties as presented in Table

Table 4 .

 4 Uncertainty modeling

Uncertainty on the lift-to-drag ratio

							Uncertainty on F
	8
	6
	PDF
	4
	2
	0
	18.80	18.85	18.90	18.95	19.00	19.05	19.10	19.15
				x

Table 5 .

 5 SF C characterization parameters

Table 6 .

 6 . Reference sample characteristics

	Mean Stand. dev. Reference sample 17.49 0.57

Table 7 .

 7

	SF C (3) Relative error

Relative errors of the mean and deviation between reference SF C sample and the estimated model

[START_REF] De Rocquigny | Inverse probabilistic modelling of the sources of uncertainty: a non-parametric simulated-likelihood method with application to an industrial turbine vibration assessment[END_REF]

Table 8 .

 8 SF C characterization parameters

	Mean Stand. dev.	Reference sample from SF C 2,2 Relative error 17.49 17.470 0.11 % 0.57 0.230 59.65 %

2

17.470 0.047 0.054 0.182 0.103 0.063

Table 9 .

 9 Relative errors of the mean and standard deviation with SF C2,2

Table 10 .

 10 SF C characterization parameters from augmented fuel mass sample

	2,2 17.50 0.281 0.008 0.012 0.191 0.219
	Mean Stand. dev.	Reference sample from SF C 2,2 Relative error 17.49 17.50 0.06 % 0.57 0.404 29.12 %

Table 11 .

 11 Relative errors of the mean and standard deviation with SF C 2,2 from augmented fuel mass sample

Table 12 .

 12 Estimation of θ = (θ 1 , θ 2 ) when n = 32

	Mean Stand. dev.	Reference sample from SF C exp (n = 32) Relative error 17.49 17.52 0.17 % 0.57 0.29 49.12 %

Table 13 .

 13 Relative errors of the mean and standard deviation between reference SF C sample and SF C exp when n = 32.

	θ 1 SF C exp 16.95 2 θ 2

Table 14 .

 14 Estimation of θ = (θ 1 , θ 2 ) when n = 82

	Mean Stand. dev.	Reference sample from SF C exp (n = 82) Relative error 17.49 17.45 0.23 % 0.57 0.501 12.1 %

Table 15 .

 15 Relative errors of the mean and standard deviation between reference SF C sample and SF C exp when n = 82.