Luc Laurent

Pierre-Alain Boucard

Bruno Soulier

Generation of a cokriging metamodel using a multiparametric strategy

Keywords: Cokriging metamodel, multiparametric strategy, LATIN method, multilevel optimization

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Optimization techniques have been in use in industry for many years. Structural elements, geometries, materials... are often designed and selected by means of specific industrial optimization software. The main difficulty in carrying out an optimization is that the precise localization of the optimum of an objective function often requires a large number of calculations. In the case of the design of structural assemblies, each calculation involves the resolution of very complex nonlinear problems due to the existence of contacts and friction between the parts. This makes the computation cost a major scientific stumbling block.

A technique widely used in order to make such optimizations affordable is multilevel model optimization [START_REF] Robinson | A case for multi-level optimisation in aeronautical design[END_REF]. Our approach belongs to that category in the sense that, contrary to the surrogate-based approach [START_REF] Booker | A rigorous framework for optimization of expensive functions by surrogates[END_REF][START_REF] Queipo | Surrogate-based analysis and optimization[END_REF], we use a metamodel only to accelerate the optimization process. The main steps are: first, the construction of a metamodel from the responses of the mechanical model; then, the localization of an approximate optimum using the metamodel; and, finally, the precise determination of the optimum using the mechanical model and the approximate optimum. This strategy leads to a reduction in computation cost and, therefore, in the time it takes to complete the optimization.

This paper focuses on the first step, i.e. the construction of the metamodel using:

a cokriging metamodel; a dedicated computational strategy.

We will show that when gradients are available a cokriging metamodel leads to a better approximation than a kriging metamodel with the same number of sample points. Moreover, our computational strategy is capable of reusing previous calculations in order to solve new problems. When the parameters vary, this property enables one to reduce the cost of a new calculation, particularly when it comes to evaluating gradients using a finite difference method.

The first chapter of the paper is a review of multilevel optimization techniques in the context of this study. More precisely, it focuses on multilevel model optimization intended for parameter optimization. The second chapter introduces the multiparametric strategy and discusses its performance in detail. The third chapter proposes and develops a gradient-based formulation of the cokriging metamodel.

In the last chapter, numerical and analytical examples are used to show that the cokriging metamodel provides a better approximation than classical (e.g. kriging) metamodels. For mechanical examples, the coupling of a cokriging metamodel with our computational strategy leads to a significant reduction in the computation costs associated with the generation of the metamodel. The quality of the cokriging metamodel is studied based on several classical criteria and compared with that of the kriging metamodel.

The optimization process

Classical direct optimization, in which the optimizer is linked directly to the mechanical solver, requires huge computer resources. The cost increases not only with the complexity of the mechanical problem (i.e. the nonlinearities included in the model), but also with the number of degrees of freedom in the mechanical system and the number of design variables in the optimization problem. In order to reduce the cost of optimization in structural design, many works concentrate on three main aspects: the development of dedicated numerical methods to address increasingly complex mechanical problems, the improvement of the performance of the optimization algorithms and the development of strategies to coordinate the exchanges between the optimization algorithm and the mechanical solver (for example, see [START_REF] Queipo | Surrogate-based analysis and optimization[END_REF][START_REF] Simpson | Metamodels for computer-based engineering design: survey and recommendations[END_REF]). The use of a multilevel optimization strategy provides a solution to the last aspect. Multilevel strategies were developed over the last twenty years. These approaches can be divided into two main categories: multilevel parameter optimization (which includes sequential methods [START_REF] Kravanja | The multilevel minlp optimization approach to structural synthesis: the simultaneous topology, material, standard and rounded dimension optimization[END_REF][START_REF] Kravanja | Efficient multilevel minlp strategies for solving large combinatorial problems in engineering[END_REF][START_REF] Liu | Global-local structural optimization using response surfaces of local optimization margins[END_REF] and iterative methods [START_REF] Alexandrov | Analytical and computational aspects of collaborative optimization[END_REF][START_REF] Chen | Multidisciplinary design optimization of mechanisms[END_REF][START_REF] Conceic ¸ão António | A multilevel genetic algorithm for optimization of geometrically nonlinear stiffened composite structures[END_REF][START_REF] Bendsøe | Optimization of structural topology, shape, and material[END_REF][START_REF] Le Riche | Design of dimensionally stable composites by evolutionary optimization[END_REF][START_REF] Theocaris | Multilevel optimal design of composite structures including materials with negative poisson's ratio[END_REF]) and multilevel model optimization (which includes hierarchical multilevel optimization [START_REF] Robinson | A case for multi-level optimisation in aeronautical design[END_REF][START_REF] Keane | Aircraft wing design using ga-based multi-level strategies[END_REF] and imbricated multilevel optimization [START_REF] Engels | Implementation of a multilevel optimisation methodology within the e-design of a blended wing body* 1[END_REF]). Similar strategies have also been used in the context of multidisciplinary optimization where each level of the optimization process concerns a specific discipline [START_REF] Chen | Multidisciplinary design optimization of mechanisms[END_REF][START_REF] Chattopadhyay | Multilevel decomposition procedure for efficient design optimization of helicopter rotor blades[END_REF]. Parallel multilevel model optimization [START_REF] El-Sayed | Optimum structural design with parallel finite element analysis[END_REF][START_REF] Umesha | Optimal design of truss structures using parallel computing[END_REF] can also be used in order to reduce the optimization cost and take maximum advantage of parallel computing architectures. This type of multilevel strategy calls for specific optimization algorithms, such as genetic algorithms [START_REF] Dunham | Design by natural selection[END_REF][START_REF] El-Beltagy | A comparison of various optimization algorithms on a multilevel problem[END_REF][START_REF] Elby | Evaluation of injection island ga performance on flywheel design optimisation[END_REF].

The multilevel modeling strategy we propose is based on multilevel model optimization. We focus on the two-level strategy illustrated in Figure 1. On the first level, which consists in a metamodel defined using a limited amount of data, the zones where the optimum is to be sought are determined using, for example, a genetic algorithm [START_REF] Holland | Adaptation in natural and artificial systems[END_REF]; then, this information is transferred to the second level, where a precise search for the minimum can be carried out using a gradientbased algorithm. Thus, the second level is an optimization process based on the full mechanical model.

Simulator

Metamodel Optimizer

Simulator Optimizer

Optimization using the metamodel Direct optimization

Scope of this article

Fig. 1: Multilevel optimization

The scope of this paper is limited to the study of the cost associated with the generation of the metamodel.

Many types of metamodels can be used to find an approximate solution: they can be based on polynomial regression [START_REF] Box | On the experimental attainment of optimum conditions[END_REF][START_REF] Simpson | Comparison of response surface and kriging models for multidisciplinary design optimization[END_REF], on neural networks [START_REF] Haykin | Neural networks: a comprehensive foundation[END_REF][START_REF] Mcculloch | A logical calculus of the ideas immanent in nervous activity[END_REF], on radial basis functions [START_REF] Hardy | Multiquadric equations of topography and other irregular surfaces[END_REF], on proper orthogonal decomposition [START_REF] Chatterjee | An introduction to the proper orthogonal decomposition[END_REF], on cumulative interpolation [START_REF] Soulier | Crashworthiness optimization using a surrogate approach by stochastic response surface[END_REF], etc... We chose to use a particular class of metamodels called kriging approximations [START_REF] Cressie | The origins of kriging[END_REF] and, more precisely, a cokriging metamodel [START_REF] Wackernagel | Multivariate geostatistics: an introduction with applications[END_REF] using derivatives [START_REF] Chung | Using gradients to construct cokriging approximation models for high-dimensional design optimization problems[END_REF]. These approximations are presented in Section 4. The interested readers could also refer to works on metamodeling strategies in structural optimization and multidisciplinary design optimization [START_REF] Barthelemy | Approximation concepts for optimum structural design -a review[END_REF][START_REF] Sobieszczanski-Sobieski | Multidisciplinary aerospace design optimization: survey of recent developments[END_REF].

The multiparametric strategy

This chapter introduces the concept of multiparametric strategy, which enables one to solve similar structural problems at a greatly reduced cost. This is used for the types of structural assembly problems which are described in Section 3.1. The strategy is based on the LATIN method, which is presented in Section 3.2. The multiparametric strategy itself is presented in Section 3.3 and its performance is discussed using the example of a 2D mechanical benchmark problem with two design variables in Section 3.3.2.

The problem of structural assemblies

The structures being considered in this paper are assemblies of linear elastic structures under the assumption of small perturbations. The only nonlinearities occur between parts of the assemblies and are due to contact and/or friction phenomena. In order to solve these problems, we use a dedicated strategy based on the LATIN algorithm introduced by P. Ladevèze [START_REF] Ladevèze | Nonlinear computational structural mechanics: new approaches and non-incremental methods of calculation[END_REF].

This strategy is based on three main points:

-The structure being studied is divided into substructures and interfaces. This is a natural approach in the context of assemblies because treating each part of an assembly as a substructure and the contact zone(s) between two parts as (an) interface(s) constitutes the simplest decomposition. What makes this decomposition unique is the use of mixed force and displacement unknowns at the interfaces; -A suitable iterative algorithm is used to solve the mechanical problem; -The operators of the method remain constant and do not depend on the loading or on the parameters of the interface (friction coefficient, gap).

The resulting approach is a mixed domain decomposition method, as opposed to the primal substructuring approach [START_REF] Mandel | Balancing domain decomposition[END_REF][START_REF] Tallec | Domain decomposition methods in computational mechanics[END_REF] or dual approach [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF].

Short summary of the LATIN method for assemblies

Our resolution strategy is based on the LATIN algorithm developed by P. Ladevèze [START_REF] Ladevèze | Nonlinear computational structural mechanics: new approaches and non-incremental methods of calculation[END_REF], which consists in solving two problems alternatively: one in the substructures and one at the interfaces. These problems are described by two groups of equations denoted A d and Γ : A d contains the linear equations related to the substructures and Γ contains the local equations (which can be nonlinear) related to the interfaces. The iterations between the two groups are carried out through the use of search directions which are parameters of the method. This resolution process, shown in Figure 2, leads to the solution defined as the intersection of spaces Γ and A d . In the context of elastic assemblies, the boundary conditions (forces F i and displacements W i ) at each interface i between two substructures are sufficient to define the solution: the internal behavior of each substructure is the solution of a classical elastic problem. Thus, an approximate solution is described entirely by s = all the interfaces The friction and contact conditions with their respective approaches and laws [START_REF] Blanzé | A modular approach to structure assembly computations: application to contact problems[END_REF] are part of equation group Γ . The linear form of A d is due to the linear behavior of the substructures.

(W i , F i ). A d Γ E - E + ŝn s s n+1 s n
The main feature of the LATIN method is that at each iteration the resolution leads to an approximate solution over the whole loading path and in all the points of the structure. Each new iteration enriches this solution until convergence.

Additional details can be found in [START_REF] Ladevèze | Nonlinear computational structural mechanics: new approaches and non-incremental methods of calculation[END_REF][START_REF] Ladevèze | Sur une famille d'algorithmes en mécanique des structures[END_REF] for the LATIN method and in [START_REF] Blanzé | A modular approach to structure assembly computations: application to contact problems[END_REF][START_REF] Champaney | Modular analysis of assemblages of three-dimensional structures with unilateral contact conditions[END_REF] for its application to assemblies.

The multiparametric strategy

The construction of the metamodel requires the evaluation of what is called in optimization an objective function or cost function at certain points of what is called the design space. In the case of a parametric study in the mechanical context, these evaluations lead to the resolution of many problems which are similar in the sense that only their parameters vary. In order to accelerate these resolutions, we use an approach called the MultiParametric Strategy (MPS) which was introduced in [START_REF] Boucard | Une application de la méthode latin au calcul multirésolution de structures non linéaires[END_REF][START_REF] Boucard | Application of the latin method to the calculation of response surfaces[END_REF], then studied and applied to various types of problems in [START_REF] Allix | A new multi-solution approach suitable for structural identification problems[END_REF][START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF][START_REF] Champaney | Adaptive multianalysis strategy for contact problems with friction[END_REF]. The following presentation introduces the principle of this strategy and discusses its performance based on an academic test example.

Principle of the multiparametric strategy

The main idea of the multiparametric strategy is very simple and consists in initializing the LATIN algorithm using the converged solution of a previous problem. As mentioned in Section 3.2, each calculated solution, associated with a set of parameters, is entirely described by the values of the variables at the boundary. Thus, the initialization process consists in reloading the boundary values associated with the converged solution of the chosen similar problem.

In order to build a metamodel, the set of the sample points (each defined in the design space by a set of parameters) can be divided into two sets: the set of the sample points associated with calculated values of the objective function, and the set of the sample points at which the values of the objective function are still unknown. The latter is not always defined clearly, especially in an optimization process in which the optimizer introduces new sample points gradually, but the first set constitutes a database which is available for subsequent calculations. Thus, if a new calculation associated with different values of the parameters is required, the algorithm can be initialized using a previous converged solution taken from the database. Each new evaluation of the objective function enriches the database.

Figure 3 illustrates the use of the multiparametric strategy to obtain an approximate solution starting from a previous calculation: in order to obtain the converged solution, s 3 can be initialized with one of the two available converged solutions s 1 or s 2 . The benefit depends on the choice of the starting point. Some indications for choosing the best initialization strategy were given in [START_REF] Soulier | A multiparametric strategy for the large-scale multilevel optimization of structural assemblies[END_REF]. In our case we will use only a "closest point" strategy: the initial solution of the new problem is chosen to be the converged solution associated with the set of design parameters which is closest to the one being considered.

Since we are assuming that the only parameters which vary are the interface parameters, space Γ alone is affected by the change. Thanks to this strategy, the solution converges in fewer iter-

A d Γ 1 Γ 2 Γ 3 s 1 s 2 s 1 s 3
Variat ion of param eters Fig. 3: The multiparametric strategy using the LATIN algorithm ations of the algorithm and, therefore, in less time.

Performance of the multiparametric strategy

Let us consider the example of a quasi-static academic problem which was presented in [START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF]. Figure 4 shows the geometry of the problem, which consists of three square parts (h = 50mm, Young's modulus E = 2 • 10 5 MPa and Poisson's coefficient ν = 0.3) in contact with friction. Each part is represented by a single substructure discretized into 20×20 bilinear quadrangles. The parametric study consists in varying the friction coefficients µ 1 and µ 2 of the two contact interfaces. The loading consists of two stages: first, a progressive vertical pressure P 1 up to a maximum of 50MPa applied at the top of substructure Ω 3 (the preloading stage), then a progressive horizontal load from 0 to 30MPa applied to substructure Ω 2 . In this test, variations of the friction coefficients µ 1 and µ 2 between 0 and 0.6 were considered, and the function studied was the reaction force on the rigid wall. Figure 5 shows the response surface of this function obtained with 18×18 values of the friction coefficients.

In order to illustrate the performance of the multiparametric strategy, the two-variable design space was sampled with a 4×4 regular grid. For each sample, the force and its gradients were calculated (using a classical finite difference method for the gradients). Table 1 summarizes the characteristics of these calculations.

The gain obtained with our method compared to a classical calculation without the MPS was estimated using the expression: Gain = CPU time without the MPS CPU time using the MPS

The simple mechanical examples used in this paper make it possible to evaluate the CPU time without the MPS, but in the general case one uses the formula:

Gain = Number of calculations×Reference CPU time CPU time using the MPS (2)
where the Reference CPU time can be the CPU time of the first calculation, the maximum CPU time necessary to obtain one evaluation, or the average CPU time for evaluations at several points of the design space without the MPS. The results presented in Table 1 show that the solver can indeed accelerate the resolution of similar problems by reusing a previously converged solution. The most remarkable gain was that obtained in the evaluation of the gradients. These were calculated using a finite difference method: each gradient required the value of the response in 3 points (a sample point plus 2 neighboring points obtained with very small variations of the parameters). The results show that for the same number of evaluations the cost was lower when both the responses and the gradients were calculated than when the responses alone were calculated: in our example, 16 evaluations of the responses alone took 47.5s, but the addition of 16 gradients (32 evaluations) took only 48.9s. Thus, thanks to the ability of the method to reduce the computation time for sample points that are close together, the strategy is particularly efficient in calculating the gradients. This feature is very interesting if one wishes to use a gradient-based optimizer to achieve a significant reduction in optimization time. In this paper, we will also undertake to use these inexpensive gradients to build a gradient-enhanced metamodel.

The cokriging metamodel

The classical approach to building a metamodel consists in using the responses of an objective function (also called a primary variable) calculated at a number of sample points. The values of the objective function at the chosen sample points can be used to build a classical metamodel. However, in many cases, it may be interesting to use additional information in the form of auxiliary variables to build a richer metamodel. Multivariate geostatistics is a field of applied mathematics which supplies methods to handle these variables jointly. The cokriging metamodel we use [START_REF] Wackernagel | Multivariate geostatistics: an introduction with applications[END_REF] is an example of such methods in which the primary variable is the objective function and the auxiliary variables are its gradients.

The kriging technique was first introduced by D.G. Krige, a mining engineer [START_REF] Krige | A statistical approach to some mine valuation and allied problems on the Witwatersrand[END_REF]. Subsequently, many mathematical enhancements were proposed by G. Matheron [START_REF] Matheron | Traité de géostatistique appliquée, Tome I[END_REF][START_REF] Matheron | Traite de Geostatistique Appliquee, Tome II: Le Krigeage[END_REF][START_REF] Matheron | Principles of geostatistics[END_REF]. Used initially in geostatistics, kriging was later coupled with calculation methods for the resolution of design problems by [START_REF] Sacks | Designs for computer experiments[END_REF][START_REF] Sacks | Design and analysis of computer experiments[END_REF]. Simpson [START_REF] Simpson | Comparison of response surface and kriging models for multidisciplinary design optimization[END_REF] gave a review of various metamodels for multidisciplinary optimization. Additional information can be found in [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF].

This chapter focuses on the gradient-enhanced cokriging metamodel [START_REF] Chung | Using gradients to construct cokriging approximation models for high-dimensional design optimization problems[END_REF][START_REF] Koehler | Computer experiments[END_REF][START_REF] Morris | Bayesian design and analysis of computer experiments: use of derivatives in surface prediction[END_REF]: Section 4.2 presents the principles and fundamentals of kriging and cokriging; then, Sections 4.4 and 4.5 deal with the construction of the cokriging metamodel in the case of ordinary cokriging. Finally, Section 4.6 discusses the problem of the determination of the cokriging parameters.

Notations

We will use the following notations:

1. x (i) , i ∈ 0, n s denotes a point in the design space D.

(x (i) , i ∈ 1, n s is one of the n s sample points, and x (0) is an arbitrary point in the design space, which may or may not be a sample point.) 2. Y (x (i) ) and Y (x (i) ) denote respectively the response of the analytical function (or the response of the mechanical model) and the approximate response given by the metamodel at point

x (i) , i ∈ 0, n s . 3. R(x (i) , x ( j)
) is a correlation function expressing the correlation relation between points x (i) and x ( j) , (i, j) ∈ 0, n s 2 .

Principle

The principle of the cokriging metamodel is similar to that of the kriging metamodel. One defines a random process, associated with the deterministic response of the objective function, which is the sum of two components (Equation 3): a linear model µ which represents the trend of the data, and a departure Z from this linear model which represents the fluctuations around the trend.

Y (x (0) ) = µ(x (0) ) + Z(x (0) ), ∀x (0) ∈ D (3) with E Z(x (0) ) = 0, Var Z(x (0) ) = 0, ∀x (0) ∈ D (4) and cov Z x (i) , Z x ( j) = 0, ∀ x (i) , x ( j) ∈ D 2 (5)
E, Var and cov are the classical statistical expected value, variance and covariance. µ is a deterministic function and Z is a stationary Gaussian process with a known stationary covariance. Depending on the definition of function µ, one can build different types of kriging or cokriging metamodels (simple kriging, where µ is the average of the values of the objective function at the sample points; ordinary kriging, where µ is an unknown constant; or universal kriging, where µ is a polynomial function).

Function Z has a zero expected value and its covariance structure is a function of a generalized distance among the sampled responses. The covariance structure can be written as:

∀ x (i) , x ( j) ∈ D 2 , cov Z x (i) , Z x ( j) = σ 2 R x (i) , x ( j) = σ 2 c i j (6)
where σ 2 is the variance of process Z. The correlation function can be a Gaussian function or a Matérn function [START_REF] Matérn | Spatial Variation[END_REF][START_REF] Stein | Interpolation of Spatial Data: some theory for kriging[END_REF].

In the context of cokriging, the variables are divided into a primary variable Y and N auxiliary variables W i . The auxiliary variables used to build a gradient-based cokriging metamodel are the components of the gradients (Equation 8). The construction of the cokriging metamodel involves Equation 7 in addition to Equation 3.

∀i ∈ 0, N , ∀x (0) ∈ D, W i (x (0) ) = µ W i + Q i (x (0) ) (7)
where

∀i ∈ 0, N , ∀x (0) ∈ D, W i (x (0) ) = ∂Y ∂ x i x (0) (8)
and

∀i ∈ 0, N , ∀x (0) ∈ D, E W i (x (0) ) = µ W i (9)
E Q i (x (0) ) = 0 (10) ∀ x (i) , x ( j) ∈ D 2 , ∀(k, l) ∈ 0, N 2 , ( 11
) cov Q k x (i) , Q l x ( j) = 0 (12) cov Q k x (i) , Z x ( j) = 0 (13) Var Q k x (i) = Var ∂ Z ∂ x k x (i) (14)
Thus, in the case of a cokriging metamodel, additional covariance relations involving the different variables must be introduced [START_REF] Morris | Bayesian design and analysis of computer experiments: use of derivatives in surface prediction[END_REF]:

cov ∂ Z ∂ x k x (i) , Z x ( j) = -σ 2 ∂ R ∂ x k x (i) , x ( j) = σ 2 c ki j (15)
.

cov Z x (i) , ∂ Z ∂ x k x ( j) = σ 2 ∂ R ∂ x k x ( j) , x (i) = σ 2 c k ji (16) cov ∂ Z ∂ x k x (i) , ∂ Z ∂ x l x ( j) = -σ 2 ∂ 2 R ∂ x k ∂ x l x (i) , x ( j) = σ 2 c kli j (17)
where

∀(i, j, k, l) ∈ 0, n s 2 × 1, n d 2 , ∀ x (i) , x ( j) ∈ D 2 , c i j =R x (i) , x ( j) = c ji (18)
c ki j = ∂ R ∂ x k x (i) , x ( j) = - ∂ R ∂ x k x ( j) , x (i) = -c k ji ( 19
)
c kli j = ∂ 2 R ∂ x k ∂ x l x (i) , x ( j) = c lki j = c kl ji (20)

The correlation function

In the following examples, we consider the correlation function to be the Matérn function [START_REF] Matérn | Spatial Variation[END_REF], defined by:

∀r ∈ R + , M (r, l) = 1 + √ 5 l r + 5 3l 2 r 2 exp - √ 5 l r (21)
where l denotes the correlation length (l > 0). More details on the Matérn function can be found in [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF].

In the general case, we use the correlation function defined by:

∀(i, j, k, l) ∈ 0, n s 2 × 1, n d 2 , ∀ x (i) , x ( j) ∈ D 2 , R x (i) , x ( j) = n d ∏ k=1 M (|x (i) k -x ( j) k |, l k ) (22)
where n d is the number of design variables in the problem; x (i) k and l k are the k th component of point x (i) and the associated correlation length respectively.

The first derivative of the Matérn function is nonsymmetrical. Therefore, one has c ki j = c k ji with (k, i, j) ∈ 1, n d × 0, n s 2 (see Equation 19).

Construction of the cokriging metamodel

The method used to build the cokriging metamodel is similar to that used to build a kriging metamodel, which leads to an estimator known as the best linear unbiased predictor as follows:

1. The objective is to determine an estimator Y of the random process Y given the linear predictor defined by Equation 23with x (0) ∈ D; this is equivalent to determining the λ 0i (∀i ∈ 1, n s ) and the λ i j (∀(i, j) ∈ 1, n s × 1, n d ) ); 2. the linear predictor must satisfy the unbiasedness conditions of Equation 24; 3. the linear predictor must minimize the mean square error defined in Equation 25subject to the previous unbiasedness conditions.

Y x (0) = n s ∑ i=1 λ 0i x (0) Y x (i) + n d ∑ j=1 n s ∑ i=1 λ i j x (0) ∂Y ∂ x j x (i) (23) E Y x (0) = E Y x (0) (24) MSE Y (x (0) ) = E Y x (0) -Y x (0) 2 (25)

Ordinary CoKriging (OCK)

For simplicity's sake, as suggested in [START_REF] Dunham | Design by natural selection[END_REF], we chose to limit ourselves to an ordinary cokriging metamodel. In this case, the deterministic function µ is an unknown constant. The construction method of Section 4.4 leads to the following estimator (from Equation 26)

Y x (0) = β c µ(x (0) ) + r T c0 C -1 c (Y sc -X c β c ) Z(x (0) ) (26)
where

β c =(X T c C -1 c X c ) -1 X T c C -1 c Y sc r c0 = c 01 c 02 • • • c 0n s c 110 c 210 • • • c n d n s 0 X c = 1, 1, • • • , 1, 0, 0, • • • , 0 T C c = C C d C d T C dd
and, using the notations presented in Equations 18-19-20:

(C) i j = c i j C d =        c 111 c 211 • • • c n d 11 c 112 • • • c n d 1n s c 121 c 221 • • • c n d 21 c 122 • • • c n d 2n s c 131 c 231 • • • . . . . . . . . . c 1n s 1 • • • c n d n s n s        C dd =       c 11 dd c 12 dd • • • c 1n s dd c 21 dd c 22 dd • • • . . . . . . . . . . . . c n s 1 dd c n s 2 dd • • • c n s n s dd       with ∀(k, l) ∈ 1, n s 2 , c kl dd =       c 11kl c 12kl • • • c 1n d kl c 11kl c 12kl • • • . . . . . . . . . . . . c n d 1kl c n d 2kl • • • c n d n d kl      
In the cokriging case, vector Y sc contains both the responses and the gradients of the objective function at the sample points, whereas in the kriging case vector Y s contains only the responses of the function. Equations 15-16-17 lead to the linear estimator as a function of the correlation vector and matrices r c0 and C c . Therefore, the resulting expression does not contain σ 2 .

Thus defined, the model can supply approximate responses of the objective function at every point of the design space. In our case, the response of the function is considered to be deterministic and we end up with a cokriging interpolation model. This type of metamodel has other advantages: for example, it provides statistical information on the process (its expected value and its variance). Due to the use of the unbiasedness condition, the expected value of Y is given by the trend model µ and the mean square error of Y [START_REF] Sacks | Design and analysis of computer experiments[END_REF]:

MSE Y x (0) = σ 2 1 -1 n s r c0 0 X T c X c C c 1 n s r c0 (27)
where 1 n s is a vector which contains n s ×1. The correlation matrix and vector appear to depend on the correlation function, i.e. on the correlation lengths. A specific strategy can be used to estimate values of the correlation lengths and σ 2 .

Estimation of the parameters

The model's parameters (such as the characteristic correlation lengths l, the variance σ of the random process Z or the regression coefficients) can be determined by maximizing the likelihood function [START_REF] Mardia | Maximum likelihood estimation of models for residual covariance in spatial regression[END_REF]. We use this technique, which relies on the maximization of the density of the observed values Y sc , to determine l and σ . The density can be viewed as a function L of parameters l and σ :

L(l, σ 2 ; Y sc ) = (2πσ 2 ) -n s (n d +1)/2 |C c (l)| -1/2 ×exp - 1 2 (Y sc -X c β c ) T C c (l) -1 (Y sc -X c β c ) (28)
where C c is the correlation matrix. This matrix depends on the correlation lengths {l 1 , l 2 , ..., l n d }, which constitute a vector l.

The maximization of the likelihood function can be expressed as:

( l, σ 2 ) = argmax l,σ 2 L(l, σ 2 ; Y sc ) (29)
This problem can usually be solved by minimizing the loglikelihood.

The variance σ 2 can be determined analytically through the derivation of the likelihood function:

σ 2 = 1 n s (Y sc -X c β c ) T C c (l) -1 (Y sc -X c β c ) (30)
One can also use an optimizer to determine the correlation length numerically. This method has some drawbacks [START_REF] Mardia | On multimodality of the likelihood in the spatial linear model[END_REF][START_REF] Sasena | Flexibility and efficiency enhancements for constrained global design optimization with kriging approximations[END_REF][START_REF] Warnes | Problems with likelihood estimation of covariance functions of spatial gaussian processes[END_REF]: in many cases with very few points, the loglikelihood is monotonous; the correlation matrix often suffers from conditioning problems which make it difficult to find a minimum. When such problems arise, one sets the parameters (particularly the correlation lengths) to fixed values.

Examples of the construction of metamodels

In this chapter, several examples based on analytical functions and mechanical test cases are presented: first, a onedimensional analytical function is used to build metamodels and, in order to compare the quality of these metamodels, different criteria are presented in (Section 5.1.1). Then, kriging and cokriging metamodels are built using two two-dimensional analytical functions with and without anisotropy (Sections 5.1.2 and 5.1.3). Finally, two mechanical problems are used to study the performance of the coupled multiparametric/cokriging strategy approach (Sections 5.2 and 5.3). For simplicity's sake, only problems with two design variables are being considered.

Analytical applications

In this section, kriging and cokriging are applied to one-and two-dimensional analytical test functions. The abbreviations OK and OCK will be used to designate Ordinary Kriging and Ordinary CoKriging respectively.

Case of a one-dimensional test function

First, we applied the two types of metamodels to the analytical function y(x) = exp(-x/10) cos(x) + x/10. We used 6 sampled responses of the analytical function to build the OK metamodel, and an additional 6 sampled derivatives to build the OCK metamodel. The correlation function was the Matérn function. The sample points were obtained using Latin Hypercube Sampling (LHS) [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF]. 7 illustrate the capability of the cokriging metamodel to interpolate not only the values of the responses, but also the sampled derivatives. In this example, the cokriging metamodel performed better: Table 2 shows that the R 2 criterion is better with cokriging than with kriging and the other criteria confirm that in this case cokriging leads to a better approximation than kriging. This statement remained true as long as we worked with only a few points. One can also observe that for a relatively smooth function the kriging metamodel converges quickly toward a good approximation of the actual function when the number of points becomes large enough.

For these two metamodels the correlation length was determined by maximizing the likelihood following the strategy introduced in Section 4.6. Figure 8 shows the loglikelihood as a function of the correlation length. The characteristics of the two metamodels are summarized in Table 2 The criteria Q i , which compare the actual response and the responses of the metamodels at n c points (n c >> n s ) of a regular grid with were calculated as follows:

Q 1 = sup i∈{1,2,...,n c } e i = e ∞ Q 2 = n c ∑ i=1 e i = e 1 Q 3 = Q 2 n c (the average of Q 2 )
where ∀i ∈ {1, 2, ..., n c }, Based on the statistical information obtained from the kriging and cokriging metamodels, one can derive confidence intervals. The two diagrams of Figures 9 and10 show the 95% Confidence Intervals (CIs) obtained with Expression 31. In these types of metamodels, the size of the confidence envelopes is determined mainly by the distance between each pair of neighboring points. For the same number of sample points, the cokriging metamodel provides narrower confidence intervals, especially close to the sample points. One could use the information derived from the variance or the confidence interval of random process Y to choose additional sample points (e.g. the points where the variance is maximum) in order to enrich the database used in constructing the metamodels. Such a strategy is not considered in this paper. Fig. 10: The confidence envelopes for the cokriging metamodel

e i = Y x (i) -Y x (i)
CI ± (x) = Y (x) ± 2 MSE Y (x) (31)

Case of a two-dimensional test function

The same two types of metamodels were used to approximate an analytical function of two variables. In order to illustrate the performance of the cokriging metamodel, we chose a very irregular function: the six-hump camel back function

(∀(x 1 , x 2 ) ∈ [-2, 2]×[-1, 1], f (x 1 , x 2 ) = (4 - 2.1x 2 1 + x 4 1 )x 2 1 + x 1 x 2 + 4(x 2 2 -1)x 2 2
). The two metamodels were constructed using 16 evaluations of the function for the kriging metamodel and an additional 16 evaluations of its gradients for the cokriging metamodel. In both cases, the correlation function was the Matérn function. First, the correlation length was considered to be the same for two components of the design space. Then, different correlation lengths were considered. The characteristics of these two metamodels are given in Table 3. For this 2D test function, the cokriging metamodel led to a relatively accurate approximation of the actual function using only a few sample points. Taking into account the gradients, we were able to develop more efficient approximate models. However, for a problem involving the evaluation of mechanical responses, the computation cost of determining n s responses and n s gradients is obviously higher than for n s responses only. (In the former case, due to the use of finite differences to calculate the gradients, the construction of the metamodel requires 3n s evaluations of the mechanical model.)

Now let us take another approach which consists in constructing metamodels using the same number of evaluations. In the following example, we used the same six-hump function, but the kriging metamodel was constructed based on the responses at 27 sample points (Figure 13a), while the cokriging metamodel was still constructed using the responses and the gradients at 9 sample points (Figure 13b). Table 4 shows that the quality of the predictions given by the two metamodels was similar, but a very important advantage of the cokriging metamodel was that it led to more zones where a minimum could be found than the kriging metamodel. If the problem involved the calculation of mechanical responses, the computation cost associated with the 9 points and their gradients would be much lower than that associated with the 27 points thanks to the multiparametric strategy, as shown in Table Now, let us take into account anisotropy. In the context of using correlation functions with the correlation length as a parameter, anisotropy is taken into account by allowing different correlation lengths for each design variable. We used the same two-dimensional example as in the previous section: the kriging and cokriging metamodels were constructed from the same data with the same number of function evaluations for the two models. The results were compared with those of the corresponding metamodels with a single correlation length.

Figures 14a and14b show the kriging and cokriging metamodels whose results are to be compared with those of Figures 13a and13b. The quality of these predictions is given in Table 5 (to be compared with the values of the quality criteria of Table 4). 14a-14b)

Table 6 shows the estimated correlation lengths obtained with the maximized likelihood with a single correlation length, then with anisotropy. For simplicity's sake, the sample points and the evaluations were normalized in order to carry out the kriging and cokriging metamodel construction process with standardized variables. The correlation lengths given in Table 6 refer to the variables after this transformation. (a) Kriging metamodel (to be compared with the surface of Figure 13a) (b) Cokriging metamodel (to be compared with the surface of Figure 13b) Fig. 14: The kriging and cokriging metamodels obtained with 27 evaluations of the analytical function with anisotropy. (To be compared with Figure 13; the isolines can be seen in Figure 23 The correlation lengths for the construction of the metamodels of the two-dimensional analytical function

Discussion of the example with anisotropy

The previous results confirmed the results obtained with a single correlation length: the quality of the predictions given by the kriging and cokriging metamodels with the same amount of data (i.e. the same number of evaluations of the function) was very similar. The introduction of anisotropy improved the approximation for both types of metamodels (see Tables 4 and5). This suggests that it is advisable to use a separate correlation length for each design variable.

Application to a contact problem with friction

In this section, we study the construction of the two previous metamodels (kriging and cokriging) using the multiparametric strategy presented in the first part of the paper along with the responses of mechanical models. The test cases considered are the three-squares example discussed in Section 3.3.2 plus a shrink-fit test case. The first example, which has a very smooth response surface, enables us to construct and study simple metamodels. The second example enables us to illustrate the strategy on a more complex and realistic response surface presenting local optima.

In order to illustrate the performance of the coupled multiparametric/cokriging strategy, we carried out two studies: one with a fixed number of mechanical calls (Section 5.2.1), and the other with a fixed quality of the metamodels (Section 5.2.2).

The case of a fixed number of mechanical calls

Each metamodel was constructed using 15 mechanical evaluations. The kriging metamodel was defined using 15 values of the force; the cokriging metamodel was defined using 5 values of the force and 5 gradients (with two components each, totaling 10 evaluations). The sample points were obtained through Latin hypercube sampling. Our calculations led to the two response surfaces of Figures 15a and15b, which are to be compared with the actual response surface of Figure 5.

Tables 7 and8 give the characteristics of the calculations and of the metamodels. With this example, the cokriging metamodel led to a much better approximation of the mechanical model than the kriging metamodel. Moreover, this result was obtained at a lower computation cost. The kriging and cokriging response surfaces of the three-squares example using 15 evaluations (15 responses for the kriging case; 5 responses and 5 gradients for the cokriging case). The isolines can be seen in Figure 24 of the appendix. Table 9: The metamodels constructed using full factorial sampling calls, the cost with the cokriging metamodel was lower. The simplicity of the response surface of this test case enabled the metamodels built with the multiparametric strategy to be studied easily and visually. In order to assess the performance of the cokriging metamodel, we also introduced a more irregular response surface associated with a shrink-fit test case.

Application to a shrink-fit problem

In order to study a case with a more irregular response surface capable of illustrating the performance of the cokriging metamodel, we introduced an additional test case: a shrinkfit assembly problem in 2D. Figure 16 shows the geometry of this problem. The reference problem consists of a cylindrical shaft inserted into a perforated rectangular plate (h 1 = 0.3m, h 2 = 0.2m, Young's modulus E = 80GPa, Poisson's coefficient ν = 0.3). The shaft is made of a perfect material, but presents geometric defects: a diameter slightly larger than its nominal value and a small eccentricity. The plate was divided into 16 substructures, each meshed with six-node triangular elements. The shaft itself was a substructure, also meshed with six-node triangular elements, with the same material characteristics and a nominal diameter of 5cm. Figure 17 shows the nominal geometry and the actual geometry of the shaft. The eccentricity defect was set to e = 0.5µm and the excess diameter was set to d r = R -R = 1µm. The loading was applied in two stages: first, the shaft was mounted in the frame; then, a uniform vertical pressure was applied progressively to the shaft up to a maximum 250MPa. The mechanical solution of this problem was obtained through a quasi-static resolution. Figure 18 shows the reference response surface plotted on a 20×20 regular grid: one can observe that for this problem the objective function has several local minima and is more irregular than that of the three-squares test case.

h 2 h 1 Ω 1 Ω 2 Ω 3 Ω 4 Ω 5 Ω 6 Ω 8 Ω 9 Ω 10 Ω 11 Ω 13 Ω 14 Ω 15 Ω 16 Ω 7 Ω 12
The kriging and cokriging metamodels of this test case were constructed from evaluations of the objective function at 30 and 10 sample points respectively. These sampled functions were obtained using the Latin Hypercube Sampling procedure. Figure 19 and Table 10 show the resulting metamodels and their characteristics.

These results show that the quality of the two metamodels is very similar. The cokriging metamodel was less expensive to build than the kriging metamodel, due largely to This paper proposes a new multilevel model optimization strategy based on a multiparametric approach and a cokriging metamodel and presents a detailed study of the cost of constructing this type of metamodel. Based on two toolsthe Multiparametric strategy and a gradient-based cokriging metamodel -the proposed approach allows us to obtain accurate approximations of objective functions. Morevoer, due to the use of the M.P.S., approximations are built with a significant reduction of the computational time in comparison with classical mechanical solver (i.e. without reinitialization of the solving algorithm). On an two-variables academic test example one can reach a gain about 3 to 4 in terms of computational time. The closer the sample points, the greater the gain. Obviously, the gain would be even greater in the case of 3D problems [START_REF] Roulet | A parallel strategy for the multiparametric analysis of structures with large contact and friction surfaces[END_REF].

This conclusion led to the construction of a gradientbased cokriging metamodel. The final conclusion of this study is the following:

-With similar numbers of sample points, cokriging leads to more accurate results than kriging, but this is not very useful in the mechanical context due to its higher computation cost. (The construction of the cokriging metamodel requires three times as many calculations as the kriging metamodel in order to obtain the responses and the gradients.) -For a given quality of the results, cokriging enables a significant reduction in the number of mechanical calls compared to kriging.

-For a given quality of the results, cokriging leads to a significant reduction in computation cost compared to kriging. -Like Kriging, Cokriging provides accurate locations of the local optima. -In all cases involving more than one design variable, the use of the anisotropy principle is crucial for optimizing the quality of the metamodels.

These results pave the way for future works on larger and more complex assembly problems in order to deal with more realistic industrial structures, such as two-and threedimensional problems with many design variables. Besides, future developments will be also to compare cokriging metamodel with other classical models such as gradient-based Radial Basis Function [START_REF] Leary | Global approximation and optimization using adjoint computational fluid dynamics codes[END_REF] and Artifical Neural Network [START_REF] Kc Giannakoglou | Aerodynamic shape design using evolutionary algorithms and new gradient-assisted metamodels[END_REF]. Another objective will be to complete the development of the first and second levels of the multilevel model optimization strategy, i.e. the exploratory phase to localize the optimum approximately and the subsequent precise optimization through new mechanical calculations using the multiparametric strategy along with the results from the first level.

A Complementary figures

Fig. 2 :

 2 Fig. 2: Schematic representation of the LATIN process

Fig. 4 : 2 FFig. 5 :

 425 Fig. 4: The geometry of the problem

Fig. 7 :

 7 Fig. 6: Ordinary kriging and cokriging

Fig. 8 :

 8 Fig. 8: The log-likelihood vs. the correlation length

2 supj∈{1, 2 ,

 22 ...,n c } Y x ( j) 2 and e = e 1 e 2 • • • e n c with • 1 and • ∞ being the L 1 norm and infinity norm respectively.

Fig. 9 :

 9 Fig. 9: The confidence envelopes for the kriging metamodel

Fig. 11 :

 11 Fig. 11: The response surface of the actual function. (The isolines can be seen in Figure 20 of the appendix.)

Fig. 12 :

 12 Fig.12:The kriging and cokriging metamodels obtained with 16 evaluations (plus 16 gradients for cokriging) of the analytical function with only one correlation length. (The isolines can be seen in Figure21of the appendix.)

Fig. 13 :

 13 Fig.13: The kriging and cokriging metamodels obtained with the 27 evaluations of the analytical function (27 responses for kriging; 9 responses and 9 gradients for cokriging) with only one correlation length. (The isolines can be seen in Figure22of the appendix.)

  Fig.15:The kriging and cokriging response surfaces of the three-squares example using 15 evaluations (15 responses for the kriging case; 5 responses and 5 gradients for the cokriging case). The isolines can be seen in Figure24of the appendix.

Ω 17 µFig. 16 :

 1716 Fig. 16: The geometry of the shrink-fit assembly problem

Fig. 17 :

 17 Fig. 17: The nominal and actual geometries of the shaft

Fig. 18 :Fig. 19 :

 1819 Fig. 18: The reference response surface of the shrink-fit test case

Fig. 20 :Fig. 21 : 2 -Fig. 22 : 2 -Fig. 23 :Fig. 24 :

 202122222324 Fig. 20: Isolines and gradients of the six-hump camel back function presented in Section 5.1.2 and in Figure 11

Fig. 25 :

 25 Fig.25:The kriging and cokriging metamodels constructed with 30 mechanical evaluations associated with the response surfaces presented in Figure19

Table 1 :

 1 The mechanical calculations using the multiparametric strategy

		Eval.	Eval. gradients	Total
	Number of calculations	16	2×16	48
	Avg. CPU time per calc. Min. CPU time Max. CPU time	2.97s 0.17s 8.12s	1.53s 0.16s 4.83s	2.01s 0.16s 8.12s
	Nb. of iterations Avg. nb. of iterations (per calc.) Min. nb. of iterations Max. nb. of iterations	2,201 137 5 381	2,215 69.2 5 225	4,416 92 5 381
	Total CPU time CPU time w/o MPS Gain	47.5s 135s 2.84	48.9s 270s 5.52	96.4s 405s 4.22

  .

	OK	OCK
	R 2 RAAE RMAE 0.7162 0.8593 0.2949 Q 1 7.219 • 10 -2 2.758 • 10 -2 0.9666 0.1326 0.4427 Q 2 5.942 1.4124 Q 3 1.974 • 10 -2 4.693 • 10 -3

Table 2 :

 2 Characteristics of the two previous metamodels RAAE and RMAE designate the Relative Average Absolute Error and Relative Maximum Absolute Error criteria.

Table 3 :

 3 Characteristics of the two metamodels (Figure12a-12b) with a single correlation length

  1.

	OK	OCK
	R 2 RAAE RMAE 2.753 0.7818 0.2920 Q 1 0.3550 Q 2 16.35 Q 3 1.022 • 10 -2 1.165 • 10 -2 0.7512 0.3579 1.600 0.1199 18.64

Table 4

 4

	: Characteristics of the previous metamodels ob-
	tained with the same number of evaluations (27 sample
	points for kriging and 9 sample points for cokriging; see
	Figures 13a-13b)
	5.1.3 The case of two correlation lengths (anisotropy)

Table 5 :

 5 Characteristics of the two metamodels with anisotropy (Figures

  of the appendix.)

	Anisotropy	No	Yes	Figures
	OK (27 eval.) OCK (9 eval. & 9 grad.) 1.014 0.9564 0.7106 1.5218 0.9878 1.310	13a/14a 13b/14b

Table 6 :

 6

Table 7 :

 7 Characteristics of the mechanical calculations

Table 8 :

 8 Characteristics of the metamodels5.2.2 The case of a fixed quality of the metamodels

	Another way to study the cost of the metamodels consists
	in constructing kriging and cokriging metamodels of similar
	quality. In this case, we constructed the metamodels using
	sample points obtained through full factorial sampling. We
	chose to build the kriging metamodel with 9×9 samples and
	the cokriging metamodel with 3×3 and 5×5 samples. The
	results are shown in Table 9.
	5×5 full factorial sampling led to a cokriging meta-
	model with the same quality as the 9×9 kriging metamodel.
	As in the previous case with a fixed number of mechanical

Table 10 :

 10 Characteristics of the metamodels

	6 Conclusion

Acknowledgements This work was supported by the Agence Nationale de la Recherche as part as ANR-08-COSI-007-10: Distributed Multi-Disciplinary Optimization.