Generation of a cokriging metamodel using a multiparametric strategy

In the course of designing structural assemblies, performing a full optimization is very expensive in terms of computation time. In order or reduce this cost, we propose a multilevel model optimization approach. This paper lays the foundations of this strategy by presenting a method for constructing an approximation of an objective function. This approach consists in coupling a multiparametric mechanical strategy based on the LATIN method with a gradient-based metamodel called a cokriging metamodel. The main difficulty is to build an accurate approximation while keeping the computation cost low. Following an introduction to multiparametric and cokriging strategies, the performance of kriging and cokriging models is studied using one- and two-dimensional analytical functions; then, the performance of metamodels built from mechanical responses provided by the multiparametric strategy is analyzed based on two mechanical test examples.

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 [1]. Our approach belongs to that category in the sense that, contrary to the surrogate-based approach [2,3], 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 [3,4]). 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 [5][6][7] and iterative methods [8][9][10][11][12][13]) and multilevel model optimization (which includes hierarchical multilevel optimization [14,1] and imbricated multilevel optimization [15]). Similar strategies have also been used in the context of multidisciplinary optimization where each level of the optimization process concerns a specific discipline [16,9]. Parallel multilevel model optimization [17,18] 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 [19][20][21].
The multilevel modeling strategy we propose is based on multilevel model optimization. We focus on the two-level strategy illustrated in Fig. 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 [22]; 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 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 [23,24], on neural networks [25,26], on radial basis functions [27], on proper orthogonal decomposition [28], on cumulative interpolation [29], etc... We chose to use a particular class of metamodels called kriging approximations [30] and, more precisely, a cokriging metamodel [31] using derivatives [32]. These approximations are presented in Sect. 4. The interested readers could also refer to works on metamodeling strategies in structural optimization and multidisciplinary design optimization [33,34].

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 Sect. 3.1. The strategy is based on the LATIN method, which is presented in Sect. 3.2. The multiparametric strategy itself is presented in Sect. 3.3 and its performance is discussed using the example of a 2D mechanical benchmark problem with two design variables in Sect. 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 Ladevèze [35].
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 [36,37] or dual approach [38].

Short summary of the LATIN method for assemblies
Our resolution strategy is based on the LATIN algorithm developed by Ladevèze [35], 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 Fig. 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 [39] are part of equation group . The linear form of A d is due to the linear behavior of the substructures.
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 [40,35] for the LATIN method and in [39,41] 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 [42,43], then studied and applied to various types of problems in [44][45][46]. 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 Sect. 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  The geometry of the problem 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 [47]. 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 iterations 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 [45]. Figure 4 shows the geome- 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: 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: 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 three 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.5 s, but the addition of 16 gradients (32 evaluations) took only 48.9 s. 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 vari-ables jointly. The cokriging metamodel we use [31] 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 Krige, a mining engineer [48]. Subsequently, many mathematical enhancements were proposed by Matheron [49][50][51]. Used initially in geostatistics, kriging was later coupled with calculation methods for the resolution of design problems by [52,53]. Simpson [24] gave a review of various metamodels for multidisciplinary optimization. Additional information can be found in [54].
This chapter focuses on the gradient-enhanced cokriging metamodel [32,55,56]: Sect. 4.2 presents the principles and fundamentals of kriging and cokriging; then, Sects. 4.4 and 4.5 deal with the construction of the cokriging metamodel in the case of ordinary cokriging. Finally, Sect. 4.6 discusses the problem of the determination of the cokriging parameters.

Notations
We will use the following notations: 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 ) is a correlation function expressing the correlation relation between points x (i) and

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 (Eq. 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. and 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: where σ 2 is the variance of process Z . The correlation function can be a Gaussian function or a Matérn function [57,58].
In the context of cokriging, the variables are divided into a primary variable Y and n d auxiliary variables W i . The auxiliary variables used to build a gradient-based cokriging metamodel are the components of the gradients (Eq. 8). The construction of the cokriging metamodel involves Eq. 7 in addition to Eq. 3. where and Thus, in the case of a cokriging metamodel, additional covariance relations involving the different variables must be introduced [56]: where

The correlation function
In the following examples, we consider the correlation function to be the Matérn function [57], defined by: where l denotes the correlation length (l > 0). More details on the Matérn function can be found in [54].
In the general case, we use the correlation function defined by: where n d is the number of design variables in the problem; x (i) k and l k are the kth component of point x (i) and the associated correlation length respectively.
The first derivative of the Matérn function is nonsymmetrical. Therefore, one has

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 Eq. 23 with x (0) ∈ D; this is equivalent to determining ); (2) the linear predictor must satisfy the unbiasedness conditions of Eq. 24; (3) the linear predictor must minimize the mean square error defined in Eq. 25 subject to the previous unbiasedness conditions.

Ordinary CoKriging (OCK)
For simplicity's sake, as suggested in [19], we chose to limit ourselves to an ordinary cokriging metamodel. In this case, the deterministic function μ is an unknown constant. The construction method of Sect. 4.4 leads to the following estimator (from Eq. 26) wherê 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 and, using the notations presented in Eqs. 18, 19, 20: 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 [53]: 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 [59]. 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 σ : 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: This problem can usually be solved by minimizing the loglikelihood.
The variance σ 2 can be determined analytically through the derivation of the likelihood function: One can also use an optimizer to determine the correlation length numerically. This method has some drawbacks [60][61][62]: in many cases with very few points, the log-likelihood 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 (Sect. 5.1.1). Then, kriging and cokriging metamodels are built using two two-dimensional analytical functions with and without anisotropy (Sects. 5.1.2 and 5.1.3). Finally, two mechanical problems are used to study the performance of the coupled multiparametric/cokriging strategy approach (Sects. 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) [63]. Figures 6 and 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 For these two metamodels the correlation length was determined by maximizing the likelihood following the strategy introduced in Sect. 4.6. Figure 8 shows the log-likelihood as a function of the correlation length.
Classical numerical optimization (using function fmincon in MATLAB) led to correlation lengths equal to l O K = 0.2807 and l OC K = 0.3481 for the kriging and cokriging metamodels respectively. (In all the kriging and cokriging examples presented in this paper, the data are normalized in order to build the metamodels from reduced variables.) The characteristics of the two metamodels are summarized in Table 2. RAAE and RMAE designate the Relative Average Absolute Error and Relative Maximum Absolute Error criteria. 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: and e = e 1 e 2 · · · e n c with · 1 and · ∞ being the L 1 norm and infinity norm respectively. Based on the statistical information obtained from the kriging and cokriging metamodels, one can derive confidence intervals. The two diagrams of Figs. 9 and 10 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 ran- x 1 x 2 Fig. 11 The response surface of the actual function (The isolines can be seen in Fig. 20 of the appendix) dom 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.

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 , Figure 11 shows the actual response surface. 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 (Fig. 12). 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.).  Table 3 Characteristics of the two metamodels (Figure 12a-12b)  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 (Fig. 13a), while the cokriging metamodel was still constructed using the responses and the gradients at 9 sample points (Fig. 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  Fig. 13; the isolines can be seen in Fig. 23 of the appendix) 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 1.

The case of two correlation lengths (anisotropy)
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.   Figure 14a, b shows the kriging and cokriging metamodels whose results are to be compared with those of Fig. 13a, b. The quality of these predictions is given in Table 5 (to be compared with the values of the quality criteria of Table 4). 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.

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 and 5). 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 Sect. 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 (Sect. 5.2.1), and the other with a fixed quality of the metamodels (Sect. 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,  The isolines can be seen in Fig. 24 of the appendix which are to be compared with the actual response surface of Fig. 5. Tables 7 and 8 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 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 metamodel with the same quality as the 9 × 9 kriging metamodel. As in the previous case with a fixed number of mechanical 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  Fig. 16 The geometry of the shrink-fit assembly problem 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.3 m, h 2 = 0.2 m, Young's modulus E = 80 GPa, Poisson's coef- 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 sixnode 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 5 cm. 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 250 MPa. The mechanical solution of this problem was obtained through a quasi-static resolution. We considered two design variables: the friction coefficient μ between the shaft and the frame (μ ∈ [0.02, 0.7]) and the orientation angle of the eccentricity defect ϕ (ϕ ∈ [0, 360]). The objective function studied was the maximum Von Mises' stress in the structure. 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.
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  19 The kriging (a) and cokriging (b) metamodels response surfaces of the shrink-fit test case using 30 evaluations (30 responses for kriging and 10 responses and 10 gradients for cokriging; the isolines are presented in Fig. 25 of the appendix) to build than the kriging metamodel, due largely to the multiparametric strategy. For this test case, the response surfaces enabled a good localization of the global minimum.

Conclusion
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-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 [64]. This conclusion led to the construction of a gradient-based 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 [65] and Artifical Neural Network [66]. 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.
Acknowledgments This work was supported by the Agence Nationale de la Recherche as part as ANR-08-COSI-007-10: Distributed Multi-Disciplinary Optimization.

Fig. 25
The kriging and cokriging metamodels constructed with 30 mechanical evaluations associated with the response surfaces presented in Fig. 19. a Kriging isolines and gradients (associ-ated with the response surface of Fig. 19a. b Cokriging isolines and gradients (associated with the response surface of Fig. 19b