Omar Gaci
email: omar.gaci@univ-lehavre.fr

Stefan Balev
email: stefan.balev@univ-lehavre.fr

Node Degree Distribution in Amino Acid Interaction Networks

Keywords: interaction network, protein structure, scale-free network

A protein interaction network is a graph whose vertices are the protein's amino acids and whose edges are the interactions between them. Using a graph theory approach, we study the properties of these networks. In particular, we are interested in the degree distribution and mean degree of the vertices. The results presented in this paper constitute the first steps of a new network approach to the protein folding problem.

I. INTRODUCTION

Proteins are biological macromolecules participating in the large majority of processes which govern organisms. The roles played by proteins are varied and complex. Certain proteins, called enzymes, act as catalysts and increase several orders of magnitude, with a remarkable specificity, the speed of multiple chemical reactions essential to the organism survival. Proteins are also used for storage and transport of small molecules or ions, control the passage of molecules through the cell membranes, etc. Hormones, which transmit information and allow the regulation of complex cellular processes, are also proteins.

Genome sequencing projects generate an ever increasing number of protein sequences. For example, the Human Genome Project has identified over 30,000 genes which may encode about 100,000 proteins. One of the first tasks when annotating a new genome, is to assign functions to the proteins produced by the genes. To fully understand the biological functions of proteins, the knowledge of their structure is essential.

In their natural environment, proteins adopt a native compact three-dimensional form. This process is called folding and is not fully understood. The process is a result of interactions between the protein's amino acids which form chemical bonds. In this paper we identify some of the properties of the network of interacting amino acids. We believe that understanding these networks can help to better understand the folding process.

In this study, we treat proteins as networks of interacting amino acid pairs [START_REF] Atilgan | Small-world communication of residues and significance for protein dynamics[END_REF]. In particular, we consider the subgraph induced by the set of amino acids participating in the secondary structure also called Secondary Structure Elements (SSE). We term this graph SSE interaction network (SSE-IN). We carry out a study to identify the node distribution relying on a dataset composed by more than 18000 proteins.

The rest of the paper is organized as follows. In section II we briefly present the main types of amino acid interactions which determine the protein structure. In section III we introduce our model of amino acid interaction networks. Section IV presents three general network models defined by their cumulative degree distribution. In section V we compare protein interaction networks to a general model and empirically characterize them based on a dataset. We show how the properties of these networks are related to the structure of the corresponding proteins. Finally, in section VI we conclude and give some future research directions.

II. PROTEIN STRUCTURE

Unlike other biological macromolecules (e.g., DNA), proteins have complex, irregular structures. They are built up by amino acids that are linked by peptide bonds to form a polypeptide chain. We distinguish four levels of protein structure:

• The amino acid sequence of a protein's polypeptide chain is called its primary or one-dimensional (1D) structure. It can be considered as a word over the 20letter amino acid alphabet.

• Different elements of the sequence form local regular secondary (2D) structures, such as α-helices or βstrands.

• The tertiary (3D) structure is formed by packing such structural elements into one or several compact globular units called domains.

• The final protein may contain several polypeptide chains arranged in a quaternary structure. By formation of such tertiary and quaternary structure, amino acids far apart in the sequence are brought close together to form functional regions (active sites). The reader can find more on protein structure in [START_REF] Branden | Introduction to protein structure[END_REF].

One of the general principles of protein structure is that hydrophobic residues prefer to be inside the protein contributing to form a hydrophobic core and a hydrophilic surface. To maintain a high residue density in the hydrophobic core, proteins adopt regular secondary structures that allow non covalent hydrogen-bond and hold a rigid and stable framework. There are two main classes of secondary structure elements (SSE), α-helices and β-sheets (see Fig. 1).

An α-helix adopts a right-handed helical conformation with 3.6 residues per turn with hydrogen bonds between C'=O group of residue n and NH group of residue n + 4.

A β-sheet is build up from a combination of several regions of the polypeptide chain where hydrogen bonds can form between C'=O groups of one β strand and another NH group parallel to the first strand. There are two kinds of β-sheet formations, anti-parallel β-sheets (in which the two strands run in opposite directions) and parallel sheets (in which the two strands run in the same direction).

III. AMINO ACID INTERACTION NETWORKS

The 3D structure of a protein is determined by the coordinates of its atoms. This information is available in Protein Data Bank (PDB) [START_REF] Berman | The protein data bank[END_REF], which regroups all experimentally solved protein structures. Using the coordinates of two atoms, one can compute the distance between them. We define the distance between two amino acids as the distance between their C α atoms. Considering the C α atom as a "center" of the amino acid is an approximation, but it works well enough for our purposes. Let us denote by N the number of amino acids in the protein. A contact map matrix is a N ×N 0-1 matrix, whose element (i, j) is one if there is a contact between amino acids i and j and zero otherwise. It provides useful information about the protein. For example, the secondary structure elements can be identified using this matrix. Indeed, α-helices spread along the main diagonal, while β-sheets appear as bands parallel or perpendicular to the main diagonal [START_REF] Ghosh | Dynamics of lysozyme structure network: probing the process of unfolding[END_REF]. There are different ways to define the contact between two amino acids. Our notion is based on spacial proximity, so that the contact map can consider non-covalent interactions. We say that two amino acids are in contact iff the distance between them is below a given threshold.A commonly used threshold is 7 Å [START_REF] Dokholyan | Topological determinants of protein folding[END_REF] and this is the value we use.

Consider a graph with N vertices (each vertex corresponds to an amino acid) and the contact map matrix as incidence matrix. It is called contact map graph. The contact map graph is an abstract description of the protein structure taking into account only the interactions between the amino acids. Now let us consider the subgraph induced by the set of amino acids participating in SSE. We call this graph SSE interaction network (SSE-IN) and this is the object we study in the present paper. The reason of ignoring the amino acids not participating in SSE is simple. Evolution tends to preserve the structural core of proteins composed from SSE. In the other hand, the loops (regions between SSE) are not so important to the structure and hence, are subject to more mutations. That is why homologous proteins tend to have relatively preserved structural cores and variable loop regions. Thus, the structure determining interactions are those between amino acids belonging to the same SSE on local level and between different SSEs on global level. Fig. 2 gives an example of a protein and its SSE-IN.

In [START_REF] Muppirala | A simple approach for protein structure discrimination based on the network pattern of conserved hydrophobic residues[END_REF], [START_REF] Brinda | A network representation of protein structures: implications for protein stability[END_REF] the authors rely on similar models of amino acid interaction networks to study some of their properties, in particular concerning the role played by certain nodes or comparing the graph to general interaction networks models. Thanks to this point of view the protein folding problem can be tackled by graph theory approaches.

As we will see in the next section, there are three main models of interaction networks, extensively studied and whose properties are identified. The purpose of our work is to identify specific properties which associate the proteins SSE-IN with a general network model. Based on such a pattern description of SSE-IN, one can plan the study of their formation, dynamics and evolution.

IV. INTERACTION NETWORKS

Many systems, both natural and artificial, can be represented by networks, that is by sites or vertices bound by links [START_REF] Strogatz | Exploring complex networks[END_REF]. The study of these networks is interdisciplinary because they appear in scientific fields like physics, biology, computer science or information technology. The purpose of these studies is to explain how elements interact inside the network and what are the general laws which govern the observed network properties.

From physics and computer science to biology and the social sciences, researchers have found that a broad variety of systems can be represented as networks, and that there log Pk log k Figure 3. Degree distribution for each the three models described by Amaral [START_REF] Amaral | Classes of small-world networks[END_REF]. The red line follows a power law, a function with a relatively "fat tail" as for scale-free networks. The green line corresponds to truncated scale-free networks because it describes a power law regime followed by a sharp cut-off. The black curve has a fast decaying tail, typically exponential, and corresponds to single-scale networks.

is much to be learned by studing these networks [START_REF] Broder | Graph structure in the Web[END_REF]. Indeed, the study of the Web [START_REF] Albert | The diameter of the world wide web[END_REF], of social networks [START_REF] Wasserman | Social network analysis : methods and applications , volume 8 of Structural analysis in the social sciences[END_REF] or of metabolic networks [START_REF] Jeong | The large-scale organization of metabolic networks[END_REF] are contribute to put in light common non-trivial properties to these networks which have a priori nothing in common. The ambition is to understand how the large networks are structured, how they evolve and what are the phenomenom acting on their constitution and formation [START_REF] Watts | Collective dynamics of 'smallworld' networks[END_REF].

One of the most important network properties is the degree distribution of vertices. A degree of a vertex is the number of edges incident to it. The mean degree of a network is the mean of the degrees of all vertices. For a network with n vertices and m edges the mean degree is z = 2m/n. We will note by p k the ratio of vertices having degree k (or the probability that a vertex has a degree k). The values p k define the degree distribution of a network. The cumulative degree distribution P k = ∞ i=k p k is the probability for a vertex to have a degree at least k.

The random graphs of Erdős and Rényi [START_REF] Erdõs | On random graphs I[END_REF], [START_REF] Erdõs | On the evolution of random graphs[END_REF] are the most studied network model. They have Poisson degree distribution. However, many real networks have different degree distributions. Amaral et al [START_REF] Amaral | Classes of small-world networks[END_REF] have studied networks that can be classified it three groups according to the shape of their cumulative degree distribution, see Fig. 3. First, scale-free networks are those with power law distribution 1) , a function which decreases polynomially with k. The second class are single scale networks with exponential degree distribution P k ∼ e -k/α . This distribution decreases exponentially, much faster than the previous. The third class are broad-scale or truncated scale-free networks with distribution

p k ∼ k -α or P k ∼ k -(α-
P k ∼ k -(α-1) e -k/α (1)
This distribution is somewhere between the previous two, a power law regime followed by a sharp exponential cutoff. The common feature of these classes is that most of the vertices have low degree and there exists a small number of high degree nodes. The last are called hubs and play important role for the connectivity of the whole network.

V. EXPERIMENTAL RESULTS

The first step before studying the proteins SSE-IN is to select them according to their SSE arrangements. Thus, a protein belongs to a SCOP fold level iff all its domains are the same. We have worked with the SCOP v1.73 classification. We have computed the measures from the previous section for the four mains classes of the hierarchical classification SCOP (see Tab. I). Thus, each class provides a broad sample guarantying more general results and avoiding fluctuations. Moreover, these four classes contain proteins of very different sizes, varying from several dozens to several thousands amino acids in SSE.

To compute the cumulative degree distribution for protein SSE-IN (denoted P k , see Eq. 1), we divided our dataset into two parts whose first one is composed with 20% of the total studied proteins. Then, we fit our specific sub dataset with a function expressed as follows

P k = a k -b exp -k/c
We have realised a numerical approximation using the method of least squares. Once, we have obtained the coefficients for our sub dataset, we apply the apply them for the rest of the studied protein SSE-IN, a sample of our results is presented on Fig. 4. We can remark that the curves describe a power law regime followed by the sharp cut-off. The power law function is expressed as following:

p(k) = 213.413k -α , where α = 3.2 ± 0.6
while the distribution is approximed by the next function:

P k = 1.48347k 0.962515 exp -k/2.12615
We observe the same result for all studied proteins, that is a cumulative degree distribution approximated by the function P k . Here, we discuss about characteristics or conditions which involve a such a behavior.

First, we are interested in the degree distribution and mainly its shape, see Fig. 5. We can see that degree distribution follows a Poisson distribution whose peak is reached for a degree near z. This result provides precision about how the vertices are connected within SSE-IN. It implies that the degree of the vertices is homogenous. In other words, a major part of them has a connectivity enough close to the mean degree. Consequently, the cumulative distribution depends on the mean degree value which acts as a threshold beyond which it decreases as an exponential since it's approximed via P k .

Second, we study how the mean degree evolves through all SSE-IN. Its distribution, see Fig. 6, indicates a relative weak variation according to the size. Even if two protein SSE-IN have size ratio around 10 or 100, their mean degree ratio is esimated to 1.05 or 1.15 and remains in the same scale order.

To illustrate the mean degree homogeneity we choose two proteins, namely 1SE9 and 1AON with sizes respectively 50 and 4998. Their size ratio is approximately 100. Even if the mean degrees are slightly different, the distributions are very similar (see Fig. 7).

To recapitulate, we show that the mean degree values constitute a threshold for protein SSE-IN cumulative degree distribution. For degrees lower than the mean degree it decreases slowly and after this threshold its decrease is fast compared to an exponential one, as shown Fig. 4,5,7.

Consequently, we find a way to approximate all proteins SSE-IN cumulative degree distribution by the function P k which can be adjusted. This function describes a power law regime followed by a sharp cut-off which arises for degree values exceeding the mean degree. Proteins SSE-IN are so truncated scale-free networks.

Since the mean degree plays the role of a threshold beyond which the cumulative degree distribution decreases exponentially, it is interesting to study its evolution with the size of the network. Fig. 6 shows that the mean degree increases very slightly with the size of the network. Even for networks with size ratio of 100, the mean degree ratio is only 1.15. As an example, see Fig. 7.

Whatever the size of the network is, we observe that the mean degree is always between 5 and 8. This mean degree interval is a common property characterizing all SSE-IN. In order to explain this property, let us consider the structure of our networks. They are composed of densely connected subgraphs corresponding to secondary structure elements (see Fig. 8). The number of edges connecting different subgraphs is relatively small, but these edges are the most important, since they correspond to interactions determining the tertiary structure.

We start by computing the mean degree in each SSE subgraph. The results are shown on Fig. 9. We can see that the mean degree evolution at microscopic level is almost the same as at macroscopic level (compare to Fig. 6). Independently of the SSE size and type, the mean degree of each SSE subgraph, z SSE is always bounded:

z min < z SSE < z max (2)
when the size of the network is more than 10. In the general case z min = 5 and z max = 8, but when we consider a specific SSE size and type, finer bounds can be found (see Fig. 9). Now let us consider a whole protein. Suppose that it contains s secondary structure elements and let the element i has n i vertices and m i edges, i = 1, . . . , s. Then the total number of vertices is n = s i=1 n i and the total number of edges is m = s i=1 m i + m inter , where m inter is the number of edges connecting vertices from different SSEs. Let r = m inter /m be the ratio of inter-SSE edges. Then:

m n = s i=1 m i + m inter s i=1 n i = s i=1 m i s i=1 n i + r m n (3)
and hence for the mean degree z we have On the other hand, from (2) it follows that

z = 2m n = 2 1 -r s i=1 m i s i=1 n i (4)
z min 2 n i < m i < z max 2 n i , i = 1, . . . , s (5)
By summing up the last equation we obtain

z min 2 < s i=1 m i s i=1 n i < z max 2 (6)
which together with (4) gives

z min 1 -r < z < z max 1 -r (7)
The last equation gives finer bounds on the mean degree. It shows that the bounds on z depend not only on the bounds on z SSE , but also on the ratio of inter-SSE edges. A higher proportion of inter-SSE edges shifts up the bounds. Proteins with bigger size have more SSEs and hence more links between different SSEs. This explains the increase of the mean degree with the size of the networks. Fig. 10 shows that the number of inter-SSE edges is quite variable, but it never exceeds 20%. It is the consequence of the excluded volume effect, since the number of residues that can physically reside within a given radius is limited. This last property explains why the mean degree is homogenous.

VI. CONCLUSION

In this paper we introduce the notion of interaction network of amino acids of a protein (SSE-IN) and study some of the properties of these networks. The main advantage of this model is that it allows to cope with different biological problems related to protein structure using graph theory tools. Ignoring details, such as the type and the exact position of each amino acid, this abstract and compact description allows to focus on the interactions' structure and organization.

In this paper we show that we can approximate all proteins SSE-IN cumulative degree distribution by a unique function. This function describes a power law regime followed by a sharp cut-off which arises for degree values exceeding the mean degree. Proteins SSE-IN are so truncated scalefree networks. This node distribution implies that there exist amino acids whose degree is marginal (greater than the mean degree).

A short term perspective is to study the hubs notably to understand how they appear in the folded protein. By this way, we would be able to understand if the nature of a hub depends on its position in the protein sequence or it depends on other chimical parameters.

As a long term perspective, the characterization we propose constitutes a first step of a new approach to the protein folding problem. The properties identified here, but also other properties we studied [START_REF] Gaci | Hubs identification in amino acids interaction networks[END_REF], [START_REF] Gaci | The small-world model for amino acid interaction networks[END_REF], can give us an insight on the folding process. They can be used to guide a folding simulation in the topological pathway from unfolded to folded state.

Figure 1 .

 1 Figure 1. Left: an α-helix illustrated as ribbon diagram, there are 3.6 residues per turn corresponding to 5.4 Å. Right: A β-sheet composed by three strands.

Figure 2 .

 2 Figure 2. Protein 1DTP (left) and its SSE-IN (right).

Figure 4 .Figure 5 .

 45 Figure 4. Cumulative degree distribution for protein 1COY SSE-IN.

Figure 6 .Figure 7 .

 67 Figure 6. Mean degree distribution according to protein SSE-IN size. It evolves with values enough close, between 5 and 8.

Figure 8 .

 8 Figure 8. SSE-IN of 1DTP protein. The edges connecting different SSE are green.

Figure 9 .

 9 Figure 9. SSE subgraphs size distribution and mean degree as a function of the size.

Figure 10 .

 10 Figure 10. Ratio of inter-SSE edges (r) as a function of the network size for the four classes of studied proteins.

Table I STRUCTURAL

 I FAMILIES STUDIED. WE CHOOSE ONLY FAMILIES WHICH CONTAIN MORE THAN 100 PROTEINS, FOR A TOTAL OF 18294 PROTEINS. WE HAVE WORKED WITH THE SCOP 1.73 CLASSIFICATION.

	Class	Number of Number of
		families	proteins
	All α 12	2968
	All β	17	6372
	α/β	18	5197
	α + β 16	3757


A Serani

R Pellegrini

J Wackers

C.-E Jeanson

P Queutey

M Visonneau

M Diez

Adaptive multi-fidelity sampling for CFD-based optimisation via radial basis function metamodels

The paper presents a study on four adaptive sampling methods of a multi-fidelity (MF) metamodel, based on stochastic radial basis functions (RBF), for global design optimisation based on expen-sive CFD computer simulations and adaptive grid refinement. The MF metamodel is built as the sum of a low-fidelity-trained metamodel and an error metamodel, based on the difference between high-and low-fidelity simulations. The MF metamodel is adaptively refined using dynamic sam-pling criteria, based on the prediction uncertainty in combination with the objective optimum and the computational cost of highand low-fidelity evaluations. The adaptive sampling methods are demonstrated by four analytical benchmark and two design optimisation problems, pertaining to the resistance reduction of a NACA hydrofoil and a destroyer-type vessel. The performance of the adaptive sampling methods is assessed via objective function convergence.

Introduction

Fluid-dynamic shape design of complex industrial systems like aerial, ground, and water-borne vehicles demands the use of high-fidel ty numerical solvers with large computational grids to assess accurately the design performance and make sound design decisions. The latter can be achieved by combining computational fluid dynamics (CFD) analysis with a shape/design modific tion tool (CAD) and a minimisation algorithm into an automatic simulationbased design optimisation (SBDO). The optimisation algorithm may require a large number of function evaluations to converge to the final solution, especially if a global optimum is desired. Therefore, the resulting computationalcostoftheSBDOprocesscouldbecome v e r yh i g h ,m a k i n gS B D Oh a r d l ya ff o r d a b l ef o rm o s t users and projects, for which computational resources and time are usually limited.

To reduce the computational cost of the SBDO process, metamodeling methods have been developed and successfully applied in several engineering field [START_REF] Viana | Special Section on Multidisciplinary Design Optimization: Metamodeling in Multidisciplinary Design Optimization: How Far Have We Really Come?[END_REF]). Among other metamodels, radial basis functions (RBF) methods have demonstrated their accuracy and efficie y in engineering design [START_REF] Jin | Comparative Studies of Metamodelling Techniques Under Multiple Modelling Criteria[END_REF]a l o n g w i t ht h e i re a s eo fi m p l e m e n t a t i o n .T h ep e r f o r m a n c e of metamodels is problem-dependent and determined by several concurrent issues, such as the presence of nonlinearities, the problem dimensionality, the noisy or smooth behaviour of the function, and the approach used for its training [START_REF] Liu | ASurveyof Adaptive Sampling for Global Metamodeling in Support of Simulation-based Complex Engineering Design[END_REF].

The research of accurate and effici t methods for metamodel-based analysis and optimisation has recently moved from standard (or static) to functionadaptive approaches, also known as dynamic metamodels. In addition to auto tuning, a dynamic metamodel is able to improve its fitting capability by adaptive sampling or training. The design of experiments (DoE) used for metamodel training is not defi ed a priori but dynamically updated, exploiting the information that becomes available during the analysis process. The purpose of performing an adaptive DoE is to add training points anywhere it i sm o s tu s e f u l ,s oa st ou s ear e l a t i v e l yl o wn u m b e r of function evaluations to represent the function. For global exploration of design/uncertainty spaces (such as in global optimisation and uncertainty quantification),adynamicRBF(DRBF)hasbeenformulatedby Volpi et a l. (2015), based on the maximum prediction uncertainty, and compared to dynamic Kriging [START_REF] Zhao | Metamodeling Method Using Dynamic Kriging for Design Optimization[END_REF]. Similarly, the non-linearity of the response and the largest distance among training-set sampleshavebeenusedbyMackmanandAllen (2010) for RBF. To improve the efficie y of global optimisation procedures, sequential and parallel infil crit e r i ah a v eb e e nu s e db a l a n c i n gt h ee x p l o r a t i o na n d exploitation of the design space: the expected improvement [START_REF] Jones | Effici t Global Optimization of Expensive Black-Box Functions[END_REF]h a sb e e n proposed to update Kriging model; multiple points have been added on Kriging model based on approximate computation of the probability of improvement [START_REF] Viana | S u r r o g a t e -B a s e d Optimization with Parallel Simulations using the Probabil-ityofImprovement[END_REF]; a multi-objective infil criteria based on the estimated variance and the response value has been proposed by Yi, Kwon, and Choi (2014); similarly, a parallel sampling criterion for DRBF has been proposed by [START_REF] Diez | Simulation-Based Design Optimization by Sequential Multi-criterion Adaptive Sampling and Dynamic Radial Basis Functions[END_REF] based on the identific tion of the global minimum of the objective function and maximum of the prediction uncertainty.

In addition to dynamic metamodels, multi-fidel ty (MF, or variable-fidel ty) approximation methods have been developed with the aim of combining the accuracy of high-fidel ty solvers/evaluations with the computational cost of low-fidel ty solvers/evaluations. Thus, MF metamodels are trained with a combination of high-fidelity (accurate, expensive) simulations and low-fidel ty (approximate, less expensive) simulations. Combining metamodeling methods with MF approximations potentially leads to a further reduction of the computational cost of the SBDO procedure. Additive and/or multiplicative correction methods, also known as 'bridge functions' or 'scaling functions' [START_REF] Han | Improving Variable-fidel ty Surrogate Modeling Via Gradient-enhanced Kriging and a Generalized Hybrid Bridge Function[END_REF], are used to build MF metamodels, using high-and low-fidel ty evaluations (Ng and Eldred 2012;Z h e n ge ta l .2013; Pellegrini et al. 2016). Several metamodels have been used in the literature with MF data, such as nonintrusive polynomial chaos (Ng and Eldred 2012), co-kriging [START_REF] Baar | Uncertainty Quantification for a Sailing Yacht Hull, Using Multi-fid lity Kriging[END_REF]a n dR B F( P e l l e g r i n i et al. 2016). In SBDO based on CFD computations, high-and low-fidel ty evaluations may be obtained by varying the physical model, the size of the computational grid, the computational time step, the convergence level of the simulations, and/or combining experimental data with numerical simulations [START_REF] Kuya | Multifidel ty Surrogate Modeling of Experimental and Computational Aerodynamic Data Sets[END_REF].

Adaptive sampling methods have been recently combined with MF metamodels: Huang et al. (2006) propose an augmented expected improvement for multi-fidel ty Kriging; Pellegrini et al. (2016)u s e a DRBF, sequentially adding low-fidel ty or both high-and low-fidel ty samples, depending on which fidelity has the greatest prediction uncertainty and on a parameter based on the computational cost; [START_REF] Liu | ASequentialSampling Strategy to Improve the Global Fidelity of Metamodels in Multi-level System Design[END_REF]p r o p o s eas a m p l i n gm e t h o db a s e do n the propagation of the low-fidelity variance to the MF metamodel, updating only the fidelity that provides the maximum reduction of the MF uncertainty; [START_REF] Cai | Adaptive Radial-Basis-Function-Based Multifidel ty Metamodeling for Expensive Black-Box Problems[END_REF] sample both high-and low-fidelity at each iteration, based on the cross-validation error and a Voronoi partition.

T h eo b j e c t i v eo ft h ep r e s e n tw o r ki st oa s s e s st h e performance of a multi-fidel ty DRBF (MF-DRBF) metamodel for global optimisation, conditional to four adaptive sampling methods. These are based on: (i) the maximum prediction uncertainty, (ii) a multi-fidel ty version of the expected improvement, (iii) the maximum prediction uncertainty and the objective function through an aggregated merit factor, (iv) multicriteria sampling based on the solution of a multiobjective problem considering the maximum prediction uncertainty and the objective function value [START_REF] Diez | Simulation-Based Design Optimization by Sequential Multi-criterion Adaptive Sampling and Dynamic Radial Basis Functions[END_REF].

The sampling methods for the MF-DRBF metamodel are demonstrated by four multi-modal analytical benchmark problems with one and two dimensions. Their performance is assessed by studying the convergence of the optimisation procedure and the normalised root mean square error of the prediction. Finally, the sampling methods are applied to two CFD-based design problems pertaining to the minimisation of (a) the drag coeffici t of a four-digit NACA hydrofoil at Reynolds number equal to 8.41 The CFD-solver fidel ty is varied by using two computationalgridlevels,definedbyanadaptivegridrefinement technique (Wackers et al. 2017). Optimisations are performed with a deterministic version of the particle swarm optimisation (DPSO, Serani et al. 2016b) algorithm.

Stochastic radial basis functions

Consider an objective function f (x),w h e r ex ∈ R N i st h ed e s i g nv a r i a b l ev e c t o ra n dN the design space dimension. Let the true function value be known in a number J of training points x j with associated objective function values f (x j ). The metamodel prediction f (x) is computed as the expected value (EV) of a prediction ensemble obtained considering a stochastic tuning parameter in the RBF kernel, e.g. τ ∼ unif[1, 3]:

f (x) = EV g (x, τ ) τ ,w i t h g (x, τ ) = J j=1 w j ||x -x j || τ ,( 1 )
where w j are unknown coeffici ts, || • || is the Euclidean norm. The coeffici ts w j are determined enforcing exact interpolation at the training points g(x j , τ) = f (x j ) by solving Aw = f,w i t hw ={w j }, a ij =||x ix j || τ and f ={f (x j )}.

The uncertainty U f (x) associated with the stochastic RBF metamodel prediction is quantifi d by the 95%-confidence interval of g(x, τ),e v a l u a t e du s i n ga Monte Carlo sampling over τ [START_REF] Volpi | Development and Validation of a Dynamic Metamodel Based on Stochastic Radial Basis Functions and Uncertainty Quantification[END_REF]. SimilarlytoKennedyandO'Hagan(2000), the MF prediction f (x) is define by an additive correction to a low-fidel ty trained metamodel f L (x) as

Multi-fidelity metamodel

f (x) = fL (x) +ε (x) ,( 2 )
where the correction is provided by the metamodel of theerr o r(o rdiscr epa ncy)ε(x),definedasthediff erence between high-(HF) and low-fid lity (LF) evaluations (f H and f L ):

ε (x) = f H (x) -f L (x) .( 3 )
The training set for fL is denoted by L,w h e r e a s the training set for ε is denoted by E ⊆ L. Assuming (for the sake of simplicity and ease of implementation) that the uncertainties associated with the low-fidelity and error metamodels (U f L and U ε respectively) are uncorrelated, the uncertainty of the MF prediction is

U f (x) = U 2 f L (x) + U 2 ε (x).( 4
) A sketch of the MF metamodel concept and notation is shown in Figure 1.N otetha t'∧'indica tesMF approximation and ' ∼ 'denotestheRBFprediction.

Adaptive sampling methods

The MF metamodel is dynamically updated by adding new training points following a two-step procedure (see Figure 2): (1) the coordinates of the new training point x ⋆ arechosenusingoneoftheadaptivesampling methods described in the following subsections; (2) once x ⋆ is identifie , either L or E are refined, based on the following statement

If U fL (x ⋆ ) ≥ βU ε(x ⋆ ),a d d {x ⋆ , f L (x ⋆ )} to L, else, add {x ⋆ , f L (x ⋆ )} to L and {x ⋆ , ε(x ⋆ )} to E, (5)
where β ∈ [0, 1) is the ratio between the LF and HF computational costs. In the fi st case, only a lowfidel ty evaluation is performed, whereas the second case requires both low-and high-fidel ty evaluations at the same point x ⋆ .

Maximum uncertainty

Maximum-uncertainty adaptive sampling (MUAS) hasbeendiscussedformulti-fidelityproblemsinPellegrini et al. ( 2016). This method identifie a new training point (see Figure 3) by solving the single-objective maximisation problem

x ⋆ = argmax x [U f (x)]. (6)

Multi-fidelity expected improvement

An extension of the expected improvement [START_REF] Jones | Effici t Global Optimization of Expensive Black-Box Functions[END_REF]f o rm u l t ifi d e l i t yR B F applications (MFEI) has been proposed by [START_REF] Pellegrini | Adaptive Sampling Criteria for Multifidel ty Metamodels in CFD-based Shape Optimization[END_REF]. It is define as

MFEI (x) = EV max f min -g L (x, τ) + g ǫ (x, τ),0 τ , with f min = min f (x) ,( 7
)
where g L (x, τ) and g ǫ (x, τ) are RBF interpolation of f L and ǫ,f o l l o w i n gE q .( 1 ) .T h e r e f o r e ,t h eM F E Ii s the expected value of potential reduction considering stochastic multi-fidel ty RBF and minima provided by f at the current iteration of the sampling process.

The MFEI adaptive sampling identifie a new training point (see Figure 4) by solving the single-objective maximisation problem

x ⋆ = argmax x [MFEI(x)]. ( 8
)

Aggregate criteria

An aggregate-criteria adaptive sampling (ACAS) method is defined with the aim of finding points with large uncertainty and small objective function value. Accordingly, ACAS identifie a new training point (see Figure 5) by solving the single-objective minimisation

x ⋆ = argmin x [ f (x) -U f (x)]. (9)

Multi-criteria

Similarly to the aggregate criteria, a multi-criteria adaptive sampling (MCAS, Yi, Kwon, and Choi 2014; [START_REF] Diez | Simulation-Based Design Optimization by Sequential Multi-criterion Adaptive Sampling and Dynamic Radial Basis Functions[END_REF])identifiesN P new training points (see Figure 6) considering both the MF prediction and the associated uncertainty (multi-criteria infil sampling). Specifical y, a multi-objective optimisation problem addressing the minimisation of the objective function and the maximisation of the prediction uncertainty is solved minimise f (x) and maximise U f (x),

subject to U f (x)>U ⋆ f , (10)
where

U ⋆ f = γ R,withR = sup[f H (x)] -inf[f H (x)] < ∞,
is a constraint parameter, defi ed in the view of the fact that: (a) sampling too close to available training points does not add useful information to the analysis, (b) as the distance between training points decreases, the matrix A may become ill-conditioned, and (c) the uncertainty at the training points is zero, i.e.

lim |x-x ⋆ |→0 U f (x) = 0 (11)
The non-dominated solution set obtained by the solution of Equation ( 10) is down-sampled in order to identify N P equally spaced points along the curvilinear coordinate define by the Pareto frontier (see Figures 6 and7). N P depends on the number of non-dominated solutionsandtypicalvaluesgofrom2to10,whereas γ is generally set equal to (0.1 ÷ 1)%R.

Optimisation problems

The assessment of the adaptive sampling methods for theMF-DRBFisbasedonfouranalyticalbenchmarks and two CFD-based design optimisation problems. A deterministic single-objective formulation of the particle swarm optimisation (DPSO) algorithm (Serani et al. 2016b), is used for the metamodel-based optimisations, as well as for the solution of the minimisation/maximisation sampling problems of Equations ( 6), (8), and (9). A multi-objective extension of DPSO as presented by [START_REF] Pellegrini | Formulation and Parameter Selection of Multi-objective Deterministic Particle Swarm for Simulation-based Optimization[END_REF]i su sed for the solution of Equation (10).

Analytical benchmark problems

Four analytical benchmark problems (SK1, SSFYY2, MLF1, Far1) are selected from Huband et al. (2006). These problems have dimensionality ranging from one totwoandcontaintwoobjectivefunctions.T osimula teaMFen vir o nmen t,thefirstob jectiv efunctio nis considered as the high-fidel ty objective (f H ), whereas the second objective is considered as the error function (ε). Therefore, the low-fidel ty evaluation is computed as f L = f Hε. The analytical benchmark problems are shown in Figure 8 and details are provided in Table 1.D e s p i t et h e i rl o wd i m e n s i o n a l i t y ,t h e s e benchmark problems are challenging for the construction of both single-and multi-fidelity metamodels, since they are generally multimodal with several minima. Furthermore, also the error functions are multimodal. Moreover, high-and low-fidelity functions are define to show quite different behaviours, as could (moderately) happen when different physical models (such as potential fl w and RANS) are used to solve the same fluid dynamic problem.

(x) = x 4 + 3x 3 -10x 2 -10x -10 [-6, 4.5] ε(x) =-0.5x 4 -2x 3 -10x 2 + 10x -5 SSFYY2 f H (x) = 10 + x 2 -10 cos(xπ/2) [-16, 8] ε(x) = (x -4) 2 MLF1 f H (x) = (1 + x/20) sin(x) [0, 20] ε(x) = (1 + x/20) cos(x) Far1 f H (x) =-2e (15(-(x1-0.1) 2 -x 2 2 )) -e (20(-(x1-0.6) 2 -(x2-0.6) 2 )) +e (20(-(x1+0.6) 2 -(x2-0.6) 2 )) + e (20(-(x1-0.6) 2 -(x2+0.6) 2 )) +e (20(-(x1+0.6) 2 -(x2+0.6) 2 )) [-0.4, 1.2] ε(x) = 2e (20(-(x1+0.1) 2 -(x2-0.3) 2 )) + e (20(-(x1-0.4) 2 -(x2-0.6) 2 )) -e (20(-(x1+0.5) 2 -(x2-0.7) 2 )) -e (20(-(x1-0.5) 2 -(x2+0.7) 2 )) +e (20(-(x1+0.4) 2 -(x2+0.8) 2 ))

CFD-based design optimisation problems

Two CFD-based design optimisation problems are considered, namely the drag coeffici t reduction of a NACA hydrofoil and the total resistance reduction of the DTMB 5415 model (a destroyer-type vessel). Details of the optimisation problems are provided in the following subsections, whereas details on the CFD solver and simulations are provided in Section 6.

NACA hydrofoil

The following minimisation problem is solved minimise C D (x)

subject to C L (x) = 0.6 and to l ≤ x ≤ u (12)
where x is the design variable vector, C D and C L are respectively the drag and lift coeffici t of a four-digit NACA hydrofoil, and l and u are the lower and upper bound of x.Theequalityconstrain tontheliftcoefficient is necessary for the optimisation of lifting hydrofoils, in order to compare diff rent geometries at the same lift force, rather than the same angle of attack. The goal of the foil is to create a specifie lift (typically equal to the weight of the object), whereas the drag depends strongly on the lift and the foil geometry. Section 6.3.1 explains how to handle the equality constraint on the lift coeffici t. Thehydrofoilshape(seeFigure9)isdefinedbythe general equation for four-digit NACA foils. The upper (y u )andlower(y l ) hydrofoil surfaces are computed as

ξ u = ξ -y t sin θ ξ l = ξ + y t sin θ y u = y c + y t cos θ y l = y c -y t cos θ with y c = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ m p 2 2p ξ c - ξ c 2 ,0≤ ξ<pc m (1 -p) 2 × (1 -2p) + 2p ξ c - ξ c 2 , pc ≤ ξ ≤ c ( 13
)
where ξ is the position along the chord, c the chord length, y c the mean camber line, p the location of the maximum camber, m the maximum camber value, t the maximum thickness, and y t the half thickness distribution given by

y t = 5t 0.2969 ξ -0.1260ξ -0.3516ξ 2 + 0.2843ξ 3 -0.1015ξ 4 (14)
In this work, the design variables vector is defined as

DTMB 5415 model

The DTMB 5415 model (see Figure 10)i sa no p e nto-public early concept of a USS Arleigh Burkeclass destroyer, widely used for towing tank experim e n t s( I r v i n e ,L o n g o ,a n dS t e r n2008), CFD studies (Stern et al. 2001) and hull-form optimisation (Serani et al. 2016a).

The following minimisation problem is solved minimise R T (x)

subject to l ≤ x ≤ u (15)
where R T is the calm water total resistance for freesurface flow at Fr = 0.30 and Re = 1.18 • 10 7 based on the length between perpendiculars, with the ship fi ed in its design position (even keel).

The modifie hull designs (g) are produced by the linear superposition of N orthonormal basis functions (ψ) on the original geometry (g 0 ):

g(ξ , x) = g 0 (ξ ) + δ(ξ , x),w i t h δ(ξ , x) = N k=1 x k ψ k (ξ ), (16)
where ξ are Cartesian coordinates, whereas {x k } N k=1 and {ψ k } N k=1 a r et h ed e s i g nv a r i a b l e sa n dt h eb a s i s functions, provided by a design-space augmented dimensionality reduction (ADR) procedure described in [START_REF] Serani | Shape Optimization under Stochastic Conditions by Design-Space Augmented Dimensionality Reduction[END_REF]. Herein, for the sake of simplicity and demonstration, N = 2 design variables are used with x k ∈ [-1.25, 1.25].

CFD solver and problem setups

CFD simulations are performed with the unstructured finite-volume Navier-Stokes solver ISIS-CFD developed at ECN -CNRS [START_REF] Queutey | A nI n t e rf a c e Capturing Method for Free-Surface Hydrodynamic Flows[END_REF], available in the FINE TM /Marine computing suite from NUMECA Int. Computational grids are created through adaptive grid refineme t and mesh deformation to take into account the need for high and low fidel ty, as well as the different geometries needed for shape optimisation. The following subsection describes the treatment of computational grids for ISIS-CFD and some methods necessary for both the CFD-based design optimisation problems: (a) hydrofoil dynamic positioning; (b) free-surface meshing improvement and (c) limiting refineme t behind the stern for the vessel.

Adaptive grid refinement

The adaptive grid refineme t method adjusts the computational grid locally, during the computation, by dividing the cells of an original coarse grid, to improve the precision [START_REF] Wackers | Combined Refineme t Criteria for Anisotropic Grid Refine ment in Free-surface Flow Simulation[END_REF](Wackers et al. , 2017)). This method performs isotropic and anisotropic refineme t of unstructured hexahedral meshes, by dividing the cells of the original grid into fine cells; this division can be repeated several times until the desired cell sizes are obtained. The decision where to refin comes from a refineme t criterion, a tensor field C(x, y, z) computed from the flow. The tensor is based on the water surface position and on second derivatives of pressure and velocity. The mesh is refin d until the dimensions d i,j (j = 1, 2, 3) of each hexahedral cell i satisfy

C i d i,j =T r (17)
As shown by Wackers et al. (2017), the refineme t criterion based on the second derivatives of the flow is not very sensitive to grid refineme t, so the cell sizes everywhere are proportional to the constant threshold T r .T h ec o s to fth em e s ha d a p ta ti o ni sm od e r a t e ;th e added computation time for refinement is less than 5% for a typical simulation.

For multi-fidelity optimisation, the interest of this procedureisthathigh-andlow-fidelityresultscanbe obtained by running the same simulation with two diff rent thresholds T r (see Figure 11).

Computational grid deformation

The unstructured hexahedral original grids for the adaptive refineme t are generated using HEXPRESS from NUMECA Int. Since unstructured grids may be quite different for geometries that are nearly identical, numerical noise may appear in the simulation results: similar geometries may have different fl w solutions. To prevent this, the simulations of all the candidate geometries are performed with the same original grid (see Figure 12 to the parent one is propagated through these layers. Two smoothing mechanisms are applied: (i) the displacements are diffused over the faces of a single layer, so the grid deformation becomes more uniform; (ii) the displacements are multiplied with a weighting factor which goes from 1 on the geometry to 0 on the outer boundaries, so that the latter are not deformed.

An advantage of the combined deformation and adaptive refineme t is that the grid deformation (see Figure 12(b)) can be performed on the coarse initial grid, instead of a fine grid where small errors in the placement of the nodes can lead to inverted cells. T h efi n a lg r i di st h e nc r e a t e du s i n ga d a p t i v er e fi n ement (see Figure 12(c)). This procedure is more robust than attempting to deform fine grids. Furthermore, the adaptive refineme t ensures that the cells are placed effici tly for capturing the fl w around the object and in the wake, independently of the amount of mesh deformation.

NACA hydrofoil

The initial computational grid for both high-and lowfidel ty has 2654 cells, the refineme t threshold value T r is set equal to 0.1 and 0.4 for high-and lowfidel ty, respectively. This results in a cell size ratio of 4:1 between the refine fin and coarse grids. The actual computational grids have 11 k and 3.6 k cells, respectively (see Figure 13(a,b)).

T h ed o m a i nr u n sf r o m1 1 c i nf r o n to ft h el e a ding edge to 16c behind the hydrofoil and from -10c to 10c vertically. Dirichlet conditions on the velocity are imposed, except on the outfl w side which has an imposed pressure condition. The hydrofoil surface is treated with a wall law, which is the engineering stateof-the art for this type of simulation, with y + = 60 for the fi st layer. Turbulence is modelled with the standard kω SST model. The fl w is computed by time-integrating the time-dependent fl w equations towards a steady state; 4000 time steps are run for each A budget of 150 simulations is provided for the adaptive sampling methods, considering both HF and LF simulations. The initial training set for the problem isasetof2N + 1pointsincludingthedomaincentre and min/max coordinates for each variable.

Dynamic positioning

To maintain a constant lift (see Equation ( 12)), the angle of incidence α for the hydrofoil is adjusted dynamically during the simulations. At regular intervals in a simulation, the difference between the targ e ta n dt h ea c t u a ll i f ti se v a l u a t e d .T h i sd i ff e r e n c ei s divided by the theoretical lift slope of 2D foils ( C L = 2π α), to give a change in angle of attack α.T h i s change is applied over a few time steps, using grid deformation. Then, the fl w is allowed to settle and another α is computed. The procedure converged to within 0.1% of the target lift in all cases. Contrary to the deformations to fit different geometries, the computational grid is deformed with an analytical weighting technique [START_REF] Leroyer | N u m e r i c a lM e t h o d s for RANSE Simulations of a Self-propelled Fish-like Body[END_REF]to accommodate the rotation.

DTMB 5415 model

The initial computational grid for both high-and lowfidel ty has 130k cells, which allows to capture the base flow features. As for the NACA hydrofoil problem, a 4:1 cell size ratio between the HF (T r = 0.0145) and LF (T r = 0.0036) computational grid is selected. The actual computational grids have approximately 4.3 M and 260 k cells, respectively. This ratio appears ideal for the present problem: for higher ratios, either the HFcomputationsaretoocostlyortheLFgridsareso coarse that they cannot capture the ship geometry.

Simulations O n c e again, the hull is treated with a wall law, y + = 60 for the first layer, and turbulence is modelled with kω SST. The free-surface deformation for wavemaking is treated with a surface-capturing approach that resolves bothwaterandairflow.Asusualforfree-surfaceflows, the flow is solved through time integration; 2000 time steps are used. The forces are averaged over the fin l 250 time steps.

On a 20-core workstation the HF and LF simulations take about 24 and 1.5 hours, respectively, corresponding to a computational cost ratio β = 0.063.

The initial training set for MF-DRBF is a set of 2N + 1 points including the domain centre and with each single design variables at either +1or-1. Adaptive sampling is then run until 80 computations, HF and LF combined, are reached.

Improving the free-surface meshing

To provide adequate surface capturing even on a deformed computational grid, the fine grid around the free surface is created entirely by adaptive refinement (see Figure 12(c)); the original grid has no specifi 2012)weighting(c)Deformedwith power 3.0 weighting refineme t around the free surface at all. This implies that the cells of the original mesh have to be divided many times by the grid adaptation. However, thanks to the anisotropic nature of the refinement algorithm [START_REF] Wackers | Combined Refineme t Criteria for Anisotropic Grid Refine ment in Free-surface Flow Simulation[END_REF], only refineme t normal to the free surface is requested. On undeformed grids (see Figure 14(a)) this works well: where the surface is at rest,thecomputationalgridisrefinedonlyinvertical direction, resulting of good quality. When the computational grid is deformed (see Figure 14(b)) the cells of the original grid are stretched and rotated. As a consequence, the cells are no longer aligned with the free surface,sotheyarerefinedinmorethanonedirection, which leads to unnecessary refineme t and mediocre grid quality. This problem is alleviated by changing the weighting factor. Instead of the weighting defined by [START_REF] Durand | Light and Flexible Fluid/Structure Interaction, Application to Sailing Boats[END_REF], which is based on the linear and s q u a r e dd i s t a n c et ot h eb o d y ,aw e i g h t i n gb a s e do n the distance to the power 3.0 is chosen. Since this weighting law goes to zero rapidly when the distance tothebodyincreases,thegridfarfromthebodyisless d e f o r m e d .A sar e s u l t ,t h er e fi n e dm e s hi si m p r o v e d (see Figure 14(c)).

Limiting refinement behind the stern

When simulating ships with adaptive refinement, unrestricted application of the velocity/pressure refineme t criterion leads to the accurate resolution of thenearandthefarwake(W ackers,Guilmineau,and Visonneau 2017). However, the far wake may not be required for drag evaluation. Therefore, tests were run where refineme t is forbidden from a certain distance b e h i n dt h es t e r n ;t h u s ,t h ef a r -w a k er e g i o nr e t a i n sa coarse mesh. Figure 15 shows the grid at the undisturbed water level for three cases: the refinement is suppressed from 0.05 L, 0.25 L, and 1.75 L behind the stern respectively. To study the effect on the computed drag, such series were computed for two geometries created by varying the optimisation design variable x 2 , as well as the undisturbed case (see Figure 16). The results for the middle limit (0.25 L) are quite close to t h el a r g e s to n e( 1 . 7 5L ) ,w i t ht h eb i g g e s td i s c r e p a n cy occurring for x 2 =-1.25. The shortest limit (0.05 L) however, changes the drag signifi antly when the mesh is deformed to x 2 = 1, probably because the stern wave p a t t e r ni sm o d i fi e d ;t h ec a s ex 2 = 1h a st h es t e e p e s t stern wave of the three. As a compromise between the computational cost and the accuracy of the total resistance estimation in all deformation cases, the resolution of the near wake (0.25 L) is chosen.

Numerical results

T h ef o l l o w i n gs u b s e c t i o n sp r e s e n tt h er e s u l t so ft h e analytical benchmarks and the CFD-based design optimisation problems. For practical purposes, solving the Equations ( 6), ( 8)-( 10), a distance-based penalisation is applied if the new candidate sample has a distance lower than 0.001 (considering a normalised design variable range) from other points of the training set. The penalisation aims at preventing sample clusterisation and ill-conditioning problems. The above minimum distance between training points is also used as stopping criterion. The number of samples (N P )c o m p u t e da to n c eb yt h eM C A Ss a m p l i n g method is set to N P = 4. Furthermore, a threshold U ⋆ f = 0.5%R is used for the MCAS (see Equation ( 10)). The computational cost is normalised with the HF cost: one HF evaluation costs 1, whereas the LF cost is equal to β, for a total cost equal to N H + βN L ,where N H and N L are the numbers of high-and low-fidelity computations, respectively.

Analytical benchmark

Analytical benchmark results are assessed with two metrics: (1) the objective function value for the current optimum (high-fidel ty evaluated) and (2) the normalised root mean square error (NMRSE) of the multi-fidel ty prediction versus the analytical function (normalisation is based on the initial HF training set range). The NRMSE is computed on a set of 50 N evenly-spaced points. The available computational budget is set equal to a maximum cost of 50N. Due to the analytical nature of the test functions, the computational costs are not significant; a computational cost ratio β = 0.1 between low-and high-fidel ty is artifici lly assumed.

Figure 17 shows the convergence of the metamodelbased optimisation towards the objective minimum and the corresponding NMRSE conditional to the adaptive sampling methods for problems SK1, SSFYY2, MLF1, and Far1. All methods achieve the minimum, excluding the MFEI sampling for the MLF1 problem (note that MFEI faces premature convergence since it has activated the clustering stopping criterion). It is worth noting that the ACAS shows better convergence to the minimum than other methods. Differently, considering the global approximation, the fastest and most signific nt NMRSE reduction is achieved by t h eM U A S .T h eM C A Sh a ss i m i l a rp e r f o r m a n c et o MUAS.

S i n c eA C A Ss h o w st h eb e s tp e r f o r m a n c ef r o m the optimisation viewpoint, a parametric analysis for ACAS only is conducted considering β = {0.05, 0.1, 0.2, 0.4, 0.8}. Moreover, a single-fidel ty DRBF t r a i n e do n l yb yH Fe v a l u a t i o n si su s e dt oa s s e s sa n d compare the MF-DRBF efficie y. The corresponding results for the analytical benchmarks are shown in Figure 18.T h eM Fa p p r o a c ha c h i e v e sf a s t e rc o nvergence towards the minimum and lower NRMSE than using single HF-based metamodel for problems SK1, SSFYY2, and MLF1, regardless of β.Forproblem MLF1, the HF-based metamodel does not achieve the objective minimum. The Far1 problem shows comparable convergence towards the objective minimum of MF and HF-based DRBF, whereas in terms of NRMSE the HF-based DRBF achieved the lowest value.

NACA hydrofoil

The drag optimisation of the NACA four-digit hydrofoil Equation ( 12) is a two-parameter problem like Far1, but it only has a single, global minimum. Theref o r e ,ac o n v e r g e n c el ik ef o rth eS S F YY 2p r o b l e mca n be expected. A particular diffic ty of this problem is that, while the minimum lies on the domain boundary, it is located close to high drag values which may preventthesamplingofpointscloseby .Theadaptive MF-DRBF-based optimisation results are compared to an optimal benchmark solution for Equation ( 12) identifie by an earlier single-fidel ty metamodelbased optimisation trained by 150 HF simulations (Ploé et al. 2017).

Figure 19 shows the convergences of the maximum prediction uncertainty, the drag coeffici t minimum, and the corresponding design variables. MUAS and MCAS have achieved lower values of the maximum prediction uncertainty than MFEI and ACAS. On the contrary, only MFEI and ACAS have reached a minimum close to the benchmark. Overall, MFEI has provided the fastest convergence towards the minimum.

Figure 20 shows the MF-DRBF prediction and the corresponding training set at the final iteration of the four sampling strategies. The MUAS and MCAS strategies have provided a global exploration of the domain. On the contrary, the MFEI and ACAS methods have clustered training points in a small region, close to the global minimum. All the sampling strategies have used a similar number of HF evaluations (see Table 2), with MUAS and MCAS spreading the HF evaluations over the whole design space. Furthermore, both MUAS and MCAS have requested HF training points in the design space corners. The MFEI and ACAS strategies have focused the HF evaluations only in the minimum region.

Furthermore, Figure 19 shows that the convergence of the uncertainty prediction of MUAS is noisy. This happens since MUAS aims at the minimisation of the maximum uncertainty, therefore clusters the training points in regions where the uncertainty is high (see Figure 20). However, the MF-DRBF uncertainty is affected by a small-scale noise found in the numerical simulations, in regions with clustered training points. (This noise is due to different geometries producing different grids, so two almost similar geometries may have different numerical errors in the solution. Another source is the small variations in lift coming from the dynamic positioning.) Therefore, MUAS has added training points in regions where the uncertainty is due to the numerical noise and not to the global objective function shape. The effects of the numerical noise are evident in Figures 20(a). MUAS has clustered samples in the neighbourhood of x ={4.0, 4.0}× 10 -2 (see Figure 20 Finally, Figure 21 shows the velocity contours of the optimal NACA hydrofoil shape identified by the adaptive sampling strategies. The differences of the hydrofoil shapes are not evident, although the velocity contours are different near the leading edge.

The results of the adaptive sampling strategies are summarised in Table 2, showing that the MFEI method is the most effective adaptive sampling technique for the present problem, followed closely by ACAS.

DTMB 5415 model

Figure 22 shows the MF-DRBF prediction and the corresponding training set at the final iteration of the four sampling strategies. The figu e shows that, even if the design optimisation is a challenging problem from the CFD perspective, the objective function shape is actually simpler than for the previous problems, since it is unimodal and has a single minimum. Thus, all sampling methods produce similar MF-DRBF metamodelsandtheoptimafoundareclose.Nevertheless, the samples spreading are completely diff rent. The MUAS has concentrated samples on the corners of the design space, following almost a symmetric pattern. Apart from the 5 initial training points, HF samples have been added in the bottom corners. The MFEI s h o w ss o m ee r r a t i cs a m p l i n g ;a sn o t e db e f o r e( s e e Figure 17) the expected improvement goes to zero at a given moment, meaning that the samples placement is stopped before exhausting the budget. However, with three HF samples lying close to the optimum (one hidden by LF samples in Figure 22) the error metamodel is locally reliable. The ACAS has concentrated all the samples around the optimum. Only one HF sample has been added, close to the optimum location (hidden behind the LF samples in Figure 22). As for the NACA problem, MCAS takes a middle road, adding samples around the optimum but also in the corners. R e g a r d i n gH Fs a m p l e s ,a l la d d i t i o n a ls a m p l e sa r ei n the design space corners. Finally, considering the MF-DRBF uncertainty, only MUAS and ACAS achieve a uniform low uncertainty around the optimum. Table 3 summarises the locations of the optimum a n dt h eo p t i m i se dr e s i s t a n c ef o rt h es a m p l i n gs t r a t egies. Apart from MFEI, the optimum locations lie close together. ACAS has achieved the lowest predicted total resistance. The actual resistance improvement is about 4% for all methods: the MUAS is best overall, whereas the MFEI has the lowest validation error (about 0.1%) with respect to the original DTMB 5415 geometry. Thisca nbeexp la inedb ytheerr o rmeta model.Since β = 0.0625 is low for this problem, few HF samples are added in the training set, meaning that the error metamodel is not define very accurately close to the optimum, resulting in a discrepancy between the predicted and the actual optima. Since MFEI has three HF points around the optimum, the error metamodel is well defin d in this region and the MF-DRBF prediction with MFEI sampling is the most reliable for this problem.

Figure 23 shows the initial and the MFEI optimal geometry (this is similar to the other optimised geometries). The optimum, compared to the original geometry, has a more slender aftship which reduces the width of the stern wave, and a slight bulge aft of the bow. This bulge creates a second bow wave out of phase with the firs one. By superposition, these two waves cancel so the total bow wave is reduced. Since this interference depends on the wave lengths, which are determined by the ship's velocity, the total resistance reduction is expected to be effective only around the target velocity.

Conclusions and future work

Four adaptive sampling methods for MF-DRBF metamodels have been presented. These are assessed by four analytical benchmarks and applied to two CFD-based design optimisation problems, pertaining to the resistance reduction of a NACA hydrofoil and a destroyer-type vessel (DTMB 5415). The MF approximation is obtained as the sum of a low-fidel ty-trained metamodel and the metamodel of the error (difference) between high-and low-fidel ty evaluations. The stochastic RBF provides the metamodel prediction and the associated uncertainty. The prediction uncertainty of both the low-fidelity and the error metamodel is used for the adaptive refineme t of the low-fidel ty and the error training sets. The ratio of the computational cost of high-and low-fidel ty evaluations affects the choice of the fidelity to sample.

Thecriteriafortherefinementofthemulti-fidelity training set are: (i) the minimisation of the maximum uncertainty of the MF metamodel prediction (MUAS), (ii) the maximisation of the multi-fidel ty expected improvement (MFEI), (iii) the minimisation ofanaggregatedmeritfactorofthepredictionuncertainty and the predicted objective function (ACAS), and (iv) a multi-objective optimisation aiming at maximising the prediction uncertainty and minimising the predicted objective function (MCAS). The MF-DRBF performance has been assessed in terms of g en cet o wa r d sth eminim um.I naddi tio n,MF -D RB F maximum uncertainty and NRMSE are used for further comparisons.

Theanalyticalbenchmarksresultshaveshownthat the ACAS (closely followed by the MFEI sampling) is the most effective method for the MF-DRBF, from the optimisation point of view. For this reason, a further analysis for the ACAS conditional to the computational-cost ratio (β) has been conducted. The ACAS strategy combined with the MF formulation hasshownbetterperformancethanusingasingleHF metamodel, regardless of β.

The CFD-based design optimisation problems are challenging for the MF metamodel. The existence of numerical noise affects the RBF interpolation, resulting in large uncertainty of the MF prediction in noisy regions of the domain. Therefore, sampling methods that directly take into account the prediction uncertainty are 'trapped' in such regions. For both CFDbased problems ACAS and MFEI are found the most effective sampling strategies.

In order to overcome the difficulties arising from the presence of numerical noise (associated to CFD simulation residuals), future work will consider the use of regularisation approaches, approximation/regression methods (e.g. least square fit) as opposed to exact interpolation, as well as extensions to the use of Gaussian process. Preliminary studies by the authors can be found in [START_REF] Wackers | A d a p t i v eM u l t i fidelity Shape Optimization Based on Noisy CFD Data[END_REF]a n dA n t o gnoli et al. (2019). Moreover, the possibility to generalise the MF concept to more than two fidelities (multi-level) is currently being investigated [START_REF] Serani | A n Adaptive N-Fidelity Metamodel for Design and Operational-Uncertainty Space Exploration of Complex Industrail Problems[END_REF] and will be further developed in future studies.

Disclosure statement

No potential conflict of interest was reported by the authors.

Funding

CNR-INM is grateful to Dr. Salahuddin Ahmed and Dr. Woei-M i nL i no ft h eO ffi c eo fN a v a lR e s e a r c h ,f o rt h e i rs u p p o r t through NICOP grant N62909-18-1-2033.

  • 10 6 and (b) the total resistance of the DTMB 5415, an early and open-to-public version of the DDG-51 (a USS Arleigh Burke-class destroyer-type vessel), in calm water at Froude number equal to 0.3. CFD simulations are based on the Reynolds-averaged Navier-Stokes equation solver ISIS-CFD, developed at Ecole Centrale de Nantes/CNRS and integrated in the FINE/Marine simulation suite from NUMECA Int.

Figure 1 .

 1 Figure 1. Multi-fidelity metamodel concept and notation.

Figure 2 .

 2 Figure 2. Updating scheme for the adaptive multi-fidelity metamodel.

Figure 3 .

 3 Figure 3. MUAS sampling strategy.

Figure 4 .

 4 Figure 4. MFEI sampling strategy.

Figure 5 .

 5 Figure 5. ACAS sampling strategy.

Figure 6 .

 6 Figure6. MCAS down-sampling strategy: green zone is defined to avoid over-fitting in the neighbourhood of the minimum.

Figure 7 .

 7 Figure 7. MCAS sampling strategy.

Figure 8 .

 8 Figure 8. Analytical benchmark problems.

Figure 9 .

 9 Figure 9. NACA 4-digit hydrofoil.

x

  ={t, m} with t ∈ [0.030, 0.120] and m ∈ [0.025, 0.070]. The maximum camber position is fi ed at p = 0.4. The simulation conditions are: velocity U = 10 m/s, chord c = 1m,fluiddensityρ = 1, 026 kg/m 3 , and Reynolds number Re = 8.41 • 10 6 based on the chord length.

Figure 10 .

 10 Figure 10. A geosym replica of the DTMB 5415 (CNR-INSEAN model 2340).

  (a)), which is deformed to fit each geometry (see Figure12(b)). The deformation algorithm[START_REF] Durand | Light and Flexible Fluid/Structure Interaction, Application to Sailing Boats[END_REF])d i v i d e st h eg r i di n layers of cells around the geometry. For each geometry face, the displacement of the faces with respect

Figure 11 .Figure 12 .

 1112 Figure 11. Coarse and fine computational grid example for the DTMB 5415 problem. The figure shows the mesh in the plane z = 0, coloured by the water volume fraction. (a) Coarse grid, T r = 0.0145 (b) Fine grid, T r = 0.0072

Figure 13 .

 13 Figure 13. NACA hydrofoil computational grids. (a) Coarse grid, 3.6 k cells (b) Fine grid, 11 k cells

  are performed on the half geometry, in model scale with L = 5.72 m. The domain runs from 1.5 L (with L thevessellength)infrontofthebowto3L behindthestern,u pto2Lla terally ,andfr o m-1.5 L to 0.5 L vertically. Dirichlet conditions on the velocity a r ei m po sedo nth ei n fl o wa n ds i d ef a c e s ,p r e s s u r ei s i m p o s e do nt h et o p ,b o t t o m ,a n do u t fl o ws i d e .

Figure 14 .

 14 Figure 14. Exampleofundisturbedfree-surfacegrid(showninx-cut planes) for DTMB 5415 problem, depending on the computational grid deformation algorithm (Red is water, blue is air). (a) Undeformed (b) Deformed with Durand (2012)weighting(c)Deformedwith power 3.0 weighting

Figure 15 .Figure 16 .

 1516 Figure 15. Computational grid refinement behind the stern for DTMB 5415 problem, conditional to the refinement limit. (a) -0.05 L (b) -0.25 L (c) -1.75 L

Figure 17 .

 17 Figure 17. Analytical benchmark problems: adaptive sampling performance comparison for β = 0.1.

Figure 18 .

 18 Figure 18. Analytical benchmark problems: ACAS performance conditional to β (only HF is used when β = 1).

Figure 19 .

 19 Figure 19. NACA hydrofoil problem: convergence of adaptive MF-DRBF maximum uncertainty, drag coefficient minimum, and corresponding coordinates x 1 and x 2 .

Figure 20 .

 20 Figure 20. NACA hydrofoil problem: training sets and drag coefficient MF-DRBF prediction at the final iteration of the adaptive sampling procedures.

  (a)) mainly LF training points; this region has the highest numerical uncertainty in the

Figure 21 .

 21 Figure 21. Hydrofoil optimisation: velocity contour for the optimal configuration identified with the sampling methods. (a) MUAS (b) MFEI (c) ACAS (d) MCAS

Figure 23 .

 23 Figure 23. Initial (top) and optimised (bottom) geometry and wavepattern (z-coordinate of the free-surface position).

Table 1 .

 1 Analytical benchmark problems.

	Name	Problem	Domain
	SK1	f H

Table 2 .

 2 Summary of the adaptive sampling methods performance on the NACA hydrofoil problem. The table provides both the metamodel-predicted ĈD and an actual simulated C D in the computed optimum.

				m ¸inimum position	m ¸inimum value
		N H [-]	N L [-]	x 1 [-]	x 2 [-]	ĈD [-]	C D [-]
	MUAS	14	136	3.8878E-2	3.0000E-2	7.1759E-3	7.2582E-3
	MFEI	15	135	4.1691E-2	3.0000E-2	7.1545E-3	7.2371E-3
	ACAS	12	138	4.1290E-2	3.0009E-2	7.1754E-3	7.2403E-3
	MCAS	11	137	3.9810E-2	3.0000E-2	7.2816E-3	7.2606E-3
	Benchmark	150		4.0993E-2	3.0027E-2		7.2340E-3

Table 3 .

 3 Summary of the performance of the adaptive sampling strategies on the DTMB 5415 problem.

				m ¸inimum position	m ¸inimum value
		N H [-]	N L [-]	x 1 [-]	x 2 [-]	RT [N]	R T [N]
	MUAS	8	72	0.4425E-2	-7.8258E-1	48.4044	48.5456
	MFEI	7	73	0.3725E-2	-9.9753E-1	48.5508	48.6088
	ACAS	6	74	-1.0500E-2	-7.4843E-1	48.3062	48.6236
	MCAS	8	70	0.3750E-2	-7.3195E-1	48.4409	48.7344
	Original			0.0	0.0		50.6465