Mostafa Meliani
email: melimostafa@gmail.com

Nathalie Bartoli
email: nathalie.bartoli@onera.fr

Thierry Lefebvre
email: thierry.lefebvre@onera.fr

Mohamed-Amine Bouhlel
email: mbouhlel@umich.edu

Joaquim R R A Martins

Joseph Morlier
email: joseph.morlier@isae.fr

Multi-fidelity efficient global optimization: Methodology and application to airfoil shape design

Predictions and design engineering decisions can be made using a variety of information sources that range from experimental data to computer models. These information sources could consist of different mathematical formulations, different grid resolutions, different physics, or different modeling assumptions that simplify the problem. This leads to information sources with varying degrees of fidelity, each with an associated accuracy and querying cost. In this paper, we propose a novel and flexible way to use multi-fidelity information sources optimally in the context of airfoil shape optimization using both Xfoil and ADflow. The new developments are based on Bayesian optimization and kriging metamodeling and allow the aerodynamic optimization to be sped up. In a constrained optimization example with 15-design variables problem, the proposed approach reduces the total cost by a factor of two compared to a single Bayesian based fidelity optimization and by a factor of 1.5 compared to sequential quadratic programming.

Nomenclature

I. Introduction

T e use of multiple fidelities is particularly interesting for aircraft design decisions where the high-dimensionality of the design space and the cost of the high fidelity analyses make global optimization near impossible. Many efforts have used multiple information sources for surrogate modeling (different accuracies for the same quantity of interest) [START_REF] Winkler | Combining Probability Distributions from Dependent Information Sources[END_REF] or optimization [START_REF] Lam | Multifidelity optimization using statistical surrogate modeling for non-hierarchical information sources[END_REF][START_REF] Poloczek | Multi-Information Source Optimization[END_REF]. There have also been efforts in modeling the differences in fidelity in an additive way [START_REF] Lewis | A multigrid approach to the optimization of systems governed by differential equations[END_REF][START_REF] Huang | Sequential kriging optimization using multiple-fidelity evaluations[END_REF]. However these approaches fail to work in the general case where different assumptions about the physics of the problem exist between fidelity levels. In this context, the work of Kennedy and O'Hagan [START_REF] Kennedy | Bayesian calibration of computer models[END_REF] provided powerful tools for multi-fidelity surrogate modeling which Forrester et al. [START_REF] Forrester | Multi-fidelity optimization via surrogate modelling[END_REF] used for a multi-fidelity extension of the Efficient Global Optimization algorithm (EGO) [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF] for unconstrained problems. In 2013, Le Gratiet [START_REF] Gratiet | Multi-fidelity Gaussian process regression for computer experiments[END_REF] proposed an alternative formulation of multi-fidelity co-kriging with explicit contributions of each fidelity level to the overall model. In 2002, EGO was extended to constrained problems with the development of SEGO [START_REF] Sasena | Flexibility and efficiency enhancements for constrained global design optimization with Kriging approximations[END_REF]). More recently SEGO was coupled to MOE (mixture of experts) methodology in order to solve high dimensional variables and constrained optimization problems such as wing shape aerodynamic optimization in the SEGOMOE framework [START_REF] Bartoli | Improvement of efficient global optimization with mixture of experts: methodology developments and preliminary results in aircraft wing design[END_REF][START_REF] Bartoli | An adaptive optimization strategy based on mixture of experts for wing aerodynamic design optimization[END_REF][START_REF] Bartoli | Adaptive modeling strategy for constrained global optimization with application to aerodynamic wing design[END_REF]. SEGOMOE will be used in this work as a reference for surrogate based constrained optimizer with a single high fidelity (HF). Of course, we will also compare to the standard gradient-based optimizer SNOPT [START_REF] Gill | SNOPT: An SQP algorithm for large-scale constrained optimization[END_REF]. In practice, the SEGOMOE framework is able to treat either unconstrained and constrained problem with or without MOE for black box optimization problems.

In this paper these works are used to improve the efficiency of global optimization algorithms which suffer from the curse of dimensionality. We demonstrate that the use of multiple fidelities helps alleviate this curse and show that we can perform high-dimensional Bayesian optimization. We begin by reviewing kriging and multi-fidelity co-kriging formulations; we then present the multi-fidelity extension of EGO called MFEGO. We provide an analysis of the properties of the proposed algorithm including its global convergence. Finally, we present two kinds of results on both unconstrained and constrained optimizations of a subsonic airfoil with 9 and 15 design variables and compare the performance of MFEGO to standard methodologies such as EGO and SNOPT. The constrained version of the algorithms will be later referred to as SEGO and MFSEGO.

II. Kriging and multi-fidelity co-kriging description

There are multiple ways to build an approximation of a function, such as artificial neural networks (ANN), radial basis functions (RBF) or support vector regressions (SVR) [START_REF] Forrester | Multi-fidelity optimization via surrogate modelling[END_REF]. In all these metamodels (or surrogate models), a fundamental assumption is that the quantity of interest y(x) can be written y(x) = ŷ(x) + , where , the residuals, are independent identically distributed normal random variables, so that fitting the model ŷ(x) is performed by minimizing a measure over .

Following the same assumptions, the specificity of kriging is that the 'errors'-or more accurately, the deviations from the base model-in the predicted values ŷ, are not independent. Rather, we take the view that the errors are a systematic function of the locations of the samples. The kriging surrogate model ŷ(x) = m(x) + Z(x), is comprised of two parts: 'regression' term m(x) and a functional departure from that regression: Z(x) [START_REF] Sasena | Flexibility and efficiency enhancements for constrained global design optimization with Kriging approximations[END_REF]. We can write the regression term for the universal kriging as the following:

m(x) = k j=1 β j f j , (1)
where { f 1 , . . . , f k } are basis functions of m(x) and β j are coefficients weighting these basis functions. Ordinary kriging denotes the special case where k = 1 and f 1 = 1, meaning the regression function takes the form of a constant, leaving most of the prediction work to Z(x). Z(x) on the other hand is a Gaussian process with a covariance function:

cov(Z(w), Z(x)) = σ 2 R(w, x), (2)
where σ 2 is a scaling factor known as the process variance, w and x are two design points in R d and R(w, x) is the spatial correlation function of the process also known as the kernel of the process. The choice of the kernel influences greatly the way data is fitted. The kernel quantifies how quickly and smoothly the function changes from point x to point w.

One of the most commonly used kernels in kriging is the squared exponential correlation kernel:

R(z, x) = exp - d k=1 θ (k) (z (k) -x (k) ) 2 (3)
where x ∈ R d , z ∈ R d and θ ∈ R d is a vector of hyperparameters of the kriging model, denoting the correlation along the different axes of space. This Kernel has the particularity that it only depends on the weighted distance between points x and z and not their actual positions. It is called a stationary Kernel. Once we have defined the Kernel, and "fitted" the vector of hyperparameters θ using the vector X T = {x 1 , . . . , x n } (with x i ∈ R d ) yielding the responses Y T = {y 1 , . . . , y n } (with y i ∈ R), we can express the mean and covariance of the Gaussian Process Z(x):

µ Z (x) = r(x, X T ) R -1 (X T , X T )Y T (4) σ 2 Z (x) = σ 2 (1 -r(x, X T ) R -1 (X T , X T )r(x, X T )) (5)
where x ∈ R d is the prediction point, and X T is the locations of training set. R(X T , X T ) is the matrix of correlations among the training points. r(x, X T ), on the other hand, denotes the correlation between the prediction point and the training points. Note that Eq. ( 4) ensures the interpolation of the training points. In fact, if we make a prediction at training point x i , the vector r(x i , X T ) will correspond to the i th line of R(X T , X T ), so that (r(x i , X T ) R -1 (X T , X T )) will give the i th unit vector and so µ Z (x i ) = y i .

If we take into account the 'regression' term m(x) the prediction becomes:

µ(x) = m(x) + r(x) R -1 (Y T -m(x))
which we can also express using the basis introduced in Eq. (1):

µ(x) = f (x) β + r(x) R -1 (Y T -F β) (6)
σ 2 (x) = σ 2 1 -r(x) R -1 r(x) + f (x) -r(x) R -1 F F R -1 F f (x) -r(x) R -1 F ( 7
)
where β is the vector of coefficients β j introduced in Eq. ( 1). F is the matrix of the values of the regression basis function at the positions of the training points, whereas f (x) is the vector of values of these functions at the prediction point.

In the case of multi-fidelity surrogate models, we can make assumptions to simplify the problem and inform our model (which results in a decrease of the data needed to learn the model). Lewis and Nash [START_REF] Lewis | A multigrid approach to the optimization of systems governed by differential equations[END_REF] proposed an additive bridge to link the high and low fidelity of a multigrid approach:

f H F (x) = f LF (x) + γ(x). (8)
When using different grids for a partial differential equation, it is expected that the low and high fidelity would have the same scale and a very high correlation. The coarse and finer meshes are expected to capture the same low frequencies, the coarser mesh, however, would fail to capture high-frequency phenomena. γ(x) ∈ R can then be used to capture the difference. However, in general when modeling multiple fidelities that involve different assumptions on the physics of the problem (e.g., viscosity, compressibility, or turbulence in aerodynamics models), chances are that the low and high fidelities will have different scales and sometimes (in rare cases) poor correlations. Kennedy and O'Hagan [START_REF] Kennedy | Bayesian calibration of computer models[END_REF] proposed a new formulation that takes the correlation and scaling into account by introducing a factor ρ ∈ R in the formulation above:

f H F (x) = ρ f LF (x) + δ(x) with f LF (•) ⊥ δ(•) (9)
where δ(•) is the discrepancy function tasked with capturing the differences between the low-and high-fidelity functions beyond scaling. The addition of the term ρ increases the robustness of the model. Le Gratiet [START_REF] Gratiet | Multi-fidelity Gaussian process regression for computer experiments[END_REF] proposed an implementation using the regression term expression of universal kriging introduced in Eq. ( 1) and extends it to take the lower fidelity model as a basis function such that the regression term becomes:

m(x) = k j=1 β j f j + β ρ f LF ,
where β ρ is an estimation of ρ performed by a classic parameter estimation such as the likelihood maximization [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF][START_REF] Rasmussen | Gaussian processes for machine learning[END_REF].

Assuming the independence of the high and low-fidelity models (two levels here), the mean and variance of the high-fidelity model are expressed:

µ H F = ρ µ LF + µ δ ( 10
)
σ 2 H F = ρ 2 σ 2 LF + σ 2 δ ( 11
)
The approach can be extended to l levels of fidelity. To make this explicit, let us denote f 0 , . . . , f l the hierarchically ranked fidelity codes (from lowest f 0 = f LF to highest f l = f H F ). Using the recursive formulation, we write:

µ k = ρ k-1 µ k-1 + µ δ k (12) σ 2 k = ρ 2 k-1 σ 2 k-1 + σ 2 δ k . (13)
The formulation by Le Gratiet [START_REF] Gratiet | Multi-fidelity Gaussian process regression for computer experiments[END_REF], if satisfying the nested DOE requirement, offers explicit expressions of the contribution of fidelity levels to the uncertainty of the model. The nested DOE requirement states that X l ⊆ X l-1 . . . ⊆ X 0 where X i is the vector of training points of the fidelity i and l is the highest fidelity. With this assumption, a point computed at the highest fidelity has to be also computed at the lowest fidelities.

We introduce the following notation in Eq. ( 13):

σ 2 δ,k = σ 2 δ k for k ∈ {1, . . . , l} (14)
σ 2 δ,0 = σ 2 0 for k = 0.
We express the uncertainty contribution of the fidelity level k at design point x (corrected from page 163, [START_REF] Gratiet | Multi-fidelity Gaussian process regression for computer experiments[END_REF]) as:

σ 2 cont (k, x) = σ 2 δ,k (x) l-1 j=k ρ 2 j , (15)
which means that the variance contribution of the fidelity level k: (σ 2 δ,k ) is scaled using the recursive values of ρ j until we get to the highest fidelity l. These contributions are essential to build Sequential Design or Optimization strategies.

Recent works on linear regression multi-fidelity surrogates have been proposed in [START_REF] Zhang | Multifidelity Surrogate Based on Single Linear Regression[END_REF] where the accuracy is proved to be better than co-kriging [START_REF] Gratiet | Multi-fidelity Gaussian process regression for computer experiments[END_REF] but in the context of this research where we extend design optimization to multi-fidelity information sources, the access to an estimation of the variance is essential as it gives us an estimation of the error at each point of the design space. It thus helps us balance exploration/exploitation in the context of global optimization as explained in the following section.

A Python implementation of multi-fidelity co-kriging based on Le Gratiet's work can be found in the open source Surrogate Modeling Toolbox [START_REF] Bouhlel | A Python surrogate modeling framework with derivatives[END_REF] (SMT). *

III. MFEGO methodology

Bayesian optimization is defined by Močkus [START_REF] Močkus | On Bayesian Methods for Seeking the Extremum[END_REF] as an optimization technique based upon the minimization of the expected deviation from the exteremum of the studied function. The objective function is treated as a black-box function. A Bayesian strategy sees the objective as a random function and places a prior over it. The prior captures our beliefs about the behavior of the function. After gathering the function evaluations, which are treated as data, the prior is updated to form the posterior distribution over the objective function. The posterior distribution, in turn, is used to construct an acquisition function (often also referred to as infill sampling criterion) that determines what the next query point should be.

A. Efficient Global Optimization: EGO

We describe here the expected improvement as well as the EGO algorithm based on Jones et al. [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF]. Let F be an expensive black-box function to be minimized. We sample F at the different locations X = {x 1 , x 2 , . . . , x n } yielding the responses Y = {y 1 , y 2 , . . . , y n }. We build a kriging model (also called Gaussian process) with a mean function µ and a variance function σ 2 as presented in Section II. The next step is to compute the expected improvement (EI) criterion. To do this, let us denote:

f min = min{y 1 , y 2 , . . . , y n },
the expected improvement function (EI) can be expressed:

E[I(x)] = E[max( f min -Y, 0)] ( 16
) * https://www.github.com/SMTorg/smt
where Y is the random variable following the distribution N (µ(x), σ 2 (x)). By expressing the right-hand side of Eq. ( 16) as an integral, and applying some tedious integration by parts, one can express the expected improvement in closed form:

E[I(x)] = ( f min -µ(x))Φ f mi n -µ(x) σ(x) + σ(x)φ f mi n -µ(x) σ(x) , if σ > 0 0, if σ = 0 , (17)
where Φ(•) and φ(•) are respectively the cumulative and probability density functions of N (0, 1). Next, we determine our next sampling point as:

x n+1 = argmax x (E[I(x)]) . ( 18
)
We then test the response y n+1 of our black-box function F at x n+1 , rebuild the model taking into account the new information gained, and research the point of maximum expected improvement again. We summarize here the EGO strategy in Algorithm 1.

Algorithm 1 EGO Algorithm

f min ← min Y 5:
x i+1 ← argmax EI(mod, f min ) Choose x that maximizes EI 6:

y i+1 ← F(x i+1 )
Probe the function at most promising point x i+1 7:

X ← X ∪ {x i+1 } 8:
Y ← Y ∪ {y i+1 } Update the samples vectors 9:

i ← i + 1 10: f min ← min Y 11:
return f min This is the best known solution after n iter iterations

B. EGO extension to multi-fidelity: MFEGO

We extend the EGO algorithm to work with multiple fidelities. In practice as the SEGOMOE framework proposed by Bartoli et al. [11,[START_REF] Bartoli | An adaptive optimization strategy based on mixture of experts for wing aerodynamic design optimization[END_REF][START_REF] Bartoli | Adaptive modeling strategy for constrained global optimization with application to aerodynamic wing design[END_REF] is capable of handling both unconstrained and constrained problems and this with or without MOE, we choose to note in this paper the unconstrained version of SEGOMOE without MOE as EGO. The multi-fidelity is defined as MFEGO. Finally for constrained problems, we will use the notation SEGO for SEGOMOE without MOE and MFSEGO for its multi-fidelity adaptation.

The main idea of the proposed algorithm is that the search for the most promising sample and the choice of level of enrichment can be seen as problems to be tackled sequentially. Indeed, we can consider that given a Gaussian process (GP) or kriging model (m k GP(µ, σ 2 )) and a current best ( f min ), EI (or another infill sampling criterion) can be trusted to find the next most promising point. The choice of the fidelity level of enrichment is a different question that can be formulated thus: given the uncertainty at a chosen point, is it more interesting to query it at the highest-fidelity level or at lower ones? This choice of a two-stage decision process (fix most promising point then fidelity level of enrichment) offers the advantage of greatly reducing the time of infill sampling criterion optimization without being restrictive in any way. The choice of the fidelity level of enrichment translates an idea stating that one should favor the use of low-fidelity samples for exploration, and high-fidelity ones for exploitation, while making sure that there is no resampling at the same location.

Let f 0 , . . . , f l be the lowest-to highest-fidelity of a quantity of interest, with querying costs c 0 , . . . , c l . Using the recursive formulation [START_REF] Gratiet | Multi-fidelity Gaussian process regression for computer experiments[END_REF] with a constant ρ, we know that:

f k = ρ k-1 f k-1 + δ k for k ∈ {1, . . . , l} (19)
ρ k-1 = corr( f k , f k-1 ) std( f k ) std( f k-1 ) (20)
σ 2 k = ρ 2 k-1 σ 2 k-1 + σ 2 δ k (21)
where the notations are std(•) for the standard deviation and corr(•, •) for the correlation. Using the notations introduced in Eq. ( 14), the variance contribution of the fidelity level k at design point x * defined by Eq. ( 15) is recalled here:

σ 2 cont (k, x * ) = σ 2 δ,k (x * ) l-1 j=k ρ 2 j
Due to the necessity of nested DOEs, all lower fidelities are enriched at the same time, the uncertainty reduction becomes:

σ 2 red (k, x * ) = k i=0 σ 2 δ,i (x * ) l-1 j=i ρ 2 j
The corresponding cost to the enrichment of level of fidelities 0 through k is:

cost total (k) = k i=0 c i
We propose the level of enrichment criterion as follows:

t = argmax k ∈(0,...,l) σ 2 red (k, x * ) cost total (k) 2
where t is the highest fidelity level to be added (the nested DOE imposed to enrich all lower fidelities). It is questionable whether it is possible to always find a common unit of measurement for the cost of an observation and the variance reduction, but as the variance scales with the square of the correlation (Eqs. ( 20) and ( 21)), it is reasonable to penalize the uncertainty reduction by the square of the cost. The correlation translates the quantity of information shared between multiple functions. Note that Le Gratiet [START_REF] Gratiet | Multi-fidelity Gaussian process regression for computer experiments[END_REF] introduces a sequential design strategy where the uncertainty is simply penalized by the cost. Tests were realized with both approaches (penalization with simply the cost, or the square of the cost), the penalization with the square of the cost constantly gave better results.

C. Algorithm

The choice of level heuristic having been presented, we summarize the proposed strategy in Algorithm 2. We note {X 0 , . . . , X l } the DOEs of the fidelity levels 0 through l. {Y 0 , . . . , Y l } are the corresponding responses.

Algorithm 2 MFEGO Algorithm

crit k ← σ 2 r ed,k (x * , model)/( k i=0 costs[i]) 2 9: if crit k ≥ crit k-1 or σ 2 r ed,
f min ← min Y l 20:
x i+1 ← argmax EI(mod, f min ) Choose x that maximizes EI 21:

ENRICH_LEVEL(mod, x i+1 , costs)

22: i ← i + 1 23: f min ← min Y l 24:
return f min This is the best known HF solution after n iter iterations Note that the algorithm above does not take into account the constant cost of an iteration (surrogate model creation and maximization of EI). This can cause the optimization to last longer. A quick fix to take into account the constant cost of an iteration is to measure and to add that cost as an offset. This translates the idea that when we perform one iteration of the algorithm, the cost of an update is not only the querying time of the high-or low-fidelity code, but also the time needed to build the surrogate model and maximize the infill sampling criterion.

As can be seen in line 19 of the MFEGO algorithm, the value of f min is only updated by the best HF value. That is because other fidelity codes and datasets are only considered to help MFEGO and are not the objective of optimization. By only updating the best solution when a high fidelity sample is requested, the optimization is made more robust. Indeed, low fidelity can only be used to reduce some amount of the expected improvement, which is an amount due to uncertainty (exploration), rather than knowledge (exploitation). This property ensures that MFEGO converges to the global optimum of the high-fidelity function (in the same sense that EGO converges to the global optimum of a function). Note that the criterion integrates the correlation of the fidelity levels and thus makes the algorithm more robust to cases where there is a poor correlation between fidelity levels. Indeed, when ρ k → 0,

σ 2 red (k,x * )
cost total (k) 2 → 0, prompting the algorithm to move to the higher fidelities.

It is important that the algorithm will not resample high-fidelity points as the EI at these points is zero. It will also not resample low-fidelity points as the choice of level criterion will become zero for the low fidelity points, pushing the algorithm to query the high fidelity ones, thus enhancing the model and solution.

D. Illustration on 1-D analytic problem

The proposed strategy is illustrated on the following 1-D analytic problem [START_REF] Forrester | Multi-fidelity optimization via surrogate modelling[END_REF]:

f H F (x) = (6x -2) 2 sin(2(6x -2)) f LF (x) = 0.5 f H F + 10(x -0.5) -5
In Figure 1 we represent graphically these two fidelity levels. f H F is the expensive function and f LF is the cheap one. We make the assumption that the cost ratio between the fidelity levels is 1/1000. We start the optimization with 3 HF samples and 6 LF ones. We define the low-and high-fidelity datasets as:

X 0 = {x LF 1 , x LF 2 , x LF 3 , x LF 4 , x LF 5 , x LF 6 } and X 1 = {x H F 1 , x H F 2 , x H F 3 }.
The corresponding responses are Y 0 and Y 1 where

Y 0 = f LF (X 0 ) and Y 1 = f H F (X 1 )
. MFEGO algorithm is used to sequentially search through the space to minimize the maximum of EI. The algorithm has at each step two choices: 1) query once the low fidelity only, 2) query once the low fidelity and once the high fidelity if the EI cannot be reasonably reduced by a low-fidelity query. Figure 2 shows the reduction of EI resulting from cheap exploration, high-fidelity exploitation and model enhancement. After adding 4 LF points to explore the space and reduce EI (some low-fidelity points can be indistinguishable due to their proximity in the images above), the image (c) shows the local exploitation and enhancement of the model by querying the high-fidelity code once. The next HF sample finds the global optimum of the function. To summarize, we find the optimum after 7 iterations: 7 LF samples and 2 HF samples have been added, whereas the initial DOE was 6 LF points and 3 HF points (summary in Table 1). The classical mono-fidelity EGO approach with an initial DOE of 4 HF points required 11 additional HF points to find the global optimum, with a total cost three times higher than the MFEGO cost (15 compared To finish, we highlight that the choice of level criterion can be used with any exploration-versus-exploitation algorithm as the same idea can be applied. We use the low fidelity to reduce the uncertainty of the model and thus reduce the Exploration contribution to the infill sampling criterion. The highest fidelities are then used for the exploitation and effectively minimize the objective.

IV. Airfoil shape optimization

To validate MFEGO (resp. MFSEGO) and compare it to EGO (resp. SEGO) and later gradient-based approaches, three optimization problems have been considered: 1) 9-D unconstrained airfoil shape optimization 2) 15-D unconstrained airfoil shape optimization 3) 15-D airfoil shape optimization with 2 equality constraints A supervised (analytic) test case is also provided in Appendix (2-D Rosenbrock function).

A. Test case description

For the airfoil shape optimization test case, we used a parametrization based on a mode decomposition proposed by Li et al. [START_REF] Li | Data-based Approach for Fast Airfoil Analysis and Optimization[END_REF]. This decomposition is a Singular Value Decomposition (SVD) on the camber and thickness of an airfoil database. We use this parametrization to define an airfoil geometry for which we can compute characteristics such as the Lift Coefficient C l , Drag Coefficient C d , and Pitching Moment C m . The goal is to find the global optimal geometry of the airfoil as computed by the high-fidelity code ADflow. † ADflow has a Reynolds Averaged Navier-Stokes (RANS) multi-block flow solver developed in the MDO Lab (University of Michigan) that has been successfully applied to a variety of aerodynamic shape optimization problems [START_REF] Lyu | Aerodynamic Design Optimization Studies of a Blended-Wing-Body Aircraft[END_REF][START_REF] Lyu | Aerodynamic Shape Optimization Investigations of the Common Research Model Wing Benchmark[END_REF][START_REF] Yu | On the Influence of Optimization Algorithm and Starting Design on Wing Aerodynamic Shape Optimization[END_REF][START_REF] Bons | Multimodality in Aerodynamic Wing Design Optimization[END_REF].

To help MFEGO find the optimum at the lowest cost possible, we use a low-fidelity code Xfoil [START_REF] Drela | XFOIL: An analysis and design system for low Reynolds number airfoils[END_REF], which takes a fraction (1/200) of the HF code time to give an approximation of the result using hypotheses simplifying the problem (linearization, ignoring shocks, compression corrections...). Although Xfoil gives results that are different from ADflow, it has a strong correlation with the RANS solver. As they both represent the same physical problem, the same trends can be observed (Lift and Drag tend to increase when we increase the Angle of Attack for example). We use the correlation between the low and high fidelities to improve the accuracy of the model and perform faster optimizations.

ADflow compute derivatives efficiently thanks to a discrete adjoint method implementation [START_REF] Kenway | Effective Adjoint Approaches for Computational Fluid Dynamics[END_REF]. This has allowed us to have a gradient-based optimization reference to compare the Bayesian optimization approaches with, especially considering that 2D airfoil optimizations are unimodal [START_REF] Yu | On the Influence of Optimization Algorithm and Starting Design on Wing Aerodynamic Shape Optimization[END_REF].

To solve such optimization problems, work has been done on the noise estimation and regression/re-interpolation of the multi-fidelity surrogate model based on the work of Forrester et al. [START_REF] Forrester | Multi-fidelity optimization via surrogate modelling[END_REF]; when using multi-fidelity analyses, we can modify the interpolating co-kriging formulation in Eq. ( 4) such that each analysis can be regressed appropriately to filter any noise present in the data [START_REF] Forrester | Multi-fidelity optimization via surrogate modelling[END_REF]. In practice, this is done by adding a noise term λ to the diagonal of the covariance matrix and estimating it by maximum of likelihood [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF][START_REF] Rasmussen | Gaussian processes for machine learning[END_REF].

B. 9-D unconstrained optimization

The main characteristics of the first airfoil shape optimization problem are summarized in Table 2. The modes design variables have been normalized to fit the segment [0, 1]. The SVD gives as output basis vectors that can be scaled so that the design variable associated to a particular mode has a range of [0, 1]. In practice the scaled basis vectors were chosen not only to allow the reconstruction of the airfoil database used by Li et al. [START_REF] Li | Data-based Approach for Fast Airfoil Analysis and Optimization[END_REF], but also an additional margin was added to allow further exploration of new geometries. The Mach number (set to 0.25) and Reynolds number (set to 6 × 10 6 ) remain constant. With a cost ratio of 1/200, we obtain the results presented in Figure 3 and Table 3. The columns HF DOE and LF DOE in Table 3 denote the number of samples of each fidelity used to build the initial design of experiment (DOE). The columns HF Opt and LF Opt are respectively the number of high-fidelity (HF) and low-fidelity (LF) calls that the optimization routine used after that. The cost is a cost normalized by the HF cost, so that HF calls contribute 1 to the total cost, and LF calls contribute 1/200. It can be directly interpreted as the CPU time needed to reach the solution. We should note that EGO and MFEGO were stopped before 'convergence' (Bayesian optimization doesn't have a characterization of optimum, a budget of iterations is generally imposed as stopping criterion) and that at 'convergence', EGO and MFEGO should give the same optimum. But we can see from Table 3 that with lesser calls we are able to obtain a better best known solution. To understand the reason MFEGO performs better than EGO, let us define the series 'Gain' of a Bayesian optimization at iteration i as:

Gain i = |sol i -sol DOE |
where sol i is the best known solution after the i th iteration and sol doe is the best known solution resulting from the random sampling of the initial DOE. We use the Gain to compare the efficiency of the EGO and MFEGO algorithms. Figure 4 shows the Gain (improvement of the objective function compared to the start of the optimization) as a function of the cost. We see that EGO tends to have long plateaus where the objective does not improve, whereas MFEGO tends to improve at (almost) every HF call. This is due to the fact that EGO is an exploitation/exploration compromise. By using the lower-fidelity, MFEGO can explore the design space more cheaply, leaving HF calls for effective improvement of the objective. This is particularly important when increasing the dimension of the problem. As we increase the number of design variables, the space to explore becomes bigger. Using multiple fidelities allows MFEGO to scale better than EGO w.r.t the number of design variables.

C. 15-D unconstrained optimization

We showcase the better scaling by solving the same shape optimization problem with more design variables are described in Table 4 The Mach number (set to 0.25) and Reynolds number (set to 6 × 10 6 ) remain constant. As this particular problem is convex and unimodal, the local optimum found by SNOPT is also the global optimum of the function (red dotted line in Figure 5). We use this optimum as a reference to show the distance between the true optimum and the current solution of EGO and MFEGO. We can see from Figure 5 that the additional dimensions result in EGO spending a lot of time on the exploration of the space of the design variables. MFEGO is more immune to this trend as it can use the lower fidelity to explore the space before committing a HF computation. This observation is consistent across multiple runs. One way to make further use of the multi-fidelity kriging is to reduce the cost of the initial DOE. For example, instead of using 40 HF (for a total cost of 40), we use 16 HF and (approximately) 800 LF (for a total equivalent cost of (approximately) 20). We see in Figure 6 that by doing so, the space is better mapped, and that this improves greatly the speed of convergence at a much lower cost. For this problem we had access to the gradient-based optimum using SNOPT [START_REF] Gill | SNOPT: An SQP algorithm for large-scale constrained optimization[END_REF]. The dashed vertical blue line in Figure 6 separates the initial DOE building phase and the optimization driven by the Bayesian algorithm phase that comes after that. Figure 6 shows that MFEGO gives better results than EGO even with a less expensive initial DOE. The low fidelity, much cheaper, contributes some information to the surrogate model. MFEGO searched through the design space using 437 LF points and probed 8 HF points to attain the SNOPT solution (110.7 for SNOPT against 110.5 for MFEGO). We summarize the results in

D. Constrained optimization

We have shown that for unconstrained optimizations, MFEGO gives superior results to EGO due to the availability of cheap exploration of the design space. We show in what follows the advantages of multi-fidelity surrogate modeling for constrained optimization in an approach called MFSEGO. Not only the lower fidelities give access to cheaper exploration, but the overall model is more accurate resulting in a smaller constraints violation. First, let us introduce the measure root mean squared constraint violation (RMSCV) for an equality constraint as:

RMSCV = 1 N N j=1 val j -target 2
where N is the number of iterations of the algorithm, val j the value of the constraint function at iteration j and target is the target value of the constraint. A smaller RMSCV means that, on average, the algorithm respects the constraints more. We compare MFSEGO with SEGO and SNOPT through the following constrained optimization problem given in Table 6. The Mach number (set to 0.25) and Reynolds number (set to 6 × 10 6 ) remain constant. To test the robustness of the algorithm we add a Pitching Moment (C m ) constraint. Note that the AIAA Aerodynamic Design Optimization Discussion Group [START_REF]Aerodynamic Design Optimization Discussion Group Website[END_REF] proposes a pitching moment inequality constraint for the RAE2822 airfoil test case. We transformed this constraint into an equality constraint to make sure it is active on our problem. For a design to be accepted, we require that the constraint violation must be below 10 -3 for the absolute value of each of the two constraints.

In Table 7, the cost of SNOPT optimizations combines the number of direct problems solved [START_REF] Yu | On the Influence of Optimization Algorithm and Starting Design on Wing Aerodynamic Shape Optimization[END_REF], and the number of adjoint problems solved to compute derivatives ( Table 7 shows that MFSEGO is more robust to the addition of new constraints. This is due to the fact that the low fidelity helps in the surrogate modeling of the C m constraint as well. As a result, MFSEGO shows a 10 times lower constraint violation (RMSCV of 4.9e-3 for MFSEGO against 8.8e-2 for EGO). We should note that the small difference between SNOPT and MFSEGO in the objective should be attributed to the constraints tolerance. We should also point out that the first feasible solution found by EGO was after 52 HF iterations, whereas MFSEGO only needed 11 HF samples before finding its first feasible solution as shown in Figure 7 (vertical green dashed line). MFSEGO performed better than EGO and SNOPT in terms of cost. SNOPT needed 23 direct problems resolutions and 69 adjoint problems to find its optimal solution. In Figure 7 we show the evolution of the objective and constraint throughout the iterations of EGO and MFSEGO. The blue dashed line separates the initial DOE phase (random sampling) and the optimization phase. We show that MFSEGO has less constraints violations (less variations around the green dashed line denoting less constraint violations).

These three airfoil shape optimization problems demonstrated the good behaviour of the proposed multi-fidelity approach compared to classical Bayesian algorithms. Nevertheless, investigations will be carried on to increase both the Mach number, and the dimension of the design space.

V. Conclusion

In this paper, we formulate an extension of Bayesian optimization algorithms to work with multi-fidelity information sources. This approach is successfully coupled to EGO and SEGO to solve both unconstrained and constrained problems. In this approach, the number of high-fidelity calls is limited as the low fidelity is used to favor the exploration phase and thus costly high-fidelity simulation is only dedicated to the exploitation phase. We demonstrate the efficiency of Bayesian optimization in several test cases by enhancing the design space exploration part. For both the unconstrained or constrained optimization problems, MFEGO and MFSEGO found better results at a 50% lower computational cost compared to EGO. All problems have been treated within the SEGOMOE framework that is now extended with this new capability to handle multi-fidelity.

Future work will integrate the gradient information in the surrogate model and possibly in the devising of new criteria to allow for faster optimization. The extension to MDO workflow is also a target with studies on multi-fidelity coupling between two disciplines.

  Optimization (for unconstrained optimization) SEGO = Super Efficient Global Optimization (for constrained optimization) SEGOMOE = SEGO with Mixture of Experts MFEGO = Multi-Fidelity EGO MFSEGO = Multi-Fidelity SEGO

Fig. 1

 1 Fig. 1 Multi-fidelity of the 1-D analytic problem [7].

Fig. 2

 2 Fig. 2 Evolution of EI throughout MFEGO iterations on the 1-D analytic problem [7].

FunctionTable 2

 2 Definition of the 9-D unconstrained optimization problem.

Fig. 3

 3 Fig. 3 Comparison of EGO and MFEGO resulting airfoil shapes and C p distributions for an unconstrained L/D maximization with 9 design variables.

Fig. 4

 4 Fig. 4 Comparison of 'Gain' of EGO and MFEGO for unconstrained L/D maximization (9 design variables).

Fig. 5

 5 Fig. 5 Comparison of 'Gain' of EGO and MFEGO for unconstrained L/D maximization (15 Design Variables). The red dashed line is the optimum found by SNOPT.

Fig. 6

 6 Fig. 6 Comparison of evolution of objective as function of iterations for EGO and MFEGO for unconstrained L/D maximization with 15 design variables.

Fig. 7

 7 Fig. 7 Comparison of objective and constraints function of iterations for SEGO (a) and MFSEGO (b) for optimization of C d subject to C l and C m constraints with 15 design variables.

Table 1 Analytic problem optimization summary.

 1 to 5.013).

		HF DOE LF DOE HF Opt LF Opt Cost
	MFEGO	3	6	2	7	5.013
	EGO	4	-	11	-	15

  .

		Function/variable Description	Quantity	Range
	maximize	L/D	Lift-to-Drag ratio	1
	with respect to	α	Angle of attack	1 [0.0, 8.0] ( • )
		θ	Thickness modes	7	[0, 1]
		δ	Camber modes	7	[0, 1]
			Total variables	15

Table 4 Definition of the 15-D unconstrained optimization problem.

 4

Table 5 .

 5 HF DOE LF DOE HF Opt LF Opt Cost

	Obj

Table 5 L

 5 /D maximization

: Comparison of cost and objective for EGO and MFEGO unconstrained optimiza- tion with 15 design variables. The reference value is obtained by SNOPT.

Table 6 Definition of the 15-D constrained optimization problem.

 6

		Function/variable Description	Quantity	Range
	minimize	C d	Drag Coefficient	1
	with respect to	α	Angle of attack	1 [0.0, 8.0] ( • )
		θ	Thickness modes	7	[0, 1]
		δ	Camber modes	7	[0, 1]
			Total variables	15
	subject to	C L = 0.5	Lift coefficient	1
		C m = 0	Pitching moment	1
			Total constraints	2

Table 7 .

 7 23 × 3 = 69). The overall cost is then calculated as: cost SNOPT = time Direct + time Adjoint time Direct × N iter where time Direct is the time needed for solving the N iter direct problems, and time Adjoint is the time needed to solve N iter × N f uncs adjoint problems. N f uncs is the number of functions evaluated, here 3: C l for the objective, and C d and C m for the constraints. N iter was 23 in this case. This SNOPT cost is compared to the cost of SEGO and MFSEGO in HF DOE LF DOE HF Opt LF Opt

	Cost	Obj	Feasible RMSCV

Table 7 Comparison of SEGO and MFSEGO for constrained optimization with 15 design variables. The reference value is obtained by SNOPT.

 7

† https://github.com/mdolab/adflow

Acknowledgments

This work is part of the activities of ONERA -ISAE -ENAC joint research group and was partially supported by a ONERA internal project MUFIN dedicated to multi-fidelity. The authors would like to thank the ISAE-Supaero Foundation for its financial support and its role in making this project possible. The authors are also grateful to Rémi Lafage for his support with the SMT implementation and Rémy Priem for his support with the SEGOMOE framework. Mohamed Bouhlel and Joaquim Martins were funded by the Air Force Office of Scientific Research (AFOSR) MURI on "Managing multiple information sources of multi-physics systems," program officer Jean-Luc Cambier, award number FA9550-15-1-0038.

Appendix A. Rosenbrock optimization

We present here the results for the Rosenbrock function (See Figure 8) as it is a particularly challenging function to optimize (the function is extremely flat around the optimum). This optimization problem also allows us to compare MFEGO algorithm to some of the literature's state-of-the-art multiple information source optimization algorithms (not presented here, refer to [START_REF] Poloczek | Multi-Information Source Optimization[END_REF]). We set up the optimization to match the one presented in [START_REF] Poloczek | Multi-Information Source Optimization[END_REF] to have results that can be compared and to have a first order appreciation of the quality of the algorithm:

we assume a cost of 1000 for f H F and 1 for f LF .

The algorithm is started with a budget of 5 high-fidelity samples and 10 low-fidelity samples. Figure 9 shows MFEGO using the low fidelity to explore the space and more intensely the valley where the optimum is. It finally places a high-fidelity sample near the theoretical optimum of the function. It only needed one HF sample to do that.

This first result on a problem clearly cut for multiple information sources (assumption that only small variations exist between low and high fidelities) rather than multi-fidelities, was encouraging to carry on working on the algorithm and apply it to the aerodynamic shape optimization of an airfoil.


Franc ¸ois Coste
email: francois.coste@inria.fr

Deep learning languages: a key fundamental shift from probabilities to weights?

Recent successes in language modeling, notably with deep learning methods, coincide with a shift from probabilistic to weighted representations. We raise here the question of the importance of this evolution, in the light of the practical limitations of a classical and simple probabilistic modeling approach for the classification of protein sequences and in relation to the need for principled methods to learn nonprobabilistic models.

Introduction

Probabilistic models have been extremely important for successful language processing, from the first accomplishments in natural language processing (NLP), with n-grams and their successors (Russell and Norvig, 2010, chapters 22 and 23), to the tools routinely used today for the annotation of biological sequences in bioinformatics, based on the profile hidden Markov models (pHMM) or covariance models (CM) [START_REF] Durbin | Biological Sequence Analysis: Probabilistic Models of Proteins and Nucleic Acids[END_REF][START_REF] Franc ¸ois Coste | Learning the Language of Biological Sequences[END_REF]. Yet, after its success in computer vision and pattern recognition, deep learning is now replacing probabilistic models to become the new key player in the NLP field, showing that neural networks based on dense vector representations are very powerful for performing a large variety of NLP tasks. Impact of deep learning is until now limited in biology, but promising advances are being made (see for instance [START_REF] Ching | Opportunities and obstacles for deep learning in biology and medicine[END_REF], notably section 3 for the applications to biological sequences) and we can expect progress in NLP to be beneficial for biological sequence analysis once again.

While the main appealing feature of neural networks was initially non-linearity, the historical "S"-shape activation functions (such as the tanh function) are more and more replaced by the Rectified Linear Unit (ReLU) function, reducing nonlinearity to its simplest form: a function returning 0 or a weighted linear combination of its input values. The most important feature of deep learning in NLP seems now to be its ability to learn intermediate representations -such as the word vector representations (word embeddings) learned by word2vec [START_REF] Mikolov | Efficient estimation of word representations in vector space[END_REF] and GloVe [START_REF] Pennington | Glove: Global vectors for word representation[END_REF], or the recent deeper contextual representations learned by ELMo [START_REF] Peters | Deep contextualized word representations[END_REF], ULMFiT [START_REF] Howard | Universal language model fine-tuning for text classification[END_REF], and GPT (Radford et al., 2018a,b)-in a hierarchical manner up to the last layers of the neural networks in charge of their weighted combination for solving the task at hand. Using weighted combination of weight vectors (or matrices) optimized by gradient descent appears to be more and more the core of deep learning, as witnessed for instance by the growing importance of tensor-based libraries and of automatic differentiation [START_REF] Gunes Baydin | Automatic differentiation in machine learning: a survey[END_REF].

Hierarchical representations are common in formal grammars and gradient descent is not an original approach. The key advantage of deep learning with respect to the inference of classical grammatical representations could come from using and combining weights rather than probabilities. We think that better understanding and evaluating the importance of this difference is essential for future research in this area. As a first contribution, we present here practical problems illustrating the fundamental limitations of current probabilistic modeling approaches for language learning. The question we raise is whether they could be properly handled with more elaborate probabilistic modeling, or whether they demonstrate the interest of the deep "Probxit" towards weighted models.

To cope with inherent "noise", "variations" and "errors", languages have long been modelled with probabilistic models. In practice, these models rarely exceed the expressiveness of probabilistic context-free grammars, or even probabilistic automata. We introduce briefly the main related definitions and notations following Vidal et al. (2005a,b) to which we refer the reader for a more detailed presentation of the different probabilistic grammatical models of languages.

A probabilistic context-free grammar (PCFG) G is defined as a tuple Q, Σ, S, R, P where Q is a set of non-terminal symbols, Σ is a finite alphabet, S ∈ Q is the initial symbol, R is a set of rules A → α with α ∈ (Q ∪ Σ) * and P : R → [0, 1] define the probabilities of rules for each A ∈ Q, under the constraint:

∀A ∈ Q : (A→α)∈R P (A → α) = 1.
(1)

A successful (leftmost) derivation of a sequence x from S by a succession of rules r 1 . . . r l ∈ R l is denoted by S r 1 ...r l ⇒ x and its probability is P (r 1 . . . r l ) = l i=1 P (r i ). If a sequence x can be successfully derived from S, its probability is P (x) = r∈R + : S r ⇒x P (r), otherwise its probability is 0. Only consistent PCFGs, i.e. satisfying:

x∈Σ * P (x) = 1 (2)
defining thus a probability distribution on Σ * , are classically considered. Probabilistic regular grammars, more often depicted as probabilistic automata (PA), are PCFGs with simpler rules of the form A → aB, A → a, or A → λ with A, B ∈ Q, a ∈ Σ and λ denoting the empty sequence.

To give intuition, we will illustrate some limitations of purely probabilistic approaches on a practical task: modeling a family of protein sequences. This task is classically achieved with pHMMs, which are discrete HMMs with a predefined left-to-right topology. Illustration will focus on the core of this approach, simplified and transposed to PA to remain in the grammatical formalism. As an example of a family of protein sequences, we consider the SH3 domain, a region fairly well conserved among several protein sequences that is known to be important for binding and interacting with other proteins. This example has been realistically covered in the pHMM tutorial by Krogh (1998). We will simply use here a toy "profile PA" built from an excerpt of the alignment of the SH3 regions from some protein sequences, shown in figure 1. The choice of this example was obviously driven by our own domain of expertise, but also because it enables us to illustrate modeling problems with respect to "objective" physico-chemical considerations. The problems presented in the following sections should still be generic enough to be easily transposed to other applications and methods using probabilistic sequence models.

Limitations

We present here a list of practical limitations of classical probabilistic language modeling that should not arise when using weighted models. This does not mean that more elaborate probabilistic approaches could not solve these issues, nor that the best way to learn weighted models overcoming these limitations is known. We see these problems as critical testbeds for studying and comparing both approaches in order to better measure the contribution of the shift from probabilities to weights in the success of deep learning. Ideally, we would like these examples to also help developing principled approaches for learning weighted grammatical models, such as those available for probabilistic models. A pragmatic first step in this direction could be the elaboration of well-founded schemes for parameter estimation of "profile weighted automata" outperforming the simple maximum likelihood estimation of pHMMs' probabilities (or, better, the sophisticated Bayesian inference methods used to boost their performances by the addition of pertinent pseudocounts from Dirichlet mixture priors, as initiated for instance by [START_REF] Sjölander | Dirichlet mixtures: a method for improved detection of weak but significant protein sequence homology[END_REF]).

L1: Probability scattering at each nondeterministic derivation step Comparing the probabilities given by probabilistic grammars to sequences of different lengths is a well-know practical problem. The profile PA in figure 1 will assign smaller probabilities to longer proteins even if they have exactly the same SH3 region, because of the loss of probability in self-loop transitions consuming the rest of the protein sequence. This problem is tackled in bioinformatics by looking at the ratio between the probability of the sequence and its probability according to a null model: these Here, probability of amino-acids in left-to-right transitions is their observed frequency in the modelled column while they are assumed equiprobable in the self-loop transitions introduced at beginning and end of the PA, to parse complete proteins, and in its middle, to allow insertions between columns 4 and 5 (as in protein O45009).

ratios are usually directly integrated in the profiles by associating log-odds scores to transitions rather than probabilities. Profile are thus already used under the form of weighted models which combine two probabilistic models. More generally for PCFGs, equation 1 implies a probability loss by a factor strictly smaller than one, at each derivation step involving a nondeterministic rule. Such derivations scatter more and more probabilities between alternative suffix languages, resulting in negligible probabilities for longer derivations. Moreover, they induce a dependence of probability values to the number of non-deterministic derivations, which has often to be corrected or compensated in practical tasks, while this could be avoided by not satisfying equation 1.

L2: Dependence to number of choices for mass probability distribution within alternative rules Another consequence of equation 1 is that the probability mass has to be shared between all alternative rules. Its amount by rule depends then from the number of alternatives. For the characterization of SH3 regions, figure 1 shows the importance of having the amino-acid G at the first position. But when we look at the second conserved column, we can see that it contains only hydrophobic amino-acids and it might be as important for SH3 function to have a hydrophobic residue at the second position as a G at the first position. This cannot be modeled with PCFGs using the amino-acid alphabet: assuming that there are 14 hydrophobic acids and that they are equiprobable, the probability of transition with these amino-acids will be at most 1 14 making a small difference with the completely equiprobable 1 20 null model, in contrast to 14 15 , or any value close to 1, that can be assigned at the first position to a transition with G. We could imagine working on a larger alphabets, such as the amino-acids powerset, but the exponential growth of the alphabet size would then complicate training. Considering weights (e.g. reflecting adequacy of particular amino-acids at position, instead of occurrence) could be a more efficient solution in practice.

L3: Identical probability mass for all choice points Because of the normalisation to 1 in equation 1, all choice points have the same total weight in final probability. Yet, some choices can be more important than others. For instance, prior expert knowledge could tell us that amino acid P in the fourth position is more important for the function of SH3 than amino-acid G in the first position, but this could not be modelled with a PCFG. Authorizing a normalisation to variable values (e.g. reflecting relevance of position modelled with a particular non-terminal) would help here. Note that it could also be an indirect way to overcome limitation L2.

4 Towards language modeling Probxit?

Analysis of limitations Limitations seen above are induced by equation 1 which enables proper PCFGs to satisfy equation 2 and thus define a probability distribution over all sequences. While discarding the constraint of equation 1 is tempting, this constraint is intrinsically needed to distribute the remaining probability mass between all sequences that can be generated after each choice point, and it makes perfect sense to model the probability of sequence occurrence in the language with respect to the uncertainty brought by the non-deterministic rules.

Actually, what is illustrated here is more a misused feature of probabilistic grammars rather than a real limitation of these models: probability of sequences under constraint of equation 1 specifies the likelihood of their occurrence inside the language, not their probability of being in the language. Considering the SH3 example, probabilities tell us how often we can expect to see a given SH3 sequence comparatively to other SH3 sequences, not to estimate the uncertainty that the sequence can or not perform the SH3 function. If we are interested in membership to language, it is clearly given by the non-probabilistic part of the grammar: a sequence w is in the language L if it can be successfully derived by rules which have non-zero probability. If we were interested in the uncertainty of membership, we should consider learning a probability P (x ∈ L), with P (x ∈ L) + P (x ∈ L) = 1 instead of equation 2. Likelihood of occurrence in the language can thus be considered at most as an imperfect proxy for membership prediction.

Research directions

The first research direction to overcome the enumerated limitations is to work on learning the topology of the grammars. Note that this could be completed by learning also a membership probability, but the meaning of this uncertainty, and how it could be learned, would have to be clarified. It could be a measure of the uncertainty of decision with respect to current knowledge (e.g. based on the number of available examples supporting the decision). In the example of SH3, we can also imagine, for instance, that the binding affinity of different SH3 regions could influence their probability of binding to another protein and to perform their function: this could be formalized as a kind of membership probability, but would require information about the probability of performing the function in the training set. We will not consider these extensions here and will continue to focus on membership problem.

Learning a good grammar topology can be hard. In contrast, adding probabilities (and a recognition threshold) to simpler topologies has been shown successful to improve sufficiently their expressiveness for many classification tasks, while maintaining the number of parameters low enough to be estimable from training sets. This strategy has been shown efficient in practice, despite the limitations that we have identified. A second research direction is to study if more elaborate probabilistic approaches could be used more successfully here. The third research direction is to get rid of the constraint of equation 1, which is not really needed for membership prediction, to learn this way weights instead of probabilities of grammar rules.

In this third approach, the weights of sequences could be divided by a partition function (if it exists) to be considered as probabilities taking values in the interval [0, 1] and satisfy equation 2. Besides the problem of the existence and the computation cost of such a partition function, we think that the constraint of equation 2 has also to be discarded, since it restricts the expressiveness of the approach. Indeed, (positive) weights are useful if a threshold greater than zero is chosen so as sequences accepted in the language are those whose weight is above this threshold (otherwise, successful derivation by transitions with non-zero weight is sufficient). By transforming weights into probabilities, the choice of a threshold t would then only enable to define finite languages (with at most 1 t sequences since total probability mass is 1), restricting thus unnecessarily the languages that can be represented.

From probabilities to weights Getting rid of both equations means shifting from probabilistic to weighted representations. This move is visible in research on learning weighted grammars (see for instance contrastive estimation from [START_REF] Smith | Contrastive estimation: Training log-linear models on unlabeled data[END_REF], still motivated by mass distribution constraint) and can be considered an intrinsic feature of deep learning approaches. While practical successes show the interest of these approaches, better understanding their fundamental contribution and being able to develop principled approaches for learning weighted grammars, such as those developed for probabilistic grammars, remains a challenge. A key we identified is not (mis)using occurrence probabilities anymore and we exhibited some advantages of using weighted rather than probabilistic models, but theoretical and practical studies are still needed. The problem of training weighted, instead of probabilistic, profile automata on proteins seems a good support for these studies since it is simple and enables objective evaluation of the approaches. To address these questions, we think that we will need first to find the equivalent for weighted grammars of maximum likelihood or Bayesian approaches (especially for the incorporation of prior knowledge, an essential feature of using pHMMs on proteins) and the nice compositional properties of probabilistic approaches, and would really like to discuss this during the workshop.

Figure 1

 1 Figure 1: pHMM-like probabilistic automata (right) from a toy alignment of 15 SH3 sequences with 7 conserved columns (left). Height of amino-acid letters is proportional to the probability of the transition labeled by the letter.Here, probability of amino-acids in left-to-right transitions is their observed frequency in the modelled column while they are assumed equiprobable in the self-loop transitions introduced at beginning and end of the PA, to parse complete proteins, and in its middle, to allow insertions between columns 4 and 5 (as in protein O45009).

Acknowledgments

I am grateful to Hugo Talibart and Witold Dyrka for fruitful discussions that motivated this work submitted as a four pages opinion paper to ACL 2019 Workshop "Deep Learning and Formal Languages: Building Bridges".


Lucie Rouleau
email: lucie.rouleau@lecnam.net

Boris Lossouarn
email: boris.lossouarn@lecnam.net

Robin Darleux
email: robin.darleux@lecnam.net

Comparison of damping performances of constrained viscoelastic layers and passive piezoelectric networks

Keywords: Dynamic response, Viscoelastic damping, Piezoelectric damping I-INCE Classification of Subject Number:47

This work aims at comparing the damping performances of two passive treatments based on piezoelectric or viscoelastic patches. The motivation for such a comparison stems from the fact that the two damping treatments have been developed fairly independently, and are rarely compared. Firstly, the dynamic response of a simply supported metallic plate, equipped with constrained viscoelastic patches or piezoelectric patches connected to an electrical network, is measured experimentally. Secondly, a numerical model of the structure is set up and validated to evaluate the damping performances of both passive treatments under different configurations (for instance equal-mass and equal-thickness configurations). Finally, with regard to this work, the advantages and the limitations in using viscoelastic or piezoelectric treatments are discussed.

INTRODUCTION

In the context of noise and vibration mitigation, several damping technologies are proposed in the literature. Among them, purely passive treatments generally lead to low-cost and robust structural vibration control. These can be achieved by bonding constrained viscoelastic patches to the structure or by using piezoelectric patches connected to an electrical network. With piezoelectric patches, the vibrational energy is transfered to an electric circuit and dissipated in resistive components ; while with viscoelastic patches, it is transfered to the viscoelastic layer (undergoing shear deformations) and is dissipated in the damping material. Resonant piezoelectric shunts are commonly tuned to a single mechanical natural frequency to be controlled [START_REF] Thomas | Performance of piezoelectric shunts for vibration reduction[END_REF]. However, broadband vibration reduction can be obtained with multimodal damping, where the structure is coupled to a multi-resonant electrical network [START_REF] Lossouarn | Multimodal vibration damping of a plate by piezoelectric coupling to its analogous electrical network[END_REF]. On the other hand, due to the frequency dependent damping properties of viscoelastic materials, viscoelastic patches naturally provide broadband damping [START_REF] Rouleau | Characterization and modeling of the viscoelastic behavior of a self-adhesive rubber using dynamic mechanical analysis tests[END_REF]. Both passive damping treatments have been developed fairly independently and are rarely compared. The goal of this work is to provide an experimental and a numerical comparison of the damping performances of constrained viscoelastic layers and passive piezoelectric networks. A simply supported metallic plate, equipped with constrained viscoelastic patches or piezoelectric patches connected to an electrical network is under study. The next section describes the experimental set-up and reports the experimental results obtained using these panels. The efficiency of both passive treatments usually depends on the volume of material used, so that the performance is generally limited by weight and size constraints. Therefore, numerical models of the damped panels are implemented and validated in order to test the two damping technologies under different configurations, such as equal-mass and equal thickness configurations. Finally, with regard to this work, the advantages and the limitations in using viscoelastic or piezoelectric treatments are discussed.

EXPERIMENTAL RESULTS

This section describes the geometry of the panels equipped with viscoelastic or piezoelectric patches and the experimental set-up used to measure the structural dynamic response. Experimental results are reported and compared.

Experimental set-up

The structure under study is an aluminium panel, of dimensions 420 × 360 × 3 mm 3 . 42 viscoelastic or piezoelectric patches are periodically distributed on the panel, as indicated in Figure 1. Each patch is of dimensions 50 × 50 mm 2 . Viscoelastic patches are composed of a 120 µm viscoelastic layer (Smacwrap R ) and a 1.04 mm composite laminate which acts as a constraining layer. Piezoelectric patches are made of PIC 153 (PI Ceramic) and are 0.5 mm thick. A vacuum bonding technique is used to assemble all the layers. These damping solutions lead to a 15% added mass for the viscoelastic patches and a 33% added mass for the piezoelectric patches, not taking into account the mass of the electrical components. The electrical network that interconnects the piezoelectric patches consists of 42 identical unit cells made of passive components such as inductors and transformers. Their values and interconnections are chosen so as to create, in the electrical domain, the analogue of the plate to be controlled [START_REF] Lossouarn | Design of a passive electrical analogue for piezoelectric damping of a plate[END_REF]. This generates a multi-resonant network whose natural frequencies are tuned to that of the mechanical plate in order to obtain the equivalent of a multimodal tuned mass damper.

Panels are mounted in a specific frame to approximate simply supported boundary conditions [START_REF] Robin | A plane and thin panel with representative simply supported boundary conditions for laboratory vibroacoustic tests[END_REF][START_REF] O Robin | Setting up plane and thin panels with representative simply supported boundary conditions: comparative results and applications in three laboratories[END_REF]. The whole structure is suspended and excited by a shaker, while the structural dynamic response is measured by a laser vibrometer. The location of the shaker is given by the coordinates x = 0.12 m, y = 0.15 m and z = 0.003 m, and a pseudorandom excitation is generated. The frequency response of the structure is measured at the excitation point, in the frequency range 0 -900 Hz, with a frequency step of 0.156 Hz.

Experimental comparison

Following the experimental procedure described in the previous section, three panels are tested: a bare aluminium panel, an aluminium panel with viscoelastic patches and an aluminium panel with piezoelectric patches connected to a multi-resonant network. The measured frequency response functions are plotted in Figure 3. Both damping treatments induce a drastic reduction of the dynamic response over the whole frequency range. Multiresonant piezoelectric damping is more efficient than viscoelastic damping for the first

Material

Young vibration modes. However, they offer similar damping performances above 400 Hz.

For practical reasons, the configurations considered for piezoelectric and viscoelastic damping do not correspond to equal-mass or equal-thickness configurations. In order to compare the two passive treatments in a more objective manner, numerical models of the damped panels are implemented and experimentally validated.

VALIDATION OF THE NUMERICAL MODELS

This section describes the numerical models used to predict the dynamic response of a panel with piezoelectric or viscoelastic patches. Both models are validated by the experimental results presented in the last section.

Numerical model for the panel with viscoelastic patches

A finite element model of the panel with viscoelastic patches is implemented [START_REF] Rouleau | Interface finite elements for the modelling of constrained viscoelastic layers[END_REF]. Each layer is meshed by 20-node hexahedron elements, so that the shear behaviour of the viscoelastic layer is properly represented. The material properties considered in the calculation are reported in Table 1, and the master curves of the viscoelastic layer are plotted in Figure 4. The frequency dependent properties of the viscoelastic layer are taken into account in the finite element model through tabular data.

Figure 5 compares the frequency response function computed with the implemented finite element model to the measured response on the panel equipped with viscoelastic patches. Numerical results are in very good agreement with the experimental results, which validates our finite element model. The peak at 685 Hz corresponds to a flexural mode of the frame, which explains why the applied damping treatments have little influence on the amplitude of the response at this frequency.

Numerical model for the panel with piezoelectric patches

The numerical model used to predict the electromechanical response of the panel is based on a finite difference scheme with 9 elements for each of the 42 unit cells. Details about the numerical model of the piezoelectric plate coupled to its electrical analogue can be found in [START_REF] Lossouarn | Multimodal vibration damping of a plate by piezoelectric coupling to its analogous electrical network[END_REF]. The frequency response function of the panel with piezoelectric patches computed with this model is compared to the measured response in Figure 6. Again, very good agreement is observed between numerical results and experimental results. Both numerical models can now be exploited to study various configurations, and compare the damping performances of constrained viscoelastic layers and passive piezoelectric networks.

COMPARISON AND DISCUSSION

Weight and dimensions are classical constraints in the design of a damping treatment. Therefore, this section is dedicated to the numerical study of two configurations : equal-mass and equal-thickness. In order to obtain an equal-mass configuration for the viscoelastic and the piezoelectric treatments (the mass of the electrical circuit is not taken into account), the thickness of the piezoelectric patches is modified in the numerical model: 0.23 mm instead of 0.5 mm. It should be noted that the electrical network has been tuned accordingly in order to maximize the vibration damping of the second mode of the plate. In this way, both viscoelastic and piezoelectric patches have a total mass of 182 g (to be compared to the mass of the bare aluminium panel: 1.2 kg). The computed frequency response functions for both panels are compared in Figure 7. Results show that below 400 Hz, the damping performances of the piezoelectric patches are better than those of the viscoelastic patches. Above 400 Hz, the inverse tendency is observed.

The equal-thickness configuration is obtained by modifying the thickness of the viscoelastic patches in the numerical model: 41 µm for the viscoelastic layer and 0.359 mm for the constraining layer. Then, by considering 0.1 mm for each adhesive layer, viscoelastic patches and piezoelectric patches have the same thickness: 0.6 mm. The computed frequency response functions for both panels are compared in Figure 8. This configuration identifies the piezoelectric solution as the most efficient damping technique in terms of vibration damping, over the whole frequency range.

The results of this comparative study suggests that design constraints clearly have an impact on the relative damping performance of viscoelastic and piezoelectric patches.

CONCLUSION

The goal of this work was to compare the damping performances of two passive treatments based on constrained viscoelastic layers or on piezoelectric patches connected to a multi-resonant electrical network. To this end, experimental and numerical tests have been carried out on a metallic panel, equipped with viscoelastic or piezoelectric patches. Numerical models of both panels have been experimentally validated, which enables the numerical study of two configurations: equal-mass and equal-thickness configurations. The main conclusions of this work are :

-The damping performances of piezoelectric networks are generally more advantageous than those of viscoelastic materials to damped the first vibration modes.

-To damp higher modes of vibrations, the most efficient damping technique depends on the considered mass and dimension constraints.

The impact of other design constraints, such as cost, integrability or temperature, should be further investigated for an objective comparison of the damping performances of both passive treatments. For instance, when piezoelectric networks and constrained viscoelastic layers offer similar performances, the use of viscoelastic materials is preferable if the operating temperature remains constant, as it constitutes a low-cost solution with high integrability.

Figure 1 :

 1 Figure 1: Dimensions and placement of patches on the aluminium panel.

Figure 2 :

 2 Figure 2: Experimental set-up for the panel with viscoelastic patches (a) and the panel with piezoelectric patches connected to an electrical network (b).

Figure 3 :

 3 Figure 3: Frequency response functions measured on the bare aluminium panel, the panel with viscoelastic patches and the panel with piezoelectric patches.

Figure 4 :

 4 Figure 4: Master curves of the viscoelastic material (Smacwrap R ) at the reference temperature T = 20 o C.

Figure 5 :

 5 Figure 5: Experimental and numerical comparison of the frequency response functions of the panel with viscoelastic patches.

Figure 6 :

 6 Figure 6: Experimental and numerical comparison of the frequency response functions of the panel with piezoelectric patches connected to a multi-resonant network.

Figure 7 :

 7 Figure 7: Comparison of frequency response functions in an equal-mass configuration.

Figure 8 :

 8 Figure 8: Comparison of frequency response functions in an equal-thickness configuration.

Table 1 :

 1 modulus [Pa] Poisson ratio Density [kg/m 3 ] Material properties considered in the finite element model of the sandwich panel.

	Aluminium	6.9 10 10	0.346	2700
	Viscoelastic (static)	2.98 10 6	0.49	867
	Laminate composite	5.94 10 10	0.273	1580

ACKNOWLEDGEMENTS

We acknowledge gratefully Safran Aircraft Engines for their contribution to the realization of the panel with viscoelastic patches.


Robin Darleux
email: robin.darleux@lecnam.net

Boris Lossouarn

Jean-François Deü

Plate electrical analogue for multimodal energy transfer

This work presents the implementation of a plate electrical analogue and its validation by comparison between measurements and simulations. We also highlight the possible use of this network as an energy harvester.

I. INTRODUCTION

Bonding piezoelectric patches on a vibrating structure creates an electromechanical coupling that enables the transfer of the vibration energy to an electrical circuit. The extracted energy can then be either dissipated in a resistive component, or recovered through suitable power electronics. In any case, the design of the electrical circuit has a significant influence on the energy conversion performance. More specifically, piezoelectric shunt damping has been studied to damp one [START_REF] Thomas | Performance of piezoelectric shunts for vibration reduction[END_REF] or several [START_REF] Wu | Method for multiple mode piezoelectric shunting with single PZT transducer[END_REF] modes. These solutions might be effective to extract energy from a few resonances, however the complexity of the electrical circuit increases when considering a large number of modes. Other solutions including interconnected piezoelectric patches coupling the vibrating structure to its electrical analogue have been developed [START_REF] Alessandroni | A revival of electrical analogs for vibrating mechanical systems aimed to their efficient control by PZT actuators[END_REF][START_REF] Lossouarn | Design of a passive electrical analogue for piezoelectric damping of a plate[END_REF]. These solutions can help scavenging energy over a wider frequency range but require bonding several piezoelectric transducers on the mechanical structure.

In this paper, we consider a plate as the mechanical structure from which vibration energy should be extracted. Based on a direct electromechanical analogy [START_REF]The solution of partial differential equations by means of electrical networks[END_REF], the unit cell of the electrical network, analog to the unit cell of the mechanical structure, is derived from a discretization of the Kirchhoff-Love theory for plates. A model of the network is developed, in which the electronical damping induced by resistive components is taken into account. The computed currents distribution is then compared to measurements. The comparison validates the fully-passive electrical network that has been realized. Perspectives on further extension to energy harvesting are highlighted as well.

II. MODELING OF A PLATE ELECTRICAL ANALOGUE

A. Unit cell of the electrical network

To define the analogous electrical network of a plate, a discrete model of the mechanical structure is derived from a finite difference model of the Kirchhoff-Love equation in harmonic motion. We then use a direct electromechanical analogy [START_REF]The solution of partial differential equations by means of electrical networks[END_REF], which is summed up in TABLE 1. The discrete mechanical equations are then transformed in electrical equations, from which the equivalent unit cell of the analogous electrical network is identified. Complete details are given in [START_REF] Lossouarn | Design of a passive electrical analogue for piezoelectric damping of a plate[END_REF]. The unit cell of the electrical network we obtain is represented in FIGURE 1. The inductance , the capacitance and the transformer of ratio 1:1 represent the mass of the mechanical cell, the bending stiffness, and the coupling between the two directions of the mechanical structure, respectively. Transformers of ratio â/2:1 represent mechanical levers that appear in the finite difference model, where â is equivalent to the length of the square cell. Additional components, such as the internal resistances of the inductance and of the transformer windings, are taken into account as well.

B. Numerical implantation of the finite difference model

The unit cell previously defined is reproduced in two orthogonal directions to form the electrical network. Following an assembly process like in the finite element method, the electrical currents and voltages that appear in FIGURE 1 are then calculated in the case of external voltages applied at nodes of the network. The frequency response function (FRF) of the voltage drop across the inductors is computed, and will be compared to experimental results in the next section.

III. MEASUREMENTS ON A FULLY-PASSIVE PLATE ELECTRICAL ANALOGUE

A. Experimental setup

A fully passive network has been realized and is shown in FIGURE 3. Inductive components have been made following the method proposed in [START_REF] Lossouarn | Design of inductors with high inductance values for resonant piezoelectric damping[END_REF] and piezoelectric capacitances have been replaced by ceramic capacitors. A voltage excitation induces an electrical current distribution through the network. Those currents are measured by a differential voltage measurement across the inductors.

B. FRFs and currents distribution in the network

The FRF of the current flowing through the inductors divided by the excitation voltage amplitude has been measured at all 42 centers of the unit cells of the network. A comparison of the results obtained by simulation and measurements at one specific unit cell is shown in FIGURE 3. The measured operational current modes shapes, which mimic the vibration modes of the mechanical structure such as in FIGURE 4, can be compared to numerical results as well. These comparisons validate the electrical network, which has been developed as an electrical analogue of a mechanical plate. In this work, a plate electrical analogue dedicated to multimodal energy transfer is proposed. The analogue network is designed by discretizing the Kirchhoff-Love theory equations and using a direct electromechanical analogy. The unit cell we obtain is periodically repeated to create the network. Simulations and measurements of the currents distribution across the inductors of the network validate the concept of a plate electrical analogue. Now validated, the network would be able to extract vibration energy from the considered plate over a frequency range of nearly 1 kHz that includes 10 modes. The connection to the vibrating structure will be done thanks to piezoelectric patches, first for a broadband vibration damping purpose, and then for an energy harvesting purpose. As it stands, the main limitation to develop a piezoelectric energy harvester with such an electrical analogous network is that the inductors and transformers of the network are too resistive. The ensuing energy dissipation leads to damped resonances (see FIGURE 3), which limits the energy that could be recovered. The electrical components should thus be optimized as to reduce the internal resistances of the network, and so enable scavenging of the energy transferred from the mechanical structure.

FIGURE 2 .

 2 FIGURE 2. ELECTRICAL NETWORK ANALOGOUS TO A SIMPLY SUPPORTED PLATE. THE NETWORK IS MADE OF 6X7 IDENTICAL UNIT CELLS.

FIGURE 3 .

 3 FIGURE 3. COMPARISON BETWEEN EXPERIMENTS AND SIMULATIONS OF THE CURRENT FRF IN THE ELECTRICAL NETWORK

TABLE 1 .

 1 DIRECT ELECTROMECHANICAL ANALOGY

	Mechanical quantity	Electrical quantity
	Force	Voltage
	Velocity	Electrical current
	Compliance	Capacitance
	Mass	Inductance
	Viscous damping	Resistance

FIGURE 1 : UNIT CELL OF THE ELECTRICAL ANALOGUE OF A PLATE (B. LOSSOUARN ET AL. [4])

Excitation voltage

Unit cell of the network


Paul Baksic
email: p.baksic@unistra.fr

Hadrien Courtecuisse

Matthieu Chabanas

Bernard Bayle

FEM-based confidence assessment of non-rigid registration

N

on-rigid registration is often used for 3D representations during surgical procedures. It needs to provide good precision in order to guide the surgeon properly. We propose here a method that allows the computation of a local upper bound of the registration confidence over the whole organ volume. Using a biomechanical model, we apply tearing forces over the whole organ to compute the upper bound of the degrees of freedom left by the registrations constraints. Confrontation of our method with experimental data shows promising results to estimate the registration confidence. Indeed, the computed maximum error appears to be a real upper bound.

Introduction

While performing a surgery, the surgeon's goal is to minimize the risk for the patient. For this purpose, minimally invasive surgery has been favored over traditional open surgery, especially in abdominal surgery. Such procedures can be very challenging for the practitioner to perform. It is mainly due to the fact that they don't see what they are doing through their own eyes. Views are often showed through a screen and captured either with an endoscope, Ultrasound probe or other devices. This is challenging because it forces the surgeon to mentally map what he sees to what he is doing. A way to help surgeons during the surgery is to add virtual information on the screen, for instance the 3D volume of a tumor into a liver. This should provide information allowing him to resect a minimum of healthy tissues.

To reach this objective, one first needs to make a registration of the organ that the surgeon is seeing through the measurement device. Combining a biomechanical model to data extracted from intraoperative images have yielded good results. This can be done either by using some points extracted from the surface [START_REF] Plantefeve | Patient-specific biomechanical modeling for guidance during minimally-invasive hepatic surgery[END_REF][START_REF] Adagolodjo | Marker-based Registration for Large Deformations-Application to Open Liver Surgery[END_REF], the whole surface [START_REF] Peterlk | Fast elastic registration of soft tissues under large deformations[END_REF] or information on the organ volume extracted from ultrasounds [START_REF] Morin | Brain-shift compensation using intraoperative ultrasound and constraint-based biomechanical simulation[END_REF]. Yet all the above methods are subject to errors. Through this method, we propose a tool providing an estimation of the confidence of a model-based registration method.

Method

Our method takes as input the Finite Element mesh of the organ in the registered state. The registration can be performed by one of the method cited above or any method combining a biomechanical model with measured data. The method consists on evaluating the maximum mobility of the model at each point of its mesh, the confidence being high when the mobility of the point is low. The main assumption of the paper is that the deformation generated by the surgeon during the procedure does not damage tissues. Therefore, the maximum mobility is computed by applying an upper bound stress taken from the literature, known to cause irreversible deformations to the organ while satisfying the registration constraints into the simulation. This stress is successively applied along multiple directions on the model, in order to test each degree of freedom of the organ. Those directions are given by D = (x, -x, y, -y, z, -z). The computation of the mobility c at each point i is formalized in the following equation :

c i = max f (||q i -C(d × f × k)||) with d ∈ D (1)
Where C is a non linear function providing the positions of the model at static equilibrium after the application the volume force f multiplied by the surrounding volume k. The rest positions of the model are considered to be the positions provided by the registration method. This may lead to more mobility of the model because of the lack of internal stress which usually rigidify the tissue. This way, the computed displacement remains an upper bound. The confidence map can be represented as shown on the figure 1.

Confrontation to real data

Our method was tested on a data set of an ex-vivo liver of a lamb. External markers were put on the surface, which were then used to perform the registration along with internal markers that were used to evaluate the accuracy of the registration method. Then the volumes of the liver and the markers were segmented from CT-Scans in five different positions of the lamb liver. Cross-registrations were performed between the different positions of the liver using external markers. We used a linear elastic model of the deformation along with a corotational finite element formulation [START_REF] Mller | Stable real-time deformations[END_REF] of the problem on tetrahedron elements. This formulation shows good results and gives good approximations of the non-linearity of the deformation while keeping low time complexity. The liver finite element model parameters are given by a Young modulus of 6kP a and a Poisson ration of 0.499 as found in the literature for healthy livers [START_REF] Hudson | Inter-and intra-operator reliability and repeatability of shear wave elastography in the liver: a study in healthy volunteers[END_REF]. After the registrations, the distance between internal markers segmented from the CT images and the one provided by the model were considered as ground truth values of the registrations' accuracy. Our method was then applied on each configuration of the liver, in order to compute the theoretical upper bound of mobility. Given the material properties of the model described above, the force f used in the equation 1 is parameterized to reach 30% of strain. Indeed, despite the stress leading to tissue tearing vary amongst patients, this value is, according to the literature [START_REF] Brown | In-Vivo and Postmortem Biomechanics of Abdominal Organs Under Com-pressive Loads: Experimental Approach in a Laparoscopic Surgery Setup[END_REF][START_REF] Kerdok | Characterizing the nonlinear mechanical response of liver to surgical manipulation Doctoral dissertation[END_REF], identified as an upper-bound value leading to irreversible deformations of the tissue.

Results

The mechanical study performed by our method takes into account the complex coupling between constraints provided by image data. Indeed, as shown in figure 1, the resulting confidence map provided by our method is not directly related to the distance with external markers. Instead, it provides complex shapes that cannot be generated with simple geometric primitives. Figure 3 shows that our method predicted 96.6% of the time a real upper bound of the registration accuracy. Figure 3: Predicted error given by our method with respect to the measured error (given by the cross registrations) at each inner marker. Each label stands for a different target configuration.

Conclusion

We proposed a method allowing to evaluate the registration accuracy of internal structures. The method provides an upper bound uncertainty of positions of a biomechanical model, that can be used to discriminate and display only the reliable parts of the model in the augmented view. Future works concern additional validation study with more data and registration methods.

Figure 1 :

 1 Figure 1: High confidence (low mobility) areas are colored with shorter wavelength (blue). If the degree of confidence is lower than a threshold, the elements become transparent.

Figure 2 :

 2 Figure 2: Volume of a lamb liver with the outer markers in red and the inner ones in green.