Christelle Vergé

Jérôme Morio
email: jerome.morio@onera.fr

Pierre Del Moral

Juan Carlos

Dolado Pérez

Probabilistic safety analysis of the collision between a space debris and a satellite with an island particle algorithm

Keywords: Rare event, sequential Monte Carlo, island particle models, debris, satellite, collision, adaptive splitting technique

Collision between satellites and space debris seldom happens, but the loss of a satellite by collision may have catastrophic consequences both for the satellite mission and for the space environment. To support the decision to trigger o a collision avoidance manoeuver, an adapted tool is the determination of the collision probability between debris and satellite. This probability estimation can be performed with rare event simulation techniques when Monte Carlo techniques are not enough accurate. In this chapter, we focus on analyzing the inuence of dierent simulation parameters (such as the drag coecient) that are set for to simplify the simulation, on the collision probability estimation. A bad estimation of these simulation parameters can strongly modify rare event probability estimations. We design here a new island particle Markov chain Monte Carlo algorithm to determine the parameters that, in case of bad estimation, tend to increase the collision probability value. This algorithm also gives an estimate of the collision probability maximum taking into account the likelihood of the parameters. The principles of this statistical technique are described throughout this chapter.

Introduction

On February 10 th 2009, active commercial satellite Iridium-33 and out of order Russian satellite Cosmos-2251 collided [START_REF] Kelso | Analysis of the Iridium 33-Cosmos 2251 collision[END_REF]. The impact produced more than 2000 trackable debris. Most of them may destroy any satellite, whether in use or not, they might encounter. The safest practice for satellites that encounter space debris is to avoid collision. Avoidance maneuvers are an ecient mean to reduce the collision probability between two orbiting objects, nevertheless they consume fuel reducing the operational lifetime of the satellite and they perturb the operational mission of the satellite. Consequently, satellite safety responsible teams have to take into account the operational mission prior to the denition of a collision avoidance maneuver and try to combine, whenever possible, planned station keeping maneuvers with collision avoidance maneuvers. Avoidance maneuvers are decided, among other parameters, based on the estimated collision probability. The orbital motion of the space objects is simulated using a simplied deterministic dynamical model that may be considered as an input-output function where the random inputs are, for instance, the position and the speed of the debris and of the satellite as well as other dynamic parameters as the drag coecient, and the output is the minimum distance between the debris and the satellite. The collision probability is then estimated on this output. This input-output function can be seen as a "black-box" with random inputs. Some parameters, denoted by a vector Θ, in black-box functions are implicitly set, such as parameters of the model (the drag coecient for instance) or of the input parametric model density, and their value inuences the collision probability estimation. These hypotheses are often assumed for simplication and computational reasons. From a risk analysis point of view, it is interesting to determine the variability of the collision probability w. r. t. the uncertainty on theses input parameters Θ or w. r. t. one particular parameter, and to quantify the impact of such tuning in the realization of a collision. Of course, dierent values of Θ can strongly modify rare event probability estimation and sometimes miss a risk situation. The issue of concern in safety would be to underestimate a risk because of a bad tuning of model parameters Θ. That is why in this paper we propose to estimate the law of the parameters Θ conditionally on a collision between the debris and the satellite. We develop in this chapter the SMC 2 (Sequential Monte Carlo Square) algorithm to estimate this kind of targeted laws introduced [START_REF] Chopin | SMC 2 : an ecient algorithm for sequential analysis of state space models[END_REF] to do ltering on hidden Markov models. We apply this island particle algorithm to debris satellite collision use case and analyse its results for the system safety.

Debris satellite collision simulation

We consider two space objects (a debris and a satellite) orbiting around an Earth centered inertial reference frame. Their geometry is assumed spherical (i.e. the objects have a high tumbling motion when compared with their orbital period) and we assume that we perfectly know the radius of such sphere and the mass of the objects. We wonder about the relative position of the two satellites and ask whether the distance between the two objects could be smaller than a conict distance T during the given time span I. To model the orbital motion of both space objects, we consider a general perturbation approach where the original equations of motion are replaced with an analytical approximation that captures the essential character of the motion over some limited time interval, which also enables analytical integration of the equations. SGP4 model [START_REF] Miura | Comparison and design of Simplied General Pertur-University[END_REF] is used to propagate the trajectories of debris and satellite according to the time. At time t, the space objects will be represented by their 6-dimensional state vectors s 1 (t) and s 2 (t), i.e. their 3-dimensional position vectors r 1 (t) and r 2 (t) and their 3-dimensional speed vectors v 1 (t) and v 2 (t) such that s i = ( r i , v i ). The initial conditions in the proposed example is dened in terms of two line elements (TLE), similar to those provided by NORAD (North American Aerospace Defense Command), as the SGP4 model is used for the orbital propagation of the considered objects. The initial condition value is denoted s m i at a given time t m i . SGP4 enables us to propagate the orbit of both space objects through time, denoted by a scalar continuous function ν such that ∀i ∈ {1, 2}, ∀t ∈ I, s i (t) = ν( s m i , t m i , t),

δ = min t∈I { r 2 -r 1 (t)}.
The function of time t ∈ I → r 2 -r 1 (t) makes δ available through numerical optimisation in a deterministic approach. In fact, the position and velocity of space objects are estimated from more or less imprecise measurements. While the measurement means used for satellites (e.g. GPS, laser) result in a reasonable orbital accuracy (e.g. several tens of meters) the measurement means used for debris and uncooperative space-objects (e.g. mainly radar and telescopes) could result in quite imprecise orbits (e.g. several hundred of meters or few kilometers). This lack of accuracy will depend on a great number of factors. TLEs sum up this information and feed the models with the couple ( s m i , t m i ) for i = 1, 2, but to cope with their uncertainty, we have added independent and identically distributed Gaussian noises to the model inputs s m i . Debris satellite conict may be modelled as an input-output function where:

• the input X represents the position and the speed of the debris space (the position and the speed of the satellite are assumed to be known). X is a 6dimensional multivariate normal random vector of mean (Θ 1 , Θ 2 , Θ 3 , Θ 4 , Θ 5 , Θ 6 ) t and covariance matrix is equal to the identity matrix dened on a measurable space (X, X ). The means corresponds to the debris measurement errors on its position and speed; • the input-output function φ enables to propagate the debris and satellite trajectories with the SGP4 model during I. The input-output code includes the transformation that allows to switch from the standard space of the input to the physical space in which evolve the satellite and debris position and speed. The function φ is a continuous positive scalar function φ : R 6 → R and is static; • the error on the drag coecient which is considered inside the function φ is also random and follows a normal distribution with mean Θ 7 and variance 1; • the output Y is the minimum distance between the debris and the satellite during I. We assume that it is a positive random variable.

The complete set of model parameters is summed up in the vector

Θ = (Θ 1 , Θ 2 , Θ 3 , Θ 4 , Θ 5 , Θ 6 , Θ 7 ) t .
The quantity of interest on the output Y is the probability

P(Y < T ) = P(φ(X) < T ) .
When the event {φ(X) < T } is rare relatively to the available simulation budget (which is often the case in safety and reliability issues), dierent algorithms described in [START_REF] Sobol | A Primer for the Monte Carlo Method[END_REF], [START_REF] Bucklew | Introduction to Rare Event Simulation[END_REF], [START_REF] Rubinstein | The Cross-Entropy Method: a unied approach to combinatorial optimization[END_REF], [START_REF] Zhang | Nonparametric importance sampling[END_REF], [START_REF] Bjerager | Methods for structural reliability computation[END_REF], [START_REF] Botev | Ecient Monte-Carlo simulation via the generalized splitting method[END_REF], [START_REF] Cérou | Sequential Monte Carlo for rare event estimation[END_REF] have notably been proposed to estimate accurately its probability.

Basics of safety analysis

In the present chapter, one focuses on the case where the law of X is uncertain and depends upon unknown parameters. We assume that X is distributed according to a well known parametric model and its parameters, denoted by a random vector Θ, have a probability density ν. We also suppose that Θ has a density f Θ w. r. t. a dominating measure of reference λ, that is

ν(dθ) = f Θ (θ) λ(dθ) .
In the application considered here, X is a random vector with a multivariate normal distribution, and Θ describe the mean of X. It corresponds notably to realistic applications where it is not always possible to evaluate accurately the density of input parameters. This formalism enables thus to consider a large range of input probability density function.

The probability of interest P(Y < T ) depends of course on Θ and thus on the distribution ν. In safety applications, it is important to estimate a superior bound of the rare event probability P(Y < T ) taking also into account the prior on Θ. The prior on Θ is important since unrealistic bad tuning values of Θ which lead to high probabilities P(Y < T ) are not relevant. The idea of this chapter is thus to determine the distribution of Θ conditionally on the fact that Y does not exceed the threshold T . This distribution, denoted by π, will be referred to in the sequel as the target law.

In the further development, when there is no confusion, we sometimes write

P(Y < T |θ) instead of P(Y < T |Θ = θ).
Note that using the Bayes' formula, the target law can be written

π(dθ) = 1 P(Y < T ) P(Y < T |θ)ν(dθ) . (1)
We propose in this paper a SMC (Sequential Monte Carlo) algorithm which evolves according to iterative selection and mutation steps, and which approximates π when the number of particles gets large. This algorithm requires the estimation of P(φ(X) < T |Θ = θ) for dierent settings of parameter θ. For that purpose, we describe the splitting algorithm that enables us to estimate this probability with accuracy. 4 The SMC 2 algorithm

Principle

The SMC 2 algorithm is based on the use of two sets of particles to iteratively approach π. The rst set of particles is dened on the parameter Θ and the second set of particles is useful to estimate the probabilities P(Y < T |θ). The complete demonstration of interacting particles systems (IPS) convergence and the link with Feynman-Kac framework is given in [START_REF] Moral | On the Concentration Properties of Interacting Particle Processes, volume 3[END_REF]. Dene T 1 , T 2 , ..., T n = T a serie of decreasing thresholds and denote for all i ∈ 0, n

π i = 1 P(Y < T i ) P(Y < T i |θ)ν(dθ) .
The target law is of course π = π n . The probability law π n is thus proportional to

π n ∝ P(Y < T n |θ) ν(dθ) π n ∝ H n (θ) ν(dθ) ,
where H n (θ) = P(Y < T n |θ). The term H n (θ) can be expressed as a product of conditional probabilities

H n (θ) =   n-1 p=1 P(Y < T p+1 |Y < T p , θ)   × P(Y < T 1 |θ) = n-1 p=0 h p (θ) , (2)
with

     h p (θ) = P(Y < T p+1 |Y < T p , θ) h 0 (θ) = P(Y < T 1 |θ) .
In this notation, we have

π n ∝ n-1 p=0 h p (θ) ν(dθ) . (3)
One can also remark that H p = H p-1 ×h p-1 and consequently the link between π p+1 and π p can be written on the following way

π p+1 = ψ hp (π p ), (4)
where ψ hp is the so-called Boltzmann-Gibbs transformation. Let P(E) be the set of probability measures on E. For all positive bounded function G, the Boltzmann-Gibbs transformation Ψ G : P(E) → P(E) is dened for all µ ∈ P(E) such that µ(G) = G(x) µ(dx) > 0 by

Ψ G (µ)(dx) := 1 µ(G) G(x) µ(dx).
If one assumes that it is possible to determine a Markovian kernel M p that let π p invariant (which is not restrictive using, for example, a stage of the acceptance/rejection of Metropolis-Hastings algorithm) we have

π p = (π p M p )(dθ) = π p (dθ )M p (θ , dθ). (5)
This yields the evolution equation

π p+1 = ψ hp (π p )M p+1 . (6)
Equation 6 may be cast in the Feynman-Kac framework and then, each measure π p can be approximated by an IPS which evolves with selection steps related to the so called potential functions h p and mutation steps related to the Markov kernel M p . Denote by {(θ 1 p , ..., θ N 1 p )} n≥0 a system of N 1 particles.

                         θ 1 0 . . . . . . θ N 1 0                          Selection -----→                          θ 1 0 . . . . . . θ N 1 0                          Mutation ------→                          θ 1 1 . . . . . . θ N 1 1                          Selection -----→                          θ 1 1 . . . . . . θ N 1 1                          Mutation ------→
The selection stage consists in sampling { θ i p } N 1 i=1 independently according to the probability measure ψ hp , i.e. selecting the particles {θ i p } N 1 i=1 with probabilities proportional to their weights {h p (θ i p )} N 1 i=1 . The mutation stage consists in updating the selected particles conditionally independently using the Markov kernel M p+1 that let π p+1 invariant. This step enables to increase the diversity of θ p without changing its probability law, that is already close to π p+1 . The Feynman-Kac theory Moral et al. [2012] ensures that at each transition stage p :

1

N 1 N 1 i=1 δ θ i p -----→ N 1 →+∞ π p .
Thus, at the end of the n th transition stage, the system of particles converges to the target law π n so that

1 N 1 N 1 i=1 δ θ i n -----→ N 1 →+∞ π n .
Nevertheless, the knowledge of h p is required to apply the dierent selection/mutation stages. In practice, h p (θ i p ) is not analytically computable but can be estimated by dening a new set of particles {ξ i,j p } N 2 j=1 on the random variable X conditionally to the dierent thresholds T p and associated to each θ i p .

Description

Consider {θ i 0 } N 1 i=1 generated with probability law ν. At iteration k of the algorithm, with k ≥ 1, we assume that particles {θ i k } N 1 i=1 are available and then, the interacting island algorithm consists in two iterative type stages:

• Selection stage The selection stage consists in choosing randomly and independently N 1 particles amongst {θ i k } N 1 i=1 with probabilities proportional to their weights {h k (θ i k )} N 1 i=1 . Thus, the particles with low weights are killed whereas those with high weights are multiplied. The number of particles is kept constant in this stage and a new set particles { θ i k } N 1 i=1 can be dened. Remind that the potential functions h k are dened by:

     h k (θ i k ) = P(Y < T k+1 |Y < T k , θ = θ i k ), k ≥ 1 h 0 (θ i 0 ) = P(X < T 1 |θ = θ i 0 )
These quantities have to be computed.

• Mutation stage Even if the number of particles is still equal to N 1 , some particles have been duplicated, so we apply a Markov kernel to increase the diversity of the particles. Building a π k+1 -reversible transition kernel that let π k+1 invariant is the objective of mutation stage. For that purpose, the acceptance/rejection step of the Metropolis-Hastings algorithm Tierney [1994] is useful. This approach results in the exploration of Θ space set without changing the θ i k 1≤i≤N 1 distribution and the increase of the particle diversity. A new particle θ i k is proposed with a ν-reversible kernel Q. The acceptation rate of a new proposal is consequently 1 ∧

H k+1 ( θ i k ) H k+1 ( θ i k )
.

If H k+1 ( θ i k ) > H k+1 ( θ i k )
, the proposal θ i k is automatically accepted and replaces θ i k in the set of current particles. Otherwise, the proposal θ i k is accepted with probability

H k+1 ( θ i k ) H k+1 ( θ i k )
. This acceptance/rejection procedure is repeated N app times to decrease the correlation between the particles. At the end of this stage, a new set of particles {θ i k+1 } N 1 1=i can be dened.

Mutation and selection stage are applied n times until reaching the target threshold T n . At the end of the algorithm, the particles {θ i n } N 1 1=i provides an estimate of π n :

π N 1 n = 1 N 1 N 1 i=1 δ θ i n . For i ∈ [1, N 1 ] and k ∈ [0, n],
the point is to estimate each probability {h l (θ i k )} 1≤l≤k . It can be done with another interacting particle system (also called, in that case, sequential Monte Carlo, importance splitting, subset simulation or subset sampling). It is a rare event estimation technique which considers the estimation of several conditional probabilities that are easier to evaluate than estimate only one probability through a very tough simulation. Its principle is also based on selection and mutation stages. Let us dene {ξ i,j 0 } N 2 j=1 with probability density f X|θ i j w.r.t. λ X , i.e. the density of X knowing

Θ = θ i k .
At iteration l of the algorithm, we assume that particles {ξ i,j l } N 2 j=1 are available and then IPS consists in two iterative stages:

• Selection stage The selection stage consists in choosing randomly and independently N 2 particles amongst the particles {ξ i,j l } N 2 j=1 which are above T l . The particles which have not reached the threshold T l are thus killed. The number of particles is kept constant, and a new set of particles { ξ i,j l } N 2 j=1 can be dened.

• Mutation stage The mutation stage is patterned with acceptance/rejection principle using the Metropolis-Hastings algorithm [START_REF] Tierney | Markov chains for exploring posterior distributions[END_REF]. A new particle ξ i,j l is then proposed with a Markov kernel Q. If φ( ξ i,j l ) < T l , then the proposal is accepted with probability 1 ∧

f X|θ i l ( ξ i,j l ) Q( ξ i,j l , ξ i,j l ) f X|θ i l ( ξ i,j l ) Q( ξ i,j l , ξ i,j l )
and ξ i,j l replaces ξ i,j l in the set of current particles. If φ( ξ i,j l ) > T l , the proposal is automatically rejected and the particle ξ i,j l is remained. This acceptance/rejection procedure is repeated N app2 times to decrease the correlation between the particles. At the end of this stage, a new set of particles {ξ i,j l+1 } N 2 j=1 can be dened. An estimate

h l (θ i k ) of h l (θ i k ) = P(Y < T l+1 |Y < T l , Θ = θ i k )
is given by the ratio between the number of {ξ i,j l } N 2 j=1 particles such that φ(ξ i,j l ) < T l+1 and the total number of particles N 2 .

Mutation and selection stages are applied k times until reaching the target threshold T k . At the end of the algorithm,

H k+1 (θ i k ) = P(Y < T k+1 |θ i k ) is estimated by H k+1 (θ i k ) = k l=0
ĥl For a given particle θ i k , a complete set of particles {ξ i,j l } 1≤l≤k 1≤j≤N 2 is thus generated. An island particle is thus constituted of a particle θ i k and its associated {ξ i,j l } 1≤j≤N 2 1≤l≤k particle set.

The SMC 2 algorithm is described more precisely in Algorithm 1. Interacting particle system for probability estimation required in Algorithm 1 is developed in Algorithm 2.

The determination of Q and Q, in the general case, implies the use of Metropolis-Hastings algorithm. Nevertheless, if µ is a standard normal distribution, a transition from x to z dened with the following expression

x → z = √ 1 -a x + √ a W, (7)
where W ∼ N (0, 1) and a is scalar parameter such as a ∈ [0, 1], is µ-reversible.

In order to use equation 7 instead of Metropolis-Hastings algorithm, it is also possible to apply a transformation on the variables X or Θ so that they follow a standard normal PDF. Depending on the available information on the PDF of X, several transformations can be proposed [START_REF] Nataf | Distribution des distributions dont les marges sont données[END_REF], Pei-Ling and Kiureghian [1991], [START_REF] Lebrun | An innovating analysis of the Nataf transformation from the copula viewpoint[END_REF], [START_REF] Rosenblatt | Remarks on a multivariate transformation[END_REF] and Lebrun and Dutfoy [2009a].

Algorithm 1 The SMC 2 algorithm 1: Setting denition: 2: Dene the thresholds T 1 , ..., T n , the sample sizes N 1 , N 2 and the number of applications N app of Markov kernel Q.

3: Initialization: 4: Sample (θ i 0 ) 1≤i≤N 1 with probability law ν.

5: for i from 1 to N 1 do 6:

Sample ξ i,j 0 1≤j≤N 2 according to the probability density f X|θ i 0 .

7: end for 8: Transition:

9: for k from 0 to n do 10:

Associate a system of particles ξ i,j l

1≤j≤N 2 1≤l≤k
to each θ i k in order to estimate h k (θ i k ) and H k+1 (θ i k ) with Algorithm 2.

11:

Selection of the θ-particles:

12:

Sample I k = (I i k ) 1≤i≤N 1 multinomially with probability proportional to {h k (θ i k )} N 1 i=1 .
13:

Set θ i k = θ I i k k .
14:

Mutation of the θ-particles:

15:

for m from 1 to N app do 16:

for i from 1 to N 1 do 17:

Sample θ i k with a ν reversible kernel Q.

18:

Sample u with a uniform random variable.

19:

if u < 1 ∧ H k+1 ( θ i k ) H k+1 ( θ i k ) then set θ i k+1 = θ i
N 1 n = 1 N 1 N 1 i=1 δ θ i n

Estimation of collision probability between orbiting objects

The SMC 2 algorithm has been applied on the debris satellite collision test case in order to estimate π, the conditional law of Θ given φ(X) < T , with the following parameters: N 1 = 1000, N 2 = 50, N app = 1, N app2 = 1. The intermediate thresholds T i on the output distance are expressed in meters with {200, 100, 66, 50, 40, 33, 28, 25, 22, 20}. The estimators of the dierent marginals of π, obtained with the SMC 2 algorithm, are given in Figure 1, where the rst marginal is related to the rst parameter and so on.

Algorithm 2 Interacting particle system for probability estimation 1: For a given value θ i k , we build an IPS which allows to estimate both

h k (θ i k ) = P (Y < T k+1 |Y < T k , Θ = θ i k ) and H k+1 = k p=0 h p (θ i k ).
2: Setting denition: 3: Dene the number of applications N app2 of Markov kernel Q and recall the iteration parameter k and the particle value θ i k , the thresholds T 1 , ..., T k+1 and the sample size N 2 , that have been dened or obtained in Algorithm 1.

4: Initialisation: 5: Sample ξ i,j

P(Y < 20|Θ = N θ i=1 θ i m /N θ ) is equal to 0.034.
The question is how to analyze the estimated density of π for the tuning of Θ.

A possible approach is to consider the Kullback-Leibler distance between the estimated marginal density of π for the parameter Θ i and the initial marginal density of ν for parameter Θ i . If the Kullback-Leibler distance is signicant for Θ i , then one can assume that Θ i has to be nely tuned and conversely. In that case, a misestimation of Θ i will indeed tend to increase the failure probability. Table 1 summaries the dierent Kullback-Leibler distance obtained for the dierent components of Θ. The rst error component Θ 1 of the position vector seems to be the most inuent parameter on P(φ(X) < T ). On the contrary, the second error component of position and speed vector, that are Θ 2 and Θ 5 require a lower accuracy since the considered values for these parameters lead the maximum of the collision probability; an error on these parameters will thus tend to decrease the failure probability. In the same way, the density parameter Θ 7 of the drag coecient does not require also a too ne tuning in the proposed example. It may be also interesting in practice to transform the six rst components of Θ into usual orbital parameters (the semi-major axis a, eccentricity e, inclination i, argument of perigee ω, longitude of the ascending node Ω, the mean anomaly m) and then to evaluate π in that case. The estimation of the marginals of π for the dierent orbital parameters is proposed in Figure 3. The corresponding Kullback-Leibler analysis is given in Table 2. The mean anomaly is on this use case the orbital parameter that has to be most nely tuned. There is indeed a higher chance that the collision probability increases if the mean anomaly is not correctly set. Kullback-Leibler distance between marginal density π and ν for parameters Θ i .

failure probability. The proposed SMC 2 algorithm has been described in the case of a general problem where the model is a black-box system with random inputs. This algorithm has been applied with success for the analysis of collision probability between space debris and satellite. The set model parameters inuence strongly the value of the collision probability and their value has to be carefully investigated to avoid collision probability underestimation. The complete interpretation of target law π remains complicated and has to be continued. The analysis of the particles obtained by the SMC 2 algorithm with Sobol indices [START_REF] Sobol | Sensitivity estimates for non linear mathematical models[END_REF]] is a potential perspective to this work.

  Fig. 1. Estimations of the marginals of π using the SMC 2 algorithm. The red curve corresponds to the standard normal density that is the initial marginal of the dierent parameters.

Fig. 2 .

 2 Fig. 2. Estimates of P(Y < T |θ) with Θ following ν and π.

Fig. 3 .

 3 Fig. 3. Estimations of the marginals of π using the SMC 2 algorithm on the orbital parameters. The red curve corresponds to the initial density of the orbital parameters.

Table 1

 1 Component of Θ Kullback-Leibler distance with the marginal of π

	Θ 1	0.46
	Θ 2	0.13
	Θ 3	0.30
	Θ 4	0.24
	Θ 5	0.11
	Θ 6	0.25
	Θ 7	0.10

ConclusionIn this chapter, we have proposed an original methodology to analyze the inuence of parameter model that are set for the sake of simplicity, on a rare

Acknowledgments

The work presented in this chapter is part of a CNES/ONERA PhD thesis.

6: Transition:

7: for l from 0 to k -1 do 8:

Selection of the ξ particles:

9:

for j from 1 to N 2 do 10:

11:

else Sample ξ i,j l randomly and uniformly among particles which are below the threshold T l+1 .

12: end if 13: end for 14:

Mutation of the ξ particles: 15: for r from 1 to N app2 do 16:

19:

else Sample u with a uniform random variable. 20:

then set ξ i,j l+1 = ξ i,j end for 25:

if r < N app2 then set ξ i,j l = ξ i,j l+1 .

26:

end if 27: end for 28:

29: end for

The estimated probabilities P(Y < T |Θ = θ), when Θ follows ν and π are represented in Figure 2. The mean probability P(Y < T |Θ = θ) when Θ follows ν is estimated to 3.9 10 -4 . When θ = N θ i=1 θ i m /N θ , the probability


L Laurent
email: laurent@lmt.ens-cachan.fr

P.-A Boucard
email: boucard@lmt.ens-cachan.fr

B Soulier
email: soulier@lmt.ens-cachan.fr

A Dedicated Multiparametric Strategy for the Fast Construction of a Cokriging Metamodel

Keywords: metamodel, cokriging, multiparametric strategy, LATIN method, assemblies, two-levels optimization

This paper deals with the building of a gradient-based metamodel using a dedicated strategy for solving structural assemblies problems. This work is the first part of a two-levels global optimization strategy. The general objective is to reduce computation costs; here, we focus on the costs which are associated with the generation of the metamodel. Our goal is achieved through the introduction of two main elements: what we call a "multiparametric strategy" based on the LATIN method, which reduces the computation costs when the parameters vary, and the use of a cokriging metamodel taking gradients into account. Several examples illustrate the efficiency of these two elements.

Introduction

Structural optimization involves two very different domains of expertise: numerical simulation and optimization. Of the two, simulation is the only one which generates high computation costs. The main difficulty in the context of optimization is that the optimizer often requires a large number of simulations to locate the optimum of an objective function. (In contact and friction problems, these simulations are nonlinear and, thus, even more costly.) Therefore, new methods for the resolution of complex optimization problems are needed. In this context we proposed to use a two-levels model optimization.

Here, we focus on the first level. Thus, this paper deals with the aim of reducing the computational cost associated to the building of kriging metamodels. Therefore we introduced a gradient-based kriging metamodel and a MultiParametric Strategy.

In the first part of this paper, the context of two-levels model optimization is reviewed. In the second part, a specific mechanical resolution process is introduced and its performance is discussed in detail. The third part concerns the kriging class of metamodels. Finally, the last part presents the application of the coupled multiparametric-metamodel approach to a mechanical test case.

Two-level model optimization

The proposed multilevel model strategy is based on the approach developed in [START_REF] Robinson | A case for multi-level optimisation in aeronautical design[END_REF]. In our case, we consider the two-levels strategy illustrated in Figure 1. The first level consists in a metamodel defined using a limited amount of data. The zones where the optimum can be located are determined using, for example, a genetic algorithm. This information is transferred to the second level, where a precise search of the minimum is carried out. Therefore, the second level consists in an optimization process based on the full mechanical model. Two elements contribute to the computation cost: the construction of the metamodel from the responses of the full mechanical model, and the direct optimization based on the full mechanical model. Both elements are involved both levels of the two-levels model strategy. Therefore, the strategy proposed here leads to a reduction in the computation cost of both phases.

Sim ula tor

Me tam ode l Op tim ize r Sim ula tor Op tim ize r

Optimization using the metamodel Direct optimization The first step of our two-level optimization is a classical approach that is often designated as surrogate-based optimization and, overall, our work fits in the context of parametric optimization. More details on the use of metamodels in the context of parametric optimization can be found in [START_REF] Barthelemy | Approximation concepts for optimum structural design -a review[END_REF][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][START_REF] Barton | Metamodel-based simulation optimization[END_REF].

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

Many types of metamodels can be used to find an approximate solution: polynomial regression [START_REF] Box | On the Experimental Attainment of Optimum Conditions[END_REF], methods based on neural networks [START_REF] Mcculloch | A logical calculus of the ideas immanent in nervous activity[END_REF], radial basis functions [START_REF] Hardy | Multiquadric Equations of Topography and Other Irregular Surfaces[END_REF], etc... We chose to use a particular category of metamodels called kriging approximations and, more precisely, a cokriging metamodel using derivatives. These approximations are presented in Section 4.

The multiparametric strategy

This chapter introduces the dedicated strategy called multiparametric strategy for solving problems on structural assemblies. Due to a judicious reinitialisation of the LATIN algorithm -described in Sect. 3.2 -multiparametric strategy enables one to reduce significantly the computation cost. The multiparametric strategy is presented in Sect 3.3.1 and example of 2D mechanical benchmark with two, three and four design variables is used for showing the performance of the multiparametric strategy (See Sect. 3.3.3).

The context

We are considering assemblies of linear elastic structures under the assumption of small perturbations. The only nonlinearities are considered to be due to contact or friction between parts. The strategy we use to solve this type of problem, based on twenty years of development at LMT-Cachan, relies on an iterative algorithm introduced by P. Ladevève [START_REF] Ladevèze | Sur une famille d'algorithmes en mécanique des structures[END_REF][START_REF] Ladevèze | Nonlinear computational structural mechanics: new approaches and non-incremental methods of calculation[END_REF].

In the context of assemblies, this strategy remains on three main points:

• the structure being studied is divided into substructures and interfaces. In the context of assembly problem, the natural way is to considering substructures as the parts of the assembly and interfaces as the contact zones among parts. Moreover a mixed domain decomposition is considered (force and displacement are unknows at the interfaces)

• a dedicated iterative algorithm is used to solve the mechanical problem;

• the operators of the method remain constant and depend neither on the loading nor on the parameters (friction coefficient, gap) of the interface.

The resolution process

The resolution strategy, known as the LATIN (LArge Time INcrement) approach [START_REF] Ladevèze | Nonlinear computational structural mechanics: new approaches and non-incremental methods of calculation[END_REF], consists of alternative resolutions of two groups of equations. The first group contains the local equations Γ (which may be nonlinear) related to the interfaces, and the second group contains the linear equations A d related to the substructures. The iterations between the two groups are handled 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 . Friction and contact problems along with their specific 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 schematic linear representation of A d is due to the linear behavior of the substructures. 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 is the result of a classical elastic problem. Thus, an approximate solution is described entirely by s = all the interfaces

A d Γ E - E + ŝn+1/2 s s n+1 s n
(W i , F i ).
Readers can found additional details in [START_REF] Ladevèze | Sur une famille d'algorithmes en mécanique des structures[END_REF][START_REF] Ladevèze | Nonlinear computational structural mechanics: new approaches and non-incremental methods of calculation[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.

A method for reducing computation costs

A parametric optimization involves many calculations, each carried out with a different set of the parameters of the problem (design variables such as friction coefficients, preloads, gaps, etc...). This leads to the resolution of many similar problems in the sense that only the parameters of the problems vary. Therefore, it is essential to find a method to accelerate these calculations. The method we use, called the Multiparametric Strategy (MPS), was introduced in [START_REF] Boucard | A multiple solution method for non-linear structural mechanics[END_REF][START_REF] Boucard | Application of the LATIN method to the calculation of response surfaces[END_REF], then studied and applied to various sample 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 multi-analysis strategy for contact problems with friction[END_REF].

Principle of the method

The main idea of the multiparametric strategy is to take advantage of a feature of the LATIN algorithm: at each iteration, the solver yields an approximate solution over the whole loading path and at all points in the structure. Then, if a new calculation associated with other values of the parameters is requested, the algorithm can be reinitialized using a previously converged solution. Here, we assume that the only parameters which are different are the interface parameters, so space Γ alone is affected by the change.

Each new converged solution enriches a database of possible restart solutions. Figure 3 illustrates the use of the multiparametric strategy to obtain an approximate solution starting from a previous calculation: the schematic representation shows that the required number of iterations to carry out the converged solution s 1 is bigger than the required number of iterations to obtain the other converged solutions s 2 or s 3 . Moreover, in order to obtain the solution s 3 , the LATIN algorithm can be reinitialized with one of the two available converged solutions s 1 or s 2 . We observe that the benefit depends on the choice of the starting point. Here, we will use only a "closest point" strategy: for each new requested solution the LATIN algorithm is reinitialized with the converged solution associated with the set of design parameters which is closest to the one being considered. 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].

A d Γ 1 Γ 2 Γ 3 s 1 s 2 s 1 s 3
Varia tion of para mete rs With this approach, the solution converges in fewer iterations of the algorithm and, therefore, in less time.

One Reduced Order Model: the Proper Generalized Decomposition (PGD)

Each converged solution used for reinitializing the LATIN algorithm can be considered as a reduced order model (ROM) of the solution in the sense that, as we deals here with elastic assemblies, the whole solution can be only described by the boundary conditions (more details can be found in [START_REF] Boucard | A multiple solution method for non-linear structural mechanics[END_REF]). In case of problems with nonlinear material behaviour for example this idea must be extended with the use of a PGD basis as ROM.

The Proper Generalized Decomposition is initially based on the space-time separated representations as a part of the LATIN method which was designed to deal with non-linearities such as plasticity or viscoplasticity. Then many developments has been achieved during the past three decades. Readers can found details in [START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF][START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF] for the developpments of the PGD approximations and applications. Recently the PGD was used on parametric models and this approach is initiated the use of the PGD to deal with parametric problems such as parametric optimization problems, inverse identification problems or real time simulation. In order to achieve these parametric problems the classical time-space separated approximation presented by the Equation 1 can be rewritten as presented by Equation 2where k denoted a parameter of the model. The PGD can also be used as space-parameter decomposition as presented by Equation 3.

u(x, t) ≈ Q i=1 X i (x) • T i (t) (1)
u(x, t, k) ≈ Q i=1 X i (x) • T i (t) • K i (k) (2) u(x, k) ≈ Q i=1 X i (x) • K i (k) (3)
Such a decomposition allows to reduce significantly the CPU time associated to the computation of the responses of the model on requested values of the time, space and/or parameter. Examples of the use of the parametric modelling can be found in [START_REF] Lamari | Nonconventional numerical strategies in the advanced simulation of materials and processes[END_REF][START_REF] Lamari | Routes for efficient computational homogenization of nonlinear materials using the proper generalized decompositions[END_REF][START_REF] Ammar | Solving parametric complex fluids models in rheometric flows[END_REF][START_REF] Pruliere | On the deterministic solution of multidimensional parametric models using the Proper Generalized Decomposition[END_REF][START_REF] Bognet | Advanced simulation of models defined in plate geometries: 3D solutions with 2D computational complexity[END_REF]. PGD was also carried out in the context of Multiparametric Strategy [START_REF] Heyberger | Multiparametric analysis within the proper generalized decomposition framework[END_REF]. Finally the PGD was also used to build response surface models [START_REF] Halabi | Multiparametric Response Surface Construction by means of Proper Generalized Decomposition: An extension of the PARAFAC procedure[END_REF]. Actually, wether it is for optimization, the building of a reponse surface or, more generally, the building of a metamodel, it requires values of an objective function obtained for different values of the parameters of the model. The PGD enables to reduce significantly the computation cost associated to these required values.

Illustration of the performance of the method A two design variables example.

The example considered here is 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 M P a and Poisson's coefficient ν = 0.3) which are in contact with friction. Each part was represented by a single substructure discretized into 24 × 24 bilinear quadrangles. The parametric study consisted in varying the friction coefficients µ 1 and µ 2 of the two contact interfaces. The loading consisted of two stages: first, a progressive vertical pressure P 1 up to a maximum of 50M P a applied at the top of substructure Ω 3 (the preloading stage); then, a progressive horizontal load from 0 to 30M P a applied to substructure Ω 2 . In this test case, 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 4 shows the response surface of this function obtained with 18 × 18 values of the friction coefficients.

h h j Ω 1 Ω 2 Ω 3 µ 1 µ 2 P 2 (MPa)
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 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 multiparametric strategy) was estimated using the following expression: Gain = Number of calculations × CPU time of the first calculation CPU time using the multiparametric strategy (4)

We considered that the time of each calculation without the multiparametric strategy was almost constant and equal to the time of the first calculation (8.2s The results presented in Table 1 show that the solver can reuse a previously converged solution to accelerate the resolution of similar problems. The most remarkable point is the gain obtained in the evaluation of the gradients. These were calculated using a finite difference method: each gradient required the value of the response at 3 points (a sample point plus 2 points obtained by very small variations of the parameters). The results show that for the same number of evaluations the cost was less when the gradients were calculated along with the responses than when the responses alone were calculated: in our example, 16 evaluations took 42.3s whereas 16 gradients (32 evaluations) took 58s.

A four design variables example.

In this paragraph we proposed a modified version of the previous two-variables example. By modifying the substructure Ω 1 and replace it with two new substrucures, a four-variables is introduced. Figure 5 shows the geometry of the problem. The material characteristics and the loading are retained and two new interfaces are introduced: a gap between substructures Ω 1 and Ω 4 and a new interface with friction between substructures Ω 2 and Ω 4 are introduced. In this context, two new design variables are added to the two existing design variables: thus we consider three friction coefficients µ 1 , µ 2 and µ 3 (between 0 and 0.6) and one gap j 1 (between -28µm and 48µm). Although we use a first-order finite difference method to compute the gradients and the fact that n d responses (in addition to the first calculation) of the objective function are needed to obtain the gradient on one sample point, the required time to evaluate the gradients is the same than the time required to evaluate all responses. For this example we obtain a gain around 4.5 for computing the responses and the gradients. Thus, the MPS is particularly efficient in calculating the gradients, thanks to the ability of the method to reduce the computation time for sample points which are close together. This led to our decision to use a dedicated gradients metamodel.

h h/2 j 2 Ω 1 Ω 2 Ω 3 Ω 4 µ 1 µ 2 µ 3 j 1 P 2 (M P a)

The cokriging metamodel

This type of metamodel is similar to the kriging metamodel proposed by [START_REF] Krige | A statistical approach to some mine valuation and allied problems on the Witwatersrand[END_REF] and developed by [START_REF] Matheron | Traité de géostatistique appliquée[END_REF][START_REF] Matheron | Principles of geostatistics[END_REF]. The cokriging strategy stems from multivariate geostatistics [START_REF] Wackernagel | Multivariate geostatistics: an introduction with applications[END_REF]. In our case, the cokriging metamodel is built using not only the responses, but also their gradients [START_REF] Morris | Bayesian design and analysis of computer experiments: use of derivatives in surface prediction[END_REF][START_REF] Koehler | Computer experiments[END_REF][START_REF] Chung | Using gradients to construct cokriging approximation models for high-dimensional design optimization problems[END_REF].

Notations

In this section, we will use the following notations: x (i) denotes a point in the design space D (x (i) is one of the n s sample points, while x (0) is a non-sample point); x i denotes the i th coordinate of a point x; Y (x) and Y (x) denote respectively the response of the analytical function (or the response of the mechanical model) and the approximate response given by the metamodel; finally, R(x (i) , x (j) ) is a correlation function.

Principle

The model can be viewed as the sum of two components: a linear model (which represents the trend of the data) and a departure from the linear model (which represents the fluctuations around the trend).

Y (x (0) ) = µ(x (0) ) + Z(x (0) ), ∀x (0) ∈ D (5)
E(Z) = 0, Cov[Z x (i) , Z x (i) ] = 0, ∀ x (i) , x (j) ∈ D 2 (6)
where E and Cov denote the classical statistical expected value and covariance. The main idea is to consider a covariance relation between the discrete responses of the deterministic function. This covariance depends only on the distances among the samples. Thus the departure is represented by a zero-mean, second-order stationary process. (The mean and the variance are constant, with a correlation depending on the distance.) Depending on the definition of the function µ, one can build different types of kriging or cokriging metamodels (simple kriging, where µ is the average of the response at sample design 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 sample 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 ij (7)
In our case, the correlation function chosen is 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 case of a cokriging metamodel, additional covariance relations involving the different variables are introduced [START_REF] Morris | Bayesian design and analysis of computer experiments: use of derivatives in surface prediction[END_REF]. In our case, the primary variables are the evaluations of the function being studied, and the secondary variables are the gradients:

Cov ∂Z ∂x k x (i) , Z x (j) = -σ 2 ∂R ∂x k x (i) , x (j) = σ 2 c kij (8) Cov Z x (i) , ∂Z ∂x k x (j) = -σ 2 ∂R ∂x k x (i) , x (j) = σ 2 c kij (9) Cov ∂Z ∂x k x (i) , ∂Z ∂x l x (j) = -σ 2 ∂ 2 R ∂x k ∂x l x (i) , x (j) = σ 2 c kilj (10)
With this relation, one can also build the kriging or cokriging metamodel by considering the following linear predictors of the non-sample point x (0) :

Kriging Y x (0) = ns i=1 λ i x (0) Y x (i) Cokriging Y x (0) = ns i=1 λ 0i x (0) Y x (i) + n d j=1 ns i=1 λ ij x (0) ∂Y ∂x j x (i)
Thus, the determination of the Best Linear Unbiased Predictor (BLUP) leads to the λ's which, for a kriging metamodel, minimize the Mean Square Error:

MSE Y (x (0) ) = E   ns i=1 λ i x (0) Y i -Y x (0) 2   (11)
subject to the unbiasedness condition:

E ns i=1 λ i x (0) Y x (i) = E Y x (0) (12)
The construction of a cokriging metamodel follows the same process.

In the universal kriging case, the best linear predictor of Y x (0) can be written in matrix form [START_REF] Sacks | Design and analysis of computer experiments[END_REF] as:

Y x (0) = x 0 T β + c T 0 C -1 (Y s -X β) (13)
where

β = (X T C -1 X) -1 X T C -1 Y s x 0 = f 1 x (0) f 2 x (0) • • • f k x (0) c 0 = c 01 c 02 • • • c 0ns (X) ij = f j x (j) (C) ij = c ij
The regressor vector β is the generalized least-squares approximation of β. Thus, the first part of Equation 13 is the generalized least-squares prediction at point x (0) . The second part can be viewed as a correction of the generalized least-squares response surface to obtain the interpolating kriging model.

The formulation of the cokriging metamodel is very similar to that of the kriging metamodel. For simplicity's sake, following the idea in [START_REF] Chung | Using gradients to construct cokriging approximation models for high-dimensional design optimization problems[END_REF], we propose to use only an ordinary cokriging metamodel:

Y x (0) = βc + c T c0 C -1 c (Y sc -X c βc ) (14)
where

βc =(X T c C -1 c X c ) -1 X T c C -1 c Y sc c c0 = c 01 c 02 • • • c 0ns c 110 c 210 • • • c n d ns0 X c = 1, 1, • • • , 1, 0, 0, • • • , 0 T C c = C C d C d T C dd
and where

(C) ij = c ij (C d ) ij =        c 111 c 211 • • • c n d 11 c 121 • • • c n d ns c 112 c 212 • • • c n d 12 c 122 • • • c n d ns2 c 113 c 213 • • • . . . . . . . . . c 11ns • • • c n d nsns        (C dd ) ij =            c 1111 c 1121 • • • c 11n d 1 c 1112 • • • c 11n d ns c 2111 c 2121 • • • . . . . . . . . . c n d 1n d 1 c n d 112 • • • . . . c n d ns11 • • • c n d nsn d ns           
In the kriging case, vector Y s contains only the responses of the function at the sample point, whereas in the cokriging case vector Y sc contains both the responses and the gradients.

With these formulations, the models can supply approximate responses of the actual function at all the points in the design space. In our case, we consider the response of the function to be deterministic and we obtain both of the interpolating models. This type of metamodel has other advantages: for example, it provides statistical information (the expected value and the variance of the process). 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 -x 0 c 0 0 X T X C -1 x 0 c 0 (15)
This relation also holds in the cases of simple or ordinary kriging or cokriging models, provided the appropriate forms of vectors and matrices x 0 , c 0 , X and C are used.

Estimation of the parameters

The model's parameters (such as the characteristic correlation length scale l, the variance σ of the random process Z or the regression coefficients) can be determined by maximizing the likelihood [START_REF] Mardia | Maximum likelihood estimation of models for residual covariance in spatial regression[END_REF]. In our case, we use this technique to determine l and σ:

(l, σ 2 ) = arg min l,σ 2 - N 2 log(2πσ 2 ) - 1 2 log(|C(l)|) - 1 2σ 2 (Y s -X β) T C(l) -1 (Y s -X β) (16)
with

N = n s for kriging n s (n d + 1) for cokriging
The variance σ 2 can be determined analytically through the derivation of the likelihood function:

σ2 = 1 N (Y s -X β) T C(l) -1 (Y s -X β) (17)
One can also use an optimizer to determine the value of the correlation length scale. This method has some drawbacks [START_REF] Warnes | Problems with likelihood estimation of covariance functions of spatial Gaussian processes[END_REF][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]: in many cases with very few points, the log-likelihood is monotonous; often, the correlation matrix is affected by conditioning problems which make finding the minimum difficult. When such problems arise, one sets the parameters (in particular the correlation length scale) to fixed values.

Analytical applications

In this section, we apply kriging and cokriging to one-and two-dimensional analytical test functions. We will use the abbreviations OK and OCK to designate respectively Ordinary Kriging and Ordinary CoKriging.

The one-dimensional test function

First, we applied the two types of metamodels to an analytical function, chosen to be y(x) = exp(-x/10) cos(x) + x/10. We used 5 sample responses of the analytical function to build the OK metamodel, and an additional 5 sample derivatives to build the OCK metamodel. The correlation function was the Matérn function. The sample points were obtained using Latin Hypercube Sampling. Figure 7 illustrates the capability of the cokriging metamodel to interpolate not only the values of the responses, but also the sample derivatives. In this example, the quality of the cokriging metamodel was better. This statement remained true as long as we worked with only a few points. We can also observe that for a relatively smooth function the kriging metamodel converges quickly toward a good approximation of the real function when the number of points becomes large enough.

The characteristics of the two metamodels are summarized in Table 2.

R 2 RAAE RMAE Q 1 Q 2 Q 3 OK 0.
Q 1 = sup i∈{1,2,...,nc} e i Q 2 = nc i=1 e i Q 3 = Q 2 n c
where ∀i ∈ {1, 2, ..., n c },

e i = (Y x (i) -Y x (i) ) 2 sup j∈{1,2,...,nc} Y x (j) 2
Based on the statistical information from the kriging and cokriging metamodels, one can derive confidence intervals. The two diagrams of Figure 8 show the 95% Confidence Intervals (CIs) obtained with Expression 18. In these types of metamodels, the size of the confidence envelopes is determined mainly by the distance between each pair of neighboring points.

CI ± (x) = Y (x) ± 2 M SE Y (x) (18)

The 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 bumpy function: the six-hump camel back function

(∀(x 1 , x 2 ) ∈ [-3, 3] × [-2, 2], 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 ). Figure 9 shows the response surface of this function. The two metamodels were constructed using 16 evaluations of the function for the kriging metamodel, and an additional 16 evaluations of the gradients for the cokriging metamodel. In both cases the correlation function was the Matérn function. The characteristics of these two metamodels are given in Table 3. Table 3: Characteristics of the two previous metamodels (Figure [START_REF] Ladevèze | Nonlinear computational structural mechanics: new approaches and non-incremental methods of calculation[END_REF][START_REF] Blanzé | A modular approach to structure assembly computations: application to contact problems[END_REF] 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. But for a problem involving the evaluation of mechanical responses the computation cost for the determination of n s responses and n s gradients is, of course, higher than that required to obtain n s responses alone. (In the former case, due to the use of finite differences to obtain gradients, the construction of the metamodel requires 3n s evaluations of the mechanical model.)

R 2 RAAE RMAE Q 1 Q 2 Q 3 OK 0.
Now let us take another approach: the idea is to build metamodels using the same number of evaluations. In the following case, we used the same six-hump function, but the kriging metamodel was constructed based on the responses at 27 sample points (Figure 12), while the cokriging metamodel was still constructed using the responses and the gradients at 9 sample points (Figure 13). [START_REF] Champaney | Modular analysis of assemblages of three-dimensional structures with unilateral contact conditions[END_REF][START_REF] Boucard | A multiple solution method for non-linear structural mechanics[END_REF] Table 4 shows that the quality of the prediction achieved with 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, thanks to the multiparametric strategy and as shown in Table 1, would be much lower than that associated with the 27 points.

R 2 RAAE RMAE Q 1 Q 2 Q 3 OK 0.

The five-and ten-dimensional test functions

In this section, we consider two analytical functions defined for n d variables:

• the Rosenbrock function, which is a smooth function, can be defined as follow:

∀x ∈ [-2, 2] n d , f (x) = n d -1 i=1 [100(x 2 i -x i+1 ) 2 + (x i -1) 2 ] (19)
• the Schwefel function, which is a very irregular function, can be defined as follow:

∀x ∈ [-500, 500] n d , f (x) = 418.9829 × n d - n d i=1 [x i sin( |x i |)] (20)
Here, we only use these analytical function with 5 or 10 variables (i.e. n d = 5 or n d = 10). Using a Latin Hypercube Sampling 10 (for n d = 5) and 2 (for n d = 10) sets of sample points are generated for each number of sample points. Here we consider the number of sample points in the horizontal axis, thus for n s sample points the kriging metamodel is built from n s responses and the cokriging metamodels is built from n s responses and n s gradients (each gradient comprises n d components). The mean of the criteria R 2 and Q 3 is calculated for each number of sample points: metamodels are built with responses (and gradients) for each set of sample points obtained for each number of sample points; then for each metamodel the quality criteria are computed and finally for each number of sample points we compute the means of the two criteria. Figures 14 and15 show the means of the quality criteria versus the number of sample points for the kriging and cokriging metamodels built using responses and gradients provided by the Rosenbrock and Schwefel functions with 5 variables. Figures 16 and17 show the same results obtained with the ten-dimensional functions. Moreover, if we consider the case of the five-dimensional Schwefel function with the same total number of evaluations (i.e., if we consider n s responses and 5 × n s gradients for the cokringing, we can use n s + 5 × n s responses for the kriging), the quality is always better with the cokriging metamodel as soon as n s is greater than 30. As we have shown before that the MPS is particularly efficient in calculating the gradients, coupling a cokriging metamodel with the MPS should lead to a very efficient strategy to build metamodels. This is what we illustrate in the next part.

Application to a mechanical problem

In this section, we present the construction of the two previous metamodels (kriging and cokriging) using the multiparametric strategy presented in the first part of this paper and the actual responses of a mechanical model. The test case being considered is the three-squares example presented in Page 5. We carried out two studies: one with a fixed number of mechanical calls, and the other with a fixed quality of the metamodels.

Case of a fixed number of mechanical calls

Both metamodels were constructed using 15 mechanical evaluations. The kriging metamodel was determined using 15 values of the force, and the cokriging metamodel was determined with 5 values of the force and 5 gradients. The sample points were obtained through Latin Hypercube Sampling. Our calculations led to the following two response surfaces (Figures 18 and19) to be compared with the actual response surface of Figure 4: Tables 5 and6 show the characteristics of the calculations and of the metamodels. In this case, 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.

Gradients

Case of a fixed metamodel quality

Another way to study the cost of the metamodels consists in constructing kriging and cokriging metamodels with the same quality. For this study, 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 7.

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. In the case of more complex problems and more parameters, one can expect even better gains.

OK

Conclusion

In this paper, we associated a cokriging metamodel with a dedicated multiparametric strategy. These two elements enabled us to achieve a significant reduction in the computation cost of constructing the metamodel.

Our multiparametric strategy takes advantage of a property of the LATIN method which is that the result of a previous calculation constitutes an efficient starting point for a subsequent calculation. This property is particularly advantageous for the evaluation of gradients, which are obtained using finite differences (i.e. with very small variations of the parameters). Thus, our strategy lends itself naturally to the use of a gradient-based metamodel such as the cokriging model presented in the paper.

Moreover, since this metamodel is part of a multilevel model optimization strategy, besides reducing the computation cost associated with the construction of the cokriging metamodel, one can a lso reduce the cost of the second level, in which algorithms based on gradient descent are used. This will be the subject of an upcoming work.

Figure 1 :

 1 Figure 1: The multilevel optimization

Figure 2 :

 2 Figure 2: Schematic representation of the LATIN process

Figure 3 :

 3 Figure 3: The multiparametric strategy using the LATIN algorithm

2 FFigure 4 :

 24 Figure 4: The geometry of the problem and the response surface

Figure 5 :

 5 Figure 5: The geometry of the problem of the four-variables example

  Total (w/o grad.) Total (w/-grad.) Responses Gradients 5Number of sample pointsGain of the MPS

Figure 6 :

 6 Figure 6: CPU time and gain of the MPS associated two different sets of sample points for the four-variables test-case

Figure 7 :

 7 Figure 7: The ordinary kriging and cokriging metamodels and their derivatives

Figure 8 :

 8 Figure 8: The confidence envelopes for the kriging and cokriging metamodels

Figure 9 :

 9 Figure 9: The actual function and its gradients (blue arrows)

Figure 10 :Figure 11 :

 1011 Figure 10: The kriging function and its gradients (blue arrows)

Figure 12 :Figure 13 :

 1213 Figure 12: The kriging surface obtained with 27 responses (27 evaluations of the function)

Q 3 Figure 14 :

 314 Figure 14: Quality of kriging and cokriging metamodels applied on five-dimensional Rosenbrock function

Figure 15 :Q 3 Figure 16 :Q 3 Figure 17 :

 15316317 Figure 15: Quality of kriging and cokriging metamodels applied on five-dimensional Schwefel function

Figure 18 :Figure 19 :

 1819 Figure 18: The kriging response surface of the three-squares test case (using 15 responses)

Table 1 :

 1 ). The mechanical calculations using the multiparametric strategy

		Eval. Eval. gradients	Total
	Number of evaluations	16	2 × 16	48
	Total CPU time	42.3s	58.1s	100.4s
	Avg CPU time per calc.	2.7s	1.8s	2.1s
	Min CPU time	0.2s	0.2s	0.2s
	Max CPU time	8.2s	6s	8.2s
	Nb. of iterations	1,896	2,580	4,476
	Avg nb. of iterations (per calc.) 118.5	80.6	93.3
	Min nb. of iterations	5	5	5
	Max nb. of iterations	381	260	381
	Gain	3.11	4.54	3.94

Table 2 :

 2 Characteristics of the two previous metamodelsCriteria RAAE and RMAE stand for Relative Average Absolute Error and Relative Maximum Absolute Error. Criteria Q i , which compare the actual response and the responses of the metamodels at n c points with n c >> n s , were calculated as follows:

  OCK 0.9958 3.027 • 10 -2 5.330 • 10 -1 1.391 • 10 -2 0.1824 2.027• 10 -4

	6991	0.4042	2.845	0.3946	13.192 1.466 • 10 -2

  8175 0.2566 2.933 0.4195 8.001 8.890 • 10 -3 OCK 0.8160 0.2580 2.541 0.3150 8.066 8.9627 • 10 -3

Table 4 :

 4 Characteristics of the two previous metamodels (Figure

Table 5 :

 5 Characteristics of the mechanical calculations OK 0.9865 8.751 • 10 -2 0.6624 2.746 • 10 -2 0.3385 8.463 • 10 -4 OCK 0.9822 0.1015 0.4024 1.014 • 10 -2 0.4441 1.110 • 10 -3

	R 2	RAAE	RMAE	Q 1	Q 2	Q 3

Table 6 :

 6 Characteristics of the metamodels (OK: ordinary kriging/OCK: ordinary cokriging)

  10 -5 1.860 • 10 -3 3.626 • 10 -4

			OCK	OCK
		9 × 9	3 × 3	5 × 5
	Total CPU Time	154.7s	61.6s	127s
	Avg CPU Time (/calc)	1.9s	2.3s	1.7s
	Min CPU Time	0.2s	0.2s	0.2s
	Max CPU Time	8.9s	8.1s	9.1s
	LATIN iterations	6,576	2,731	5,436
	Avg iterations (/calc)	81.2	101	72.5
	Min iterations	5	5	5
	Max iterations	381	381	381
	Gain	4.65	3.53	5.39
	R 2	0.9878	0.9702	0.9942
	RAAE	2.511 • 10 -2	0.1234	5.147 • 10 -2
	RMAE	0.1442	0.5703	0.3243
	Q 1	1.301 • 10 -3 2.036 • 10 -2 6.585 • 10 -3
	Q 2	3.044 • 10 -2	0.7439	0.1450
	Q 3	7.610 •

Table 7 :

 7 The metamodels constructed using full factorial sampling

Acknowledgment

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