Janis Janusevskis
email: janusevskis@emse.fr

Rodolphe Le Riche
email: leriche@emse.fr

Simultaneous kriging-based sampling for optimization and uncertainty propagation

Robust analysis and optimization is typically based on repeated calls to a deterministic simulator that aim at propagating uncertainties and finding optimal design variables. Without loss of generality a double set of simulation parameters can be assumed: x are deterministic optimization variables, u are random parameters of known probability density function and f (x, u) is the objective function attached to the simulator. Most robust optimization methods involve two imbricated tasks, the u's uncertainty propagation (e.g., Monte Carlo simulations, reliability index calculation) which is recurcively performed inside optimization iterations on the x's. In practice, f is often calculated through a computationally expensive software. This makes the computational cost one of the principal obstacle to optimization in the presence of uncertainties.

This report proposes a new efficient method for minimizing the mean objective function, minx∈I R n EU f (x, U ). The efficiency stems from the simultaneous sampling of f for uncertainty propagation and optimization, i.e., the hierarchical imbrication is avoided. Y t (x,u) (ω), a kriging (Gaussian process conditioned on t past calculations of f ) model of f (x, u) is built and the mean process,

is analytically derived from it. The sampling criterion that yields both x and u is the one-step ahead minimum variance of the mean process Z at the maximizer of the expected improvement. The method is compared with Monte Carlo and kriging-based approaches on analytical test functions in two, four and six dimensions.

* to account for manufacturing uncertainties or errors that are introduced during the manufacturing process.

.

* in the field of Multidisciplinary Design Optimization, where the disciplines are optimized without full information about the outputs of other disciplines (which are inputs of the subsystem to be optimized).

Introduction

Modern design analysis and optimization tasks are based on deterministic simulation models. In cases where these models involve time-consuming evaluations the methods based on Response Surface models (metamodels, surrogate models) are often used. A popular method for creating such metamodels is kriging or Gaussian Process (GP) regression [START_REF] Sacks | Design and analysis of computer experiments[END_REF], [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning[END_REF], which provides probabilistic estimates of the underlying model. These estimates can be used in optimization procedures to formulate the so called Infill Sampling Criteria (ISC) [START_REF] James | Flexibility and Efficiency Enhancements for Constrained Global Design Optimization with Kriging Approximations[END_REF] that allow efficiently sample for points.

However, real life problems are almost always non-deterministic. The procedures of robust design attempt to incorporate and account for some of the uncertainties and find designs whose quality degrades in a controlled, gradual way due to perturbations in parameters.

Robust design optimization is of interest:

* to account for model errors. Conventional deterministic optimization depends on the objective function and the constraints. However, these functions are approximations of the real world. As long as one does not have a detailed knowledge of the error function of the model, one cannot be certain that the model optimum can be mapped to the true optimum [START_REF] Beyer | Robust optimization -a comprehensive survey[END_REF].

Extensive survey on robust optimization can be found in [START_REF] Park | Robust design : An overview[END_REF], [START_REF] Beyer | Robust optimization -a comprehensive survey[END_REF]. Different formulations for robust design optimization may be found in [START_REF] Beyer | Robust optimization -a comprehensive survey[END_REF], [START_REF] Jin | The use of metamodeling techniques for optimization under uncertainty[END_REF], [START_REF] Salazar | An empirical study of the use of confidence levels in RBDO with Monte Carlo simulations[END_REF], [START_REF] Du | Towards a better understanding of modeling feasibility robustness in engineering design[END_REF], [START_REF] Taflanidis | Stochastic system design and applications to stochastically robust structural control[END_REF], [START_REF] Nikolaos | Optimization under uncertainty: State-of-the-art and opportunities[END_REF]. Often the robust design optimization is expressed as the simultaneous minimizer of the mean and variation of the objective function while satisfying probabilistic constraints. As stated in [START_REF] Beyer | Robust optimization -a comprehensive survey[END_REF] "modeling robustness aspects is only one side of the coin, solving the related optimization problems is the other often computationally demanding side, and only a few computational tools have been developed so far" [START_REF] Beiqing | A robust design method using variable transformation and gauss-hermite integration[END_REF].

In the simplified strategies the optimization problem is reformulated so as to be solved using standard techniques of mathematical programming. These techniques of robust optimization are restricted to cases where the original problem can be well approximated by linear or quadratic functions [START_REF] Beyer | Robust optimization -a comprehensive survey[END_REF]. Keeping aside special engineering problems where the robustness measures can be obtained explicitly and thus leading to a classical (deterministic) optimization problem, usually some sort of estimates of the robustness measures are used [START_REF] Du | Towards a better understanding of modeling feasibility robustness in engineering design[END_REF], [START_REF] Beiqing | A robust design method using variable transformation and gauss-hermite integration[END_REF]. The simple approach involves Monte Carlo (MC) simulations, thus leading to optimization of noisy functions. Given a fixed design point x, one estimates the statistics (mean value, variance, etc.) and use them as inputs for a (derivative free) numerical optimization algorithm. The main drawback of such approach is the number of points, k, necessary to obtain accurate MC estimates. Therefore a trade off must be made between the computational cost (in terms of number of runs k) and the precision (fidelity level) of the solution. To reduce the cost of repeated MC evaluations during the optimization, the number of runs k may be determined, using, for example, a statistical test [START_REF] Cantú-Paz | Adaptive sampling for noisy problems[END_REF], [START_REF] Salazar | An empirical study of the use of confidence levels in RBDO with Monte Carlo simulations[END_REF]. A detailed presentation of previous works where the noisy optimization is based on kriging metamodels is given in section 1.1. Different approaches exist that try to break the nested relationship between MC sampling and optimization loops. When using MC simulations to obtain robustness measures, the metamodels are used in order to reduce the time associated with sampling and also to smooth and filter the noise of the obtained statistics. The surrogate models may be utilized at the level of the optimization, of the uncertainty quantification or at both levels. A review and different formulations of surrogate based approaches and case studies involving trust regions (in order to manage the accuracy of the approximations) can be found in [START_REF] Wojkiewicz | Formulations for surrogate-based optimization under uncertainty[END_REF]. Evolutionary Algorithms (EA) [START_REF] Eiben | Introduction to Evolutionary Computing (Natural Computing Series)[END_REF] are intrinsically robust to the noise in the objective function. A literature review addressing uncertainties in evolutionary optimization is given in [START_REF] Yaochu | Evolutionary optimization in uncertain environments-a survey[END_REF]. The simplest approach is based on multiple evaluations of the objective function at every population member and the mean value is utilized as fitness (direct MC implementation). In other approaches the noisy value is used directly as an input for optimization algorithm, see for example [START_REF] Beyer | Towards a theory of 'evolution strategies'. Some asymptotical results from the (1, +λ)-theory[END_REF], [START_REF] Beyer | Functions with noise-induced multimodality: a test for evolutionary robust optimization-properties and performance analysis[END_REF]. As stated in [START_REF] Beyer | Robust optimization -a comprehensive survey[END_REF] there are reasons to believe that EA are especially good at direct robust optimization (i.e., without calculating robust measures explicitly). However both approaches increase the number of function evaluations per generation by a factor between three and 100 [START_REF] Hansen | A method for handling uncertainty in evolutionary optimization with an application to feedback control of combustion[END_REF].

We consider problems where the true function is given by a complex deterministic simulator f (x, u) such as a CFD solver, which gives a constant value when x and u are fixed. In the context of robust optimization we will look at u as an instance of random variables U , therefore function f (x, U ) becomes random itself. In such cases where the optimization is based on computation intensive simulations (even for deterministic optimization problems) stochastic search algorithms are not suitable. Assuming that gradient information and analytical expressions of the robustness measures are not available, we focus on the methods that are more suited for such problems and in particular, kriging based optimization algorithms. A new methods is proposed that simultaneously minimizes the expectation function and optimally samples the random variables.

A bibliography of kriging based optimization of noisy simulators

Kriging based optimization procedures have become popular for the optimization of time consuming simulators. The first algorithm was introduced by by Schonlau, Welch and Jones [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF], [START_REF] Donald | A taxonomy of global optimization methods based on response surfaces[END_REF] and named Efficient Global Optimization (EGO). EGO uses an ISC called Expected Improvement (EI) to provide the next point for evaluation by the deterministic simulator [START_REF] James | Flexibility and Efficiency Enhancements for Constrained Global Design Optimization with Kriging Approximations[END_REF]. In recent years modifications to this approach have been proposed thus adapting it for the optimization of nosy functions and partially converged simulations. Forrester and others [START_REF] Forrester | Optimization using surrogate models and partially converged computational fluid dynamics simulations[END_REF] have used kriging optimization for the multi-fidelity models based on CFD simulations. The overall optimization is divided into two steps. In the first step the initial Design of Experiments (DOE) is run using low fidelity calculations, the metamodels are created and points are added until a criterion of approximation quality is satisfied. In the second step, additional high fidelity (fully converged) points are evaluated, using the error of approximation of the low and high fidelity solutions, the EI based optimization is performed.

In [START_REF] Huang | Global optimization of stochastic black-box systems via sequential kriging meta-models[END_REF], Huang and others extend the EGO algorithm for stochastic black box optimization tasks (Sequential Kriging Optimization). The next point in the optimization is chosen by maximizing the Augmented Expected Improvement (AEI), where the kriging is used to filter the noise at the current best point. The EI is augmented by an utility function, which accounts for the uncertainty of the model, and multiplied by a factor, which allows to account for the diminishing improvement of additional replicates as the prediction becomes more accurate. The algorithm is compared with alternative methods (Revised Simplex Search -RSS-, Simultaneous Perturbation Stochastic Approximation -SPSA-, DIRECT).

In [START_REF] Huang | Sequential kriging optimization using multiple-fidelity evaluations[END_REF], the authors use multiple metamodels to describe the systemic error of the corresponding fidelity level (for modelling multiple fidelity models, see [START_REF] Kennedy | Predicting the output from a complex computer code when fast approximations are available[END_REF]). Authors extend the AEI function by connecting it to the evaluation costs. The extended AEI is then employed to select both the location and the fidelity level for the next Finite Element (FE) simulation. Vazquez et al. [START_REF] Villemonteix | Global optimization of expensiveto-evaluate functions: an empirical comparison of two sampling criteria[END_REF] empirically compared several EI criteria (using best observations, minimum of kriging mean, AEI and Informational Approach to Global Optimization [START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF]) for the optimization of noisy functions. The tests were generated by GP realizations. The kriging is used to filter white noise.

In [START_REF] Picheny | Noisy expected improvement and on-line computation time allocation for the optimization of simulators with tunable fidelity[END_REF], a kriging based method accounts for uncertainty (fidelity level) at every data point, filters the noise and selects next data point. The approach is based on modified EI that describes quantiles. The next point is chosen by maximizing the modified EI. The proposed improvement function is also used to estimate the needed fidelity level. Simulator time slots are added iteratively to increase fidelity level and decrease noise at the next point. The benefit of the new time slot is measured by the EI decrease. The algorithm takes into account limited resources.

In [START_REF] Garnett | Gaussian processes for global optimization[END_REF], a multi-step lookahead optimization procedure is proposed using Bayesian Monte Carlo integration to obtain marginalized hyper-parameters of the process model [START_REF] Michael A Osborne | Towards real-time information processing of sensor network data using computationally efficient multi-output gaussian processes[END_REF]. The Bayesian expected loss criterion is introduced which replaces EI. The next point minimizes the expected loss. The proposed method is extended to the noisy case where the objective is corrupted by Gaussian noise. The one-step and two-steps ahead techniques are compared for nosy and non noisy test beds with other methods, such as SNOBFIT [START_REF] Huyer | SNOBFIT -stable noisy optimization by branch and fit[END_REF] and Implicit Filtering (IF) for noisy functions.

Conditioning problems are avoided by requiring that observations are sufficiently far from each other. When one observation is very close to an existing one, instead of directly adding the point, the authors add derivative information on the line that connects the two points.

Another approach has been introduced in the Machine Learning community, where the GP is constructed stimulating linear filters with white noise. Such procedure is equivalent to constructing GP using kernel convolution with white noise. It allows to obtain covariance kernels without directly specifying and parametrizing new covariance kernels. For example, it can model multiple outputs and their dependency, and be applied for models of different fidelity levels [START_REF] Boyle | Dependent gaussian processes[END_REF][START_REF] Boyle | Gaussian processes for regression and optimisation[END_REF].

Previously discussed methods deal with optimization problem of the form min

x f (x) + ǫ x ∈ R m
where ǫ introduces some noise. More generally, robust optimization addresses the problem of min

x f (x, u),
where the source of the noise in the objective is due to uncertainties in some of the model variables and/or parameters u.

Note that only some of the most sophisticated approaches (e.g. [START_REF] Huang | Global optimization of stochastic black-box systems via sequential kriging meta-models[END_REF] and [START_REF] Picheny | Noisy expected improvement and on-line computation time allocation for the optimization of simulators with tunable fidelity[END_REF]) provide means for choosing the precision (number of MC runs) of the simulations. And most of the methods dealing with noisy functions assume normally distributed and independent noise sources. In this report, a method is presented that solves the robust optimization problem min x E U [f (x, U )]. A kriging model of the f (x, u) is created and analytically translated into a kriging model of E U [f (x, U )], called the projected model. An EI and a sampling criteria define next (x, u) points for simulations. An important feature of the method is therefore that the optimization and the sampling are handled simultaneously.

GP for noisy inputs and convolutions

There is a limited amount of papers that are devoted to building GP while accounting for uncertain inputs.

In the work of Girard [START_REF] Girard | Approximate methods for propagation of uncertainty with gaussian process models[END_REF] the uncertainty on the inputs is propagated to the predictive (conditional) mean and covariance from a Bayesian perspective. The predictive mean and variance at uncertain points and also the predictive mean and covariance of conditional Gaussian Process for uncertain inputs are obtained in the case of Gaussian kernels and normally and independently distributed noise on the input parameters. In the case of non Gaussian kernels the author suggest to use Taylor series approximations. The results are applied for multi step ahead time series analysis and control of non-linear dynamic systems [START_REF] Girard | Gaussian process priors with uncertain inputs -application to multiple-step ahead time series forecasting[END_REF], [START_REF] Quiñonero-Candela | Propagation of uncertainty in bayesian kernel models -application to multiple-step ahead forecasting[END_REF]. Some comments on the construction of covariance kernels in case of noisy input can be found in [START_REF] Seeger | Bayesian Gaussian Process Models: PAC-Bayesian Generalisation Error Bounds and Sparse Approximations[END_REF] (section 2.1.8.2). It is noted that when the inputs contain independent and normally distributed noise, it is possible to calculate the mean and variance of the resulting process. However in [START_REF] Girard | Approximate methods for propagation of uncertainty with gaussian process models[END_REF] and [START_REF] Seeger | Bayesian Gaussian Process Models: PAC-Bayesian Generalisation Error Bounds and Sparse Approximations[END_REF] authors wrongly state that the resulting process in not Gaussian.

In [START_REF] Higdon | Space and space-time modeling using process convolutions[END_REF] a summary of possible applications of kernel convolutions is given. Some of those are non-parametric covariance modelling, dependent Gaussian processes, non-stationary covariance construction, reduction in dimension, etc.

In the work of Paciorek [START_REF] Paciorek | Nonstationary Gaussian Processes for Regression and Spatial Modelling[END_REF], a generalization of the kernel convolution covariance of Higdon et al. (1999) is used to construct new non-stationary kernels. The smoothing properties of the new non-stationary kernels depend on the smoothness of the underlying stationary kernels. A particular attention is given to the non-stationary version of Matèrn correlation function and its smoothness properties.

Problem formulation

We minimize mean response of the simulator

min x E U [f (x, U )] (1)
where x ∈ R n and u ∈ R m is an instance of multivariate random variable U with joint probability distribution and E U is the expectation measure. We assume that each evaluation of f (x, u) ∈ R involves time consuming call to a deterministic simulator. We say that X = {x 1 , ..., x t } and U = {u 1 , ..., u t }, where t ∈ N is the number of initial design points (number of calls to the simulator f (x, u)) and Y = {f (x 1 , u 1 ), .., f (x t , u t )} is the simulator responses at the design points. By Y (x,u) (ω) we denote Gaussian process with given mean m Y (x, u) and given covariance function

cov Y (x, u; x ′ , u ′ ) E[Y (x,u) (ω)] = m Y (x, u), COV(Y (x,u) (ω), Y (x ′ ,u ′ ) (ω)) = cov Y (x, u; x ′ , u ′ ), (2)
Y (x,u) (ω) ∼ GP (m Y (x, u), cov Y (x, u; x ′ , u ′ )) . (3)
We assume that f (x, u) can be approximated by y(x, u) which is an instance of Y (x,u) (ω) and f (x, U ) can be represented by Y (x,U ) (ω). Let Y t (x,u) (ω) denote a GP conditioned on the t observations

Y t (x,u) (ω) = Y (x,u) (ω)|Y (X,U) (ω) = Y . (4)
The conditional process Y t (x,u) (ω) utilizes data of the observations. The conditional expectation and the conditional covariance [START_REF] Sacks | Design and analysis of computer experiments[END_REF], [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning[END_REF] are

E Y t (x,u) (ω) = m Y (x, u) + cov Y (x, u; X, U)cov Y (X, U; X, U) -1 (Y -m Y (X, U)), (5)
COV Y t (x,u) (ω) = cov Y (x, u; x, u) -cov Y (x, u; X, U)cov Y (X, U; X, U) -1 cov Y (X, U; x, u). (6)
We can therefore estimate the mean (predicted kriging mean) and the variance (predicted kriging variance) of Y t (x,u) (ω) at any unobserved point (x, u).

Projected process

In order to minimize the expectation eq. ( 1) and take full advantage of Gaussian process metamodel eq. ( 4), it is possible to construct a new Gaussian process by taking the expectation of the conditioned process with respect to

U Z t (x) (ω) = E U [Y t (x,U ) (ω)] = R m Y t (x,u) (ω)d µ (u), (7)
where d µ (u) is a probability measure on u.

From eq. ( 7) we can see that the Z t (x) (ω) is a linear combination of Gaussian processes and therefore also a Gaussian process and we call it "projected process". The link can be made with Bayesian Monte Carlo integration [START_REF] Edward | Bayesian monte carlo[END_REF], where we have Gaussian prior on the function y t (x,u) .

Expectation

To obtain the mean of Z (x) (ω) as a function of x we calculate

m Z (x) =E Y Z (x) (ω) = E Y E U [Y (x,U ) (ω)] = Ω R m y(x, u)dµ(u)dP (ω) (8)
If we assume that Fubinis condition is satisfied

Ω R m |y(x, u)|dµ(u)dP (ω) < ∞, (9)
we can rewrite eq. ( 8) as

m Z (x) = R m Ω y(x, u)dP (ω)dµ(u) = R m E Y [Y (x,u) (ω)]dµ(u) (10)
in other terms, under Fubinis condition,

E Y Z (x) (ω) = E Y E U [Y (x,U ) (ω)] = E U E Y [Y (x,U ) (ω)] . (11)
Note that, because Y t is conditioned on points with finite values, Fubini's condition seems satisfied (we do not provide a rigorous proof). We see that the mean of Z (x) (ω) is an integral over the u space of the mean of Y (x,u) (ω)

m Z (x) = R m m Y (x, u)dµ(u). (12)
Z t (x) (ω) in eq. ( 7) is defined in terms of the conditional process Y t (x,u) (ω), therefore, the mean m Y (.) is the conditional mean in eq. ( 5). If

X * = [x 1 * , ...x * r ] ⊤ is a r × n matrix of r test points we have E Z t (X * ) (ω) = R m m Y (X * , u)dµ(u) + R m cov Y (X * , u; X, U)cov Y (X, U; X, U) -1 (Y -m Y (X, U))dµ(u) (13)
Here

cov Y (X, U; X, U) -1 is a t × t matrix and (Y -m Y (X, U)) is t × 1 vector. If we denote by c the t × 1 vector cov Y (X, U; X, U) -1 (Y -m Y (X, U)
), then eq. ( 13) becomes

E Z t (X * ) (ω) = R m m Y (X * , u)dµ(u) + t i=1 c i R m cov Y (X * , u; x i , u i )dµ(u) = R m m Y (x j * , u)dµ(u) 1≤j≤r + t i=1 c i R m cov Y (x j * , u; x i , u i )dµ(u) 1≤j≤r . (14)

Covariance

To obtain the covariance of Z (x) (ω) as a function of x and x ′ we use a similar reasoning

cov Z (x, x ′ ) = COV(Z (x) (ω), Z (x ′ ) (ω)) = E Y Z (x) (ω) -E Y Z (x) (ω ′ ) Z (x ′ ) (ω) -E Y Z (x ′ ) (ω ′′ ) = E Y E U [Y (x,U ) (ω)] -E Y E U [Y (x,U ) (ω ′ )] E U [Y (x ′ ,U ) (ω)] -E Y E U [Y (x ′ ,U ) (ω ′′ )] (15)
Using the same transformation as in eq. ( 10) and assuming that eq. ( 9) is satisfied, we can rewrite eq. ( 15) as

cov Z (x, x ′ ) = E Y E U [Y (x,U ) (ω)] -E U E Y [Y (x,U ) (ω ′ )] E U [Y (x ′ ,U ) (ω)] -E U E Y [Y (x ′ ,U ) (ω ′′ )] (16)
Using the linearity of the expectation we write

cov Z (x, x ′ ) = E Y E U Y (x,U ) (ω) -E Y [Y (x,U ) (ω ′ )] E U Y (x ′ ,U ) (ω) -E Y [Y (x ′ ,U ) (ω ′′ )] (17)
By the definition of expectation

cov Z (x, x ′ ) = Ω R m Y (x,u) (ω) - Ω Y (x,u) (ω ′ )dP (ω ′ ) dµ(u) R m Y (x ′ ,u ′ ) (ω) - Ω Y (x ′ ,u ′ ) (ω ′′ )dP (ω ′′ ) dµ(u ′ ) dP (ω) (18)
rewriting we obtain

cov Z (x, x ′ ) = Ω R m R m Y (x,u) (ω) - Ω Y (x,u) (ω ′ )dP (ω ′ ) Y (x ′ ,u ′ ) (ω) - Ω Y (x ′ ,u ′ ) (ω ′′ )dP (ω ′′ ) dµ(u)dµ(u ′ ) dP (ω) (19)
If the Fubini's condition

Ω R m R m Y (x,u) (ω) - Ω Y (x,u) (ω ′ )dP (ω ′ ) Y (x ′ ,u ′ ) (ω) - Ω Y (x ′ ,u ′ ) (ω ′′ )dP (ω ′′ ) dµ(u)dµ(u ′ ) dP (ω) < ∞ (20)
is satisfied (see related comment in section 3.1), it is possible to exchange the order of integration and we can rewrite eq. ( 19) as

cov Z (x, x ′ ) = R m R m Ω Y (x,u) (ω) - Ω Y (x,u) (ω ′ )dP (ω ′ ) Y (x ′ ,u ′ ) (ω) - Ω Y (x ′ ,u ′ ) (ω ′′ )dP (ω ′′ ) dP (ω)dµ(u)dµ(u ′ ) (21)
We can see that the covariance of

Z (x) (ω) is a double integral over the u × u space of the covariance of Y (x,u) (ω) cov Z (x; x ′ ) = R m R m cov Y (x, u; x ′ u ′ )dµ(u)dµ(u ′ ) (22) Therefore if X * = [x 1 * , ..., x r * ] ⊤ is a r × n matrix of r test points, then COV Z (X * ) (ω) = cov Z (x i ; x j ) 1≤i,j≤r = R m R m cov Y (x i , u; x j , u ′ )dµ(u)dµ(u ′ ) 1≤i,j≤r (23)
where cov Y (.) is given by conditional covariance of eq. ( 6). Let C denote t × t matrix cov Y (X, U; X, U) -1 , then we can write

COV Z t (X * ) (ω) = R m R m cov Y (x i * , u; x j * , u ′ )dµ(u)dµ(u ′ ) 1≤i,j≤r - R m R m t k t l c kl cov Y (x i * , u; x l , u l )cov Y (x k , u k ; x j * , u ′ )dµ(u)dµ(u ′ ) 1≤i,j≤r . (24)
By taking the integrals of each element in the sum and by the independence of the two integrals we rewrite eq. ( 24) as

COV Z t (X * ) (ω) = R m R m cov Y (x i * , u; x j * , u ′ )dµ(u)dµ(u ′ ) 1≤i,j≤r - t k t l c kl R m cov Y (x i * , u; x l , u l )dµ(u) R m cov Y (x k , u k ; x j * , u ′ )dµ(u ′ ) 1≤i,j≤r (25)

Solutions of the integrals

To calculate the covariance and mean using formulas eq. ( 25) and eq. ( 14), we need to solve integrals in the form

I(x, u, x ′ ) = R m cov Y (x, u; x ′ , u ′ )dµ(u ′ ). ( 26
)
The covariance is typically built as a product of covariances for each dimension

cov Y (x, u; x ′ , u ′ ) = n i=1 cov Y (x i ; x ′ i ) m i=1 cov Y (u i ; u ′ i ). ( 27
)
Since we consider stationary processes, we have

m i=1 cov Y (u i ; u ′ i ) = k(u ′ -u)
and for dµ(u ′ ) we have a probability measure with density function g(u ′ ), then using the definition of convolutions

I(x, u, x ′ ) = n i=1 cov Y (x i ; x ′ i ) R m k(u ′ -u)g(u ′ )d(u ′ ) = n i=1 cov Y (x i ; x ′ i ) g * k(u). (28)
We can view the new covariance as a weighted average of the old covariance at the point u, where the weighting function g(-u ′ ) is the density function of U . To approximately solve the integral for an arbitrary kernel and normally independently distributed U , it is possible to use Taylor approximation in a fashion similar as in [START_REF] Girard | Approximate methods for propagation of uncertainty with gaussian process models[END_REF].

Analytical solutions of the integrals for Gaussian kernel and independently normally distributed U

Here we show the analytical solution of the integrals for Gaussian kernels in the case of independently normally distributed U . For simplification we assume that the mean function m Y (x) is constant in u, therefore the solution of the first term in eq. ( 14) reduces to

R m m Y (X * , u)dµ(u) = m Y (X * )

Integral of the covariance

We further assume that cov Y (u; u ′ ) in eq. ( 26) is a Gaussian kernel of the form

cov Y (u i ; u ′ i ) = exp - 1 2 u i -u ′ i θ i 2 . ( 29
)
U is a random variable that comes from a multivariate normal distribution with density function f N (µ,Σ) , where µ = µ i is the mean and Σ = diag(σ 2 i ) is the covariance matrix. We look at the product of Gaussian covariances as a density function

m i=1 cov Y (u i ; u ′ i ) = (2π) m/2 |D| 1/2 f N (0,D) (u ′ -u) (30)
where D = diag(θ 2 i ), then using eq. ( 28) the integral eq. ( 26) becomes

I(x, u, x ′ ) = n i=1 cov Y (x i ; x ′ i ) (2π) m/2 |D| 1/2 R m f N (0,D) (u -u ′ )f N (µ,Σ) (u ′ )du ′ = n i=1 cov Y (x i ; x ′ i ) (2π) m/2 |D| 1/2 f N (µ,Σ) * f N (0,D) (u). ( 31
)
The convolution of two independent normal distributions is also normal,

f N (µ,Σ) * f N (0,D) (u) = f N (µ,Σ+D) (u) (32)
or

I(x, u, x ′ ) = n i=1 cov Y (x i ; x ′ i ) (2π) m/2 |D| 1/2 f N (µ,Σ+D) (u) = n i=1 cov Y (x i ; x ′ i ) |D| 1/2 |Σ + D| 1/2 exp - 1 2 (u -µ) T (Σ + D) -1 (u -µ) . (33)
A similar result is obtained under the name of Bayes-Hermite Quadrature [START_REF] O'hagan | Bayes-hermite quadrature[END_REF].

Double integral of the covariance

To calculate the covariance of Z using eq. ( 25), we also need the analytical solution of the integral

I(x, x ′ ) = R m R m cov Y (x, u; x ′ , u ′ )dµ(u ′ )dµ(u). ( 34
)
Using the same reasoning as in the previous section, we write

I(x, x ′ ) = n i=1 cov Y (x i ; x ′ i ) |D| 1/2 |Σ + D| 1/2 R m exp - 1 2 (u -µ) T (Σ + D) -1 (u -µ) dµ(u) = n i=1 cov Y (x i ; x ′ i ) (2π) m/2 |D| 1/2 R m f N (0,Σ+D) (µ -u)f N (µ,Σ) (u)du = n i=1 cov Y (x i ; x ′ i ) (2π) m/2 |D| 1/2 f N (0,Σ+D) * f N (µ,Σ) (µ) = n i=1 cov Y (x i ; x ′ i ) (2π) m/2 |D| 1/2 f N (µ,2Σ+D) (µ) = n i=1 cov Y (x i ; x ′ i ) |D| 1/2 |2Σ + D| 1/2 . ( 35
)
The results of eq. ( 33) and eq. ( 35) are identical to those presented in [START_REF] Girard | Approximate methods for propagation of uncertainty with gaussian process models[END_REF] (page 54). As mentioned in [START_REF] Girard | Approximate methods for propagation of uncertainty with gaussian process models[END_REF] the new kernel associated to Z t is a function of Σ + D = diag θ 2 1 + σ 2 1 , ..., θ 2 m + σ 2 m , so that each parameter of the covariance function is weighted by the corresponding uncertainty of U in the given dimension.

2D example

Our intention here is to illustrate the approximation of E U [f (x, U )] using Z t (x) (ω) and Y (x,u) (ω) when U ∼ N (µ, σ 2 ). In this case we assume that (x, u) ∈ R × R and that the simulator f (x, u) is the Camel back function

f (x, u) = (4 -2.1u 2 + 1 3 u 4 )u 2 + xu + (-4 + 4x 2 )x 2 (36)
in the interval [-1, -1] 2 (fig. 1). We use two different data sets: a 7 and a 20 points Latin Hypercube (LH) design optimized by MSE criterion [START_REF] Auzins | Direct optimization of experimental designs[END_REF]. The R 2 approximations measures of f (x, u) by E[Y t (x,u) (ω)] for both data sets are 0.6073 and 0.9842, respectively. We consider two uncertainty cases: in the first case, µ = 0.5 and σ = 0.1, and in the second case, µ = 0.05 and σ = 0.2. Instances of the function f (x, U ) are shown in fig. 1, by independently sampling u from each distribution.

The results for both cases are shown in fig. 2. The first and second rows corresponds to the first uncertainty case and the third and fourth -to the second. The targeted E U [f (x, U )] is obtained by numerical integration. Figure 2 confirms, that the prediction accuracy of Z t for E U [f (x, U )] depends on our knowledge about the underlying function f (x, u). Adding data points improves the approximation of f (x, u), i.e., takes

E(Z t (x) (ω)) closer to E U [f (x, U )] and decreases VAR(Z t (x) (ω)). (a) DOE, E(Y t (x,u) (ω)), 1 st case (b) E(Z t (x) (ω)) and E(f (x, U )) (c) VAR(Z t (x) (ω)) -1.0 -0.
(d) DOE, E(Y t (x,u) (ω)), 1 st case (e) E(Z t (x) (ω)) and E(f (x, U )) (f) VAR(Z t (x) (ω)) (g) DOE, E(Y t (x,u) (ω)), 2 nd case (h) E(Z t (x) (ω)) and E(f (x, U )) (i) VAR(Z t (x) (ω)) (j) DOE, E(Y t (x,u) (ω)), 2 nd case (k) E(Z t (x) (ω)) and E(f (x, U )) (l) VAR(Z t (x) (ω))
Figure 2: Kriging approximations of f (x, u) and E U (f (x, U )) for the Camel back test function. Left column, DOE (squares) and contour lines of predicted mean E(Y t (x,u) (ω)). Middle column, predicted E(Z t (x) (ω)) (dashed) and E U (f (x, U )) (solid). Right column, predicted variance VAR(Z t (x) (ω)). First and second rows, U ∼ N (0.5, 0.01), third and fourth rows, U ∼ N (0.05, 0.04).

Infill strategy for robust global optimization

In order to solve the optimization problem of eq. ( 1), we utilize the statistical models of f (x, u) and E U [f (x, U )] (eq. ( 4) and eq. ( 7)) to define an infill sampling criterion [START_REF] James | Flexibility and Efficiency Enhancements for Constrained Global Design Optimization with Kriging Approximations[END_REF]. We follow the EGO philosophy [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] that employs kriging models as approximations of a deterministic function and provides a measure of local uncertainty. We propose an infill strategy that consists of two steps: at first we choose x next by maximizing an Expected Improvement of the projected process Z t (x) (ω), eq. ( 7). x next is an interesting point in the x (deterministic) space. Then we choose a new point for the actual simulation (x t+1 , u t+1 ), in (x, u) (joint) space by minimizing the variance of the projected process Z t+1 (x) (ω) at the interesting point x next . Two goals are thus pursued simultaneously: a rapid minimization of the objective (mean response), and a spare use of the real (computationally costly) simulator f by a controlled sampling of the uncertain parameters.

Maximization of Expected Improvement

We define the improvement as

I t x (ω) = f t min -Z t (x) (ω) + , (37)
where [.] + denotes max(0, .). As we do not have observations of E U [f (x, U )], we say that f t min is a minimum of the predicted mean of the projected process Z t (x) (ω)

f t min = min x E Z t (x) (ω) . ( 38
)
The Expected Improvement is simply

EI Z (x) = E I t x (ω) = E f t min -Z t (x) (ω) + ( 39
)
for which we have an analytical expression in [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] EI

Z (x) = VAR(Z t (x) (ω)) [uΦ(u) + φ(u)] , (40)
where

u = f t min -E(Z t (x) (ω)) VAR(Z t (x) (ω)
) and where Φ and φ are the normal cumulative distribution function and density function, respectively.

By maximizing EI eq. ( 40) we obtain an interesting point in the x space,

x next = arg max x EI Z (x). ( 41
)

Minimization of Variance

Maximizing the EI gives some insight (x next ) on where to sample in the x space in order to optimize (min x E U [f (x, U )]). However, before the simulator f can be called, u must be chosen. In order to sample points in the (x, u) space, we propose to use the point that reduces the most the uncertainty in the estimated E U [f (x, U )] at the point of interest

x next . The process Z t approximates E U [f (x, U )].
It is a projection of Y t from the (x, u) to the x space and the variance of Z t at a given x depends on our knowledge about the simulator around our point of interest in the (x, u) space. We choose the next point for simulation (x t+1 , u t+1 ) such that

(x t+1 , u t+1 ) = arg min x,u VAR Z t+1 (x next ) (ω). ( 42
)
Note that, such sampling criterion follows the principle of deterministic EGO, where sampling at a new point reduces the uncertainty at this point to zero, thus eventually allowing sample points outside of current minimum in unexplored areas. The sampling criterion of eq. ( 42) involves a non linear optimization. With the estimation of kriging hyper parameters, the determination of f t min eq. ( 38), the maximization of EI eq. ( 41), this is the fourth non linear optimization of the overall procedure, which therefore ends up being somewhat calculation intensive. Note however, that none of these optimization tasks involve calls to the simulator f . They only work with the underlying regression models Y and Z. This is the reason why the method described in this paper, like all kriging based methods, is suitable for computationally expesive simulators. From eq. ( 25) we obtain

VAR Z t+1 (x next ) (ω) = R m R m cov Y (x next , u; x next , u ′ )dµ(u)dµ(u ′ ) - R m cov Y (x next , u; X t+1 , U t+1 )dµ(u) cov Y (X t+1 , U t+1 ; X t+1 , U t+1 ) -1 R m cov Y (X t+1 , U t+1 ; x next , u ′ )dµ(u ′ ), (43)
where X t+1 = (X t , x t+1 ) and U t+1 = (U t , u t+1 ) and X t , U t are observed data points. In VAR Z t+1 (eq. ( 42), eq. ( 43)) the hyper parameters of the underlying Gaussian process Y are estimated using only the t previously observed data points (X t , U t , Y t ).

In order to solve the minimization task of eq. ( 42) we need to calculate eq. ( 43) for different (x t+1 , u t+1 ). The first integral in eq. ( 43) is constant with respect to (x t+1 , u t+1 ) and we denote it as

σ 2 Z σ 2 Z = R m R m cov Y (x next , u; x next , u ′ )dµ(u)dµ(u ′ ) .
We also define the vector

k(x t+1 , u t+1 ) = R m cov Y (X t+1 , U t+1 ; x next , u)dµ(u), the matrix C t+1 = cov Y (X t+1 , U t+1 ; X t+1 , U t+1 )
and the upper triangular matrix T t+1 obtained by Cholesky decomposition such that T t+1′ T t+1 = C t+1 (where ′ denotes transposition). Now we can express the variance in eq. ( 43) as

VARZ t+1 (x next ) (ω) =σ 2 Z -k(x t+1 , u t+1 ) ′ (C t+1 ) -1 k(x t+1 , u t+1 ) = σ 2 Z -k(x t+1 , u t+1 ) ′ (T t+1 ) -1 (T t+1′ ) -1 k(x t+1 , u t+1 ) = σ 2 Z -((T t+1′ ) -1 k(x t+1 , u t+1 )) ′ (T t+1′ ) -1 k(x t+1 , u t+1 ). (44)
To calculate the above variance we only need to solve the linear system

(T t+1′ )b = k ′ .
Cholesky decomposition has complexity O((t+1) 3 ), and it becomes time consuming as t increases. To further decrease the computation time, we use Cholesky decomposition of block matrices of C t+1 ,

C t+1 = C t c t+1 c t+1′ s
where c t+1 = cov Y (X t , U t ; x t+1 , u t+1 ) and s = cov Y (x t+1 , u t+1 ; x t+1 , u t+1 ). Let T t be the upper triangular matrix obtained by Cholesky decomposition T t′ T t = C t . T t is independent of (x t+1 , u t+1 ). T t+1 can be expressed in terms of T t for any (x t+1 , u t+1 ) by using T t , (T t ) -1 and matrix multiplication as

T t+1 = T t (T t′ ) -1 c t+1 O s -c t+1′ (T t ) -1 (T t′ ) -1 c t+1 1/2 .
With the block matrix decomposition, the costly O(t 3 ) Cholesky decomposition needs to be performed only once per optimization (eq. ( 42)). The criterion eq. ( 42) will typically choose a point (x t+1 , u t+1 ) that is close to (x next , µ). However, once sampling near (x next , µ) has been done, it may also choose points that are away from (x next , µ). We illustrate such situation in fig. 3 and fig. 4, where we show two cases with different initial DOEs and x next . We show contour lines and surfaces plots of VAR (Z t+1 (x next ) (ω)) and the underlying VAR (Y t (x,u) (ω)). In the first case the x next = -0.7. The corresponding surface of the variance to minimize is shown in the panels (c) and (d) of fig. 3 . In this case the minimum (x t+1 , u t+1 ) (star) coincide with the (x next , µ) (dimond). In the second case the new point is x next = -0.5. The corresponding surface of the variance to minimize is shown in the panels (c) and (d) of fig. 4. As we already have some points in the (x, u) space around (x next , µ) (diamond), the minimum point (star) of VAR (Z t+1 (x next ) (ω)) is away from the (x next , µ) moving towards the direction of higher uncertainty in Y t (x,u) (ω). (x next ) (ω)) as function of (x t+1 , u t+1 ), the optimum (x t+1 , u t+1 ) (star) and point (x next , µ) (diamond). (d): 3D surface of predicted variance VAR(Z t+1 (x next ) (ω)). Notice how (x next , µ) and (x t+1 , u t+1 ) coincide.

(a) DOE, (x next , µ), (u t+1 , x t+1 ), VAR(Y t (x,u) (ω)) (b) VAR(Y t (x,u) (ω)) (x next , µ), (u t+1 , x t+1 ), VAR(Z t+1 (x next ) (x, u)) (d) VAR(Z t+1 (x next ) (x, u))

Implementation and Results

The proposed approach to solve eq. ( 1) may be summarized in an algorithm 1.

Algorithm 1 Compleate method (EI-VAR)

Create an initial DOE in the (x, u) space and obtain simulator responses. while stopping criterion not met do Create regression model Y t (x,u) (ω) of the deterministic simulator in the joint (x, u) space; Obtain regression model Z t (x) (ω) in the deterministic space (x), using covariance information of Y t (x,u) (ω); Maximize EI Z (x) to choose x next ; Minimize VAR Z t+1 (x next ) (ω) to obtain the next point (x t+1 , u t+1 ); Calculate simulator response at the next point f (x t+1 , u t+1 ) and update DOE.

end while

The proposed algorithm involves 4 optimization subtasks at every iteration:

1. Maximizing the Likelihood to obtain covariance parameters and create Y (x,u) (ω). (x next ) (ω)) as function of (x t+1 , u t+1 ), the optimum (x t+1 , u t+1 ) (star) and point (x next , µ) (diamond). (d): 3D surface of predicted variance VAR(Z t+1 (x next ) (ω)). Notice how (x next , µ) and (x t+1 , u t+1 ) differ.

(a) DOE, (x next , µ), (u t+1 , x t+1 ), VAR(Y t (x,u) (ω)) (b) VAR(Y t (x,u) (ω)) (c) (x next , µ), (u t+1 , x t+1 ), VAR(Z t+1 (x next ) (x, u)) (d) VAR(Z t+1 (x next ) (x, u))

Finding the minimum

f t min = min x E Y Z t (x) (ω) .

Maximizing the expected improvement

x next = arg max x EI Z (x).
4. Obtaining the next point in the joined space

(x t+1 , u t+1 ) = arg min x,u VAR Z t+1 (x next ) (ω).
To solve these optimization problems, which are typically multi-modal we use the CMA-ES algorithm [START_REF] Hansen | The CMA evolution strategy: a comparing review[END_REF], a state-ofthe-art global optimizer. Numerical problems may be encountered during the minimization of the variance in the step 4: when the new point (x, u) is added close to existing points, the covariance matrix becomes badly conditioned. To overcome this problem, we increase the diagonal of the covariance matrix by adding a "noise" term for the calculation of the

VAR [Z t+1 (x next ) (ω)].
In order to reduce the computational burden and eliminate numerical problems encountered solving the 4th suboptimization, we investigate also simplified approach (algorithm 2) where x t+1 = x next and u t+1 is sampled from the corresponding distribution (thus eliminating the 4th sub-optimization problem).

In the following, test cases will be shown on analytical "simulators" f . It should be kept kept in mind that these test Algorithm 2 Simplified approach, with sampling u (EI-Sample) Create an initial DOE in the (x, u) space and obtain simulator responses. while stopping criterion not met do Create regression model Y t (x,u) (ω) of the deterministic simulator in the joint (x, u) space; Obtain regression model Z t (x) (ω) in the deterministic space (x), using covariance information of Y t (x,u) (ω); Maximize EI Z (x)(ω)) to choose x next ; Sample u t+1 from the corresponding distribution; Calculate simulator response at the next point f (x next , u t+1 ) and update DOE. end while cases are meant to illustrate the method but that the real scope of application is expensive simulators.

Graphical illustration on the 2D Camel back test function

To illustrate behaviour of the proposed algorithm we use the same test case as in section 4.2, where x, u ∈ R × R, the simulator f (x, u) is the Camel back function eq. ( 36) in the interval [-1, -1] 2 . The initial DOE consists of 5 points (one point in each corner and one point in the center). We consider two uncertainty cases: µ = 0.5 and σ = 0.1 in the first case and µ = 0.05 and σ = 0.2, in the second case. Figure 5 and fig. 6 show some steps of the optimization. It can be seen that already after the 5th iteration the predicted mean E(Z (x) (ω)) captures quite well the shape of the objective E(f (x, U )) and the predicted variance VAR(Z (x) (ω)) is quite small around the optimum x * = arg min x E U (f (x, U )). During the next iterations points are added in the optimal region and the approximating surface becomes there even more precise. In order to reduce the model uncertainty, at the iterations 7, 13 and 16 the EI chooses points that are away from the current optimum. After theses steps the uncertainty about the model is very small and for the subsequent iterations EI chooses points that are in the vicinity of local minimum. It is interesting to notice that during the initial iterations the points are spread in the x space but for the u's are close to the mean µ of U . After some iterations when the variance at the point of maximum EI has been reduced, the new points may be away from the (x next , µ) and tend to spread more in the (x, u) space. Figure 7 provides a comparison between the complete method EI-VAR (algorithm 1) and simplified approach EI-Sample (algorithm 2). The first row of fig. 7 corresponds to the DOE of the EI-Sample in the first case (U ∼ N (0.5, 0.01)). The second and the third rows of fig. 7 show the DOE in the second case (U ∼ N (0.05, 0.04)) of the complete (EI-VAR) and simplified (EI-Sample) approach, respectively. As can be seein in fig. 6 and fig. 7, even though the proposed method is global, the DOE produced are focused on the region of the optimum (x * ). Although, in the current case, such early focus is a sign of efficiency, it may become an issue of premature convergence to a non global optimum on more misleading problems. The early convergence is a typical feature of the one step lookahead EI based optimizers [START_REF] Donald | A taxonomy of global optimization methods based on response surfaces[END_REF], [START_REF] Forrester | Global optimization of deceptive functions with sparse sampling[END_REF], [START_REF] Ginsbourger | Kriging is well-suited to parallelize optimization[END_REF]: if the variance far from the current optimum is small, the EI tends to zero in these regions. A good practice when using EI based method is to use an initial DOE which is large enough to capture the overall trend (complexity) of the underlying simulator. Note that the complete method (EI-VAR) spreads points more away from the (x next , µ). Such a better space filling property explains why, in fig. 8,E(Y t (x,u) (ω)) converges faster, in terms of square error, to f (x, u) with the complete approach (EI-VAR) then with the simplified one (EI-Sample).

Comparisons

In the rest of this section we illustrate performance of the proposed method on analytical test cases and compare it to the MC based optimization.

Comparison test cases

The test cases are based on the Michalewitz's function

f (x) = - d i=1 sin(x i ) sin( ix 2 i π ) 2 ,
where x is in the interval [0, π] d . The simulator is composed as an additive function Each row represents an iteration (between the 1st and the 5th). Left column: contour lines of predicted mean E(Y t (x,u) (ω)), initial DOE (squares), point (x next , µ) where x next = arg max x EI Z (x) (diamond), new point (x t+1 , u t+1 ) where (x t+1 , u t+1 ) = arg min x,u VAR Z t+1 (x next ) (ω) (star), points from previous iterations (circles). Middle column: objective E(f (x, U )) (solid) and predicted mean E(Z t (x) (ω)) (dashed). Right column: predicted variance VAR(Z t (x) (ω)) (solid) and rescaled EI Z (x) (dashed). (x next ) (ω) (star), points from previous iterations (circles). Middle column: objective E(f (x, U )) (solid) and predicted mean E(Z t (x) (ω)) (dashed). Right column: predicted variance VAR(Z t (x) (ω)) (solid) and rescaled EI Z (x) (dashed).

f (x, u) = f (x) + f (u). (a) Initial DOE, E(Y t (x,u) (ω)) (b) Initial E(Z t (x) (ω)) (c) Initial VAR(Z t (x) (ω)), EI Z (x) (d) 1st iter. E(Y t (x,u) (ω)) 1st E(Z t (x) (ω)) (f) 1st iter. VAR(Z t (x) (ω)), EI Z (x) (g) 3rd iter. E(Y t (x,u) (ω)) (h) 3rd iter. E(Z t (x) (ω)) (i) 3rd iter. VAR(Z t (x) (ω)), EI Z (x) (j) 5th iter. E(Y t (x,u) (ω)) (k) 5th iter. E(Z t (x) (ω)) (l) 5th iter. VAR(Z t (x) (ω)), EI Z (x)
(a) 7th iter. E(Y t (x,u) (ω)) (b) 7th iter. E(Z t (x) (ω)) (c) 7th iter. VAR(Z t (x) (ω)), EI Z (x) (d) 13th iter. E(Y t (x,u) (ω)) (e) 13th iter. E(Z t (x) (ω)) (f) 13th iter. VAR(Z t (x) (ω)), EI Z (x) (g) 16th iter. E(Y t (x,u) (ω)) (h) 16th iter. E(Z t (x) (ω)) (i) 16th iter. VAR(Z t (x) (ω)), EI Z (x) (j) 50th iter. E(Y t (x,u) (ω)) (k) 50th iter. E(Z t (x) (ω)) (l) iter. VAR(Z t (x) (ω)), EI Z (x)
(a) 5Oth iter. E(Y t (x,u) (ω)) (b) 5Oth iter. E(Z t (x) (ω)) (c) 5Oth iter. VAR(Z t (x) (ω)), EI Z (x) 5Oth iter. E(Y t (x,u) (ω)) (e) 5Oth iter. E(Z t (x) (ω)) (f) 5Oth iter. VAR(Z t (x) (ω)), EI Z (x) (g) 50th iter. E(Y t (x,u) (ω)) (h) 50th iter. E(Z t (x) (ω)) (i) 50th iter. VAR(Z t (x) (ω)), EI Z (x)
Figure 7: 50th iterations on the Camel back test case. 1st row: U ∼ N (0.5, 0.01), simplified approach (EI-Sample). 2nd and 3rd rows: U ∼ N (0.05, 0.04), complete (EI-VAR) and simplified (EI-Sample) approach, respectively. Left column: contour lines of predicted mean E(Y t (x,u) (ω)), initial DOE (squares), point (x next , µ) where x next = arg max x EI Z (x) (diamond), new point (x t+1 , u t+1 ) where (x t+1 , u t+1 ) = arg min x,u VAR Z t+1 (x next ) (ω) (star), points from previous iterations (circles). Middle column: objective E(f (x, U )) (solid) and predicted mean E(Z t (x) (ω)) (dashed). Right column: predicted variance VAR(Z t (x) (ω)) (solid) and rescaled EI Z (x) (dashed). The complexity (multi-modality) of the test case is dependent on the number of dimensions n and m of x and u, respectively. The components for 3 dimensions are shown in fig. 9.

For comparison we use three test cases with equal number of x and u dimensions:

• 2 dimensional test case (n = 1 and u = 1), where µ = 1.5 and σ = 0.2;

• 4 dimensional test case (n = 2 and u = 2), where µ = (1.5, 2.1) ′ ; and Σ = diag(0.2 2 , 0.2 2 );

• 6 dimensional test case (n = 3 and u = 3), where µ = (1.5, 2.1, 2) ′ ; and Σ = diag(0.2 2 , 0.2 2 , 0.3 2 ).

Settings of internal optimization

The proposed method is a deterministic algorithm which involves several optimization subtasks. To solve these tasks we use a stochastic optimization method CMA-ES [START_REF] Hansen | The CMA evolution strategy: a comparing review[END_REF], which introduces some randomness in the proposed method.

The CMA-ES optimization algorithm has three parameters: a population size λ, an initial standard deviation σ 0 and a maximum number of calls to the objective t max . For each of the four optimization subtasks, the setting are the following:

• for finding kriging hyper parameters (maximizing likelihood) t max = 200 + 100(n + m), λ = 4 + ⌊3 log(n + m)⌋ and σ 0 is one third of the domain.

• for finding the minimum of E Y Z t (x) (ω) and maximizing the expected improvement EI Z (x), t max = 300nλ, λ = 4 + ⌊3 log(n)⌋n and initial standard deviation σ 0 is one half of the domain. • for minimizing the variance VARZ t+1 (x next ) (ω),t max = 300(n + m)λ where population size λ = 4 + ⌊3 log(n + m)⌋ and σ 0 is one half of the domain.

Using such optimization strategy for 6 dimensional test case on 3.6 GHz CPU required about 30 to 200 seconds for hyperparameter estimation, 8 to 180 seconds for minimization of the mean of projected process, 8 to 30 seconds for the maximization of the EI and 17 to 50 seconds for the minimization of the variance at the next point of interest for every iteration of the proposed method. This strategy is not optimal and alternative study to explore the best optimization procedure for every sub-optimization would be needed. In fig. 10 median, 80th and 20th percentiles are shown for 6D test case using 10 independent runs of the proposed method (EI-VAR) with fixed initial design (MSE optimized LH [START_REF] Auzins | Direct optimization of experimental designs[END_REF]) for all runs. It reflects the variation in the performance of the proposed method due to random sub-optimization procedures.

Comparison methods

To illustrate performance of the proposed method we compare it to the MC based optimization. We use EI and kriging model with noisy MC estimations of the objective values E U (f (x, U )). The precision (fidelity level) of the MC estimation of the objective is limited by the inverse of the number of calls k to the underlying simulator f (x, u). For illustration and in order to keep the procedure as simple as possible, the k values are kept constant. In order to approximate the E U (f (x, U )) and filter the error of MC estimate, we use kriging model with Gaussian covariance kernel and homogeneous noise which is added on the diagonal of the covariance matrix. We assume that this noise is not known and together with the kernel hyper parameters is estimated by maximizing the Likelihood. To choose the next point for MC estimation we use maximum EI, where an improvement is defined as

I t x (ω) = f t min -Y t (x) (ω) + ,
and Y t (x) (ω) is a kriging model of the noisy MC estimate conditioned on the nosy observations Y. The current minimum is the best noisy observation f t min = min(Y). For the maximization of the Likelihood and EI the CMA-ES is used. The method is summarized in algorithm 3. Note that the kriging model is built in n dimensions as the objective E U (f (x, U )) is function only of x. The point which is chosen by EI leads to k MC calls to the simulator f (x, u) at every iteration.

Experimental procedure

To compare the performance of the various methods we use the absolute distance between the optimum given by the method and the true minimum |xx * | where x * = arg min x E U [f (x, U )] versus number of calls to the simulator f (x, u). As we optimize stochastic functions we perform 10 runs of each method and use the average. The Random Latin Hypercube design (RLHD) was used as an initial DOE for all methods creating new initial design for every run. The A 50 iterations budget is used as a stopping criterion for EI-VAR and EI-Sample and 30 iterations budget for EI-MC.

Results

The results are summarized in fig. 11, fig. 12 and fig. 13. The mean of the convergence for every dimension and total for all dimensions of x are given. The confidence intervals are illustrated using 20th and 80th percentile. The main conclusion from the results of the test cases is that the convergence rates for the proposed approach (EI-VAR) are much better in comparison to the MC based optimization (EI-MC). The precision of the solutions of the MC based approaches is limited by the number of MC runs and in some cases (low number of MC runs k) did not converge to the true optimum (see for example fig. 13). For MC based approaches (EI-MC), where higher number of MC samples k where used [START_REF] Beyer | Robust optimization -a comprehensive survey[END_REF][START_REF] Eiben | Introduction to Evolutionary Computing (Natural Computing Series)[END_REF][START_REF] Huang | Global optimization of stochastic black-box systems via sequential kriging meta-models[END_REF] (thus providing more precise estimation of the objective), the proposed method (EI-VAR) has converged even before the initial DOE for EI-MC has been obtained. By definition EI-VAR is O(k) times less expensive for iteration and it indicates the potential benefit of the proposed approach as usually MC methods require at least hundreds of runs per estimation. Notice also that the confidence intervals for the (EI-VAR and EI-Sample) are smaller compared to the MC based approach and tend to decrease with the number of sampled points. The convergence of the simplified approach (EI-Sample) seems to be similar to the complete method (EI-VAR) for 2D and 4D test cases. However with more complex model as in 6D case, the results of simplified approach are much worse than the complete approach and may even be out performed by MC based optimization (this can be observed for the 3rd component of x fig. 13), therefore the complete approach (EI-VAR) is safer and preferable to the simplified one (EI-Sample).

Finally in fig. 14 we show 20th and 80th percentile of convergence rates for 10 runs of the complete approach (EI-VAR) on the 6D test function for two cases: initial DOE is fixed MSE optimizd LH and initial DOE is RLH independently generated for every run. The fixed design leads to slightly better convergence and slightly less variance in the results. The results indicate that for the given test function the method is somehow robust to initial DOE and most of the variance in the results are due to the stochastic nature of the sub-optimization procedures.

Conclusion and Remarks

In this paper, a new robust optimization method was proposed for problems where some of the input parameters are random. The method minimizes the expected output of the simulator by, first, creating a kriging model in the joined space of deterministic and probabilistic input parameters. Then another "projected" Gaussian process is built to represent the mean output of the simulator. The maximization of the EI together with the minimization of the variance are used to define a sampling criterion for the optimization and simultaneous sampling in the uncertain parameter space. The method is compared with simplified MC based optimization on several test cases and the results indicate significant enhancement in terms of precision and number of calls to the simulator. Another advantage of the proposed method is that one does not have to choose the number of the MC simulations. Furthermore the kriging model of the simulator and expectation measure are obtained, which may be used for other analyses. For example one can use these models to numerically obtain other statistical measures such as variance. The potential limitations are associated with the fixed covariance kernel and noise models. The performance of the proposed method probably is dependant on the ability of the regression model Y to model the true simulator, therefore it may encounter problems when the Gaussian kernel is not good choice for the covariance function. The Gaussian covariance kernel yields very smooth (infinitely differentiable) trajectories. Such strong smoothness assumptions are Figure 14: The convergence of the complete approach (EI-VAR). 20th and 80th percentile of 6D test case using 10 runs with fixed initial design (dimond) and generating new random LH for every run (star).

unrealistic for modelling many physical processes [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning[END_REF]. The application to other analytical and practical test problems as well as possibility of using other kernel functions and other noise models is a subject for further research.

From practical point of view the method may be applicable only to costly simulators, as it involves four global optimization subtasks. For effective implementation the study of efficient strategies to solve these subtasks would be needed.

Figure 1 :

 1 Figure 1: Noisy Camel back test case. Plots (b) and (c) show f (x, u), where the u's are i.i.d. for each new x.

Figure 3 :

 3 Figure 3: Illustration of the sampling criterion on the Camel back function. Panel (a): contour lines of predicted variance VAR(Y t (x,u) (ω)), initial DOE (squares), point (x next , µ) (diamond) and point (x t+1 , u t+1 ) = arg min VAR (Z (x next ) (ω)) (star). (b) : 3D surface of the predicted variance VAR(Y t (x,u) (ω)). (c): contour lines of the predicted variance VAR(Z t+1(x next ) (ω)) as function of (x t+1 , u t+1 ), the optimum (x t+1 , u t+1 ) (star) and point (x next , µ) (diamond). (d): 3D surface of predicted variance VAR(Z t+1 (x next ) (ω)). Notice how (x next , µ) and (x t+1 , u t+1 ) coincide.

Figure 4 :

 4 Figure 4: Illustration of the sampling criterion on the Camel back function. Panel (a): contour lines of predicted variance VAR(Y t (x,u) (ω)), initial DOE (squares), point (x next , µ) (diamond) and point (x t+1 , u t+1 ) = arg min VAR (Z (x next ) (ω)) (star). (b): 3D surface of the predicted variance VAR(Y t (x,u) (ω)). (c): contour lines of the predicted variance VAR(Z t+1(x next ) (ω)) as function of (x t+1 , u t+1 ), the optimum (x t+1 , u t+1 ) (star) and point (x next , µ) (diamond). (d): 3D surface of predicted variance VAR(Z t+1 (x next ) (ω)). Notice how (x next , µ) and (x t+1 , u t+1 ) differ.

Figure 5 :

 5 Figure 5: Illustration of the simultaneous kriging based algorithm on the 2D Camel back function (U ∼ N (0.5, 0.01)).Each row represents an iteration (between the 1st and the 5th). Left column: contour lines of predicted mean E(Y t (x,u) (ω)), initial DOE (squares), point (x next , µ) where x next = arg max x EI Z (x) (diamond), new point (x t+1 , u t+1 ) where (x t+1 , u t+1 ) = arg min x,u VAR Z t+1 (x next ) (ω) (star), points from previous iterations (circles). Middle column: objective E(f (x, U )) (solid) and predicted mean E(Z t (x) (ω)) (dashed). Right column: predicted variance VAR(Z t (x) (ω)) (solid) and rescaled EI Z (x) (dashed).

Figure 6 :

 6 Figure 6: Illustration of the simultaneous kriging based algorithm on the 2D Camel back function (U ∼ N (0.5, 0.01)). Continuation of fig. 5. Left column: contour lines of predicted mean E(Y t(x,u) (ω)), initial DOE (squares), point (x next , µ) where x next = arg max x EI Z (x) (diamond), new point (x t+1 , u t+1 ) where (x t+1 , u t+1 ) = arg min x,u VAR Z t+1 (x next ) (ω) (star), points from previous iterations (circles). Middle column: objective E(f (x, U )) (solid) and predicted mean E(Z t (x) (ω)) (dashed). Right column: predicted variance VAR(Z t (x) (ω)) (solid) and rescaled EI Z (x) (dashed).

( a )Figure 8 :

 a8 Figure 8: MSE error of E(Y t (x,u) (ω)), compared to f (x, u), where f is the Camel back function. Solid line: complete method (EI-VAR), dashed line : simplified method (EI-Sample).

Figure 9 :

 9 Figure 9: Components of Michalewitz test function for the first, second and third dimension

Figure 10 :

 10 Figure 10: Median (solid), 80th (diamond) and 20th (star) percentiles for 6D test case using 10 independent runs with fixed initial DOE for all runs.

Algorithm 3

 3 EI based on MC estimations of E U (f (x, U )) (EI-MC) Create an initial DOE in the x space. Estimate E U (f (x, U )) by k MC runs to the f (x, u). while stopping criterion not met do Create kriging model Y t (x) (ω) of the noisy simulator in the x space, by maximizing likelihood; Maximize EI Y (x)) to choose x next ; Estimate objective E U (f (x next, U )) by k MC runs to the simulator f (x, u); Update DOE. end while number of points for initial DOE is t = 5d, where d is a number of dimensions (d = n + m for the EI-VAR and EI-Sample and d = n for MC based approaches EI-MC). When d = 1 the number of points for initial DOE t = 3 (2D test case for EI-MC).

  (a) Mean convergence results (b) convergence (zoomded in) (c) 20th and 80th percentile

Figure 11 :

 11 Figure 11: Convergence rates of different methods for 2D Michalewitz test function (n = 1, m = 1). The EI-MC with fixed nummber of runs k = 3 (diamond), k = 5 (square), k = 10 (star), k = 20 (dotted), complete method by minimizing variance (EI-VAR) (solid) and simplified approach by sampling u (EI-Sample) (cross).

  (a) Mean convergence for x 1 (b) Mean convergence for x 2 (c) Mean convergence for all dimensions (d) Mean convergence for all dimensions (zoomed in) (e) 20th and 80th percentile of convergence (f) 20th and 80th percentile of convergence (zoomed in)

Figure 12 :

 12 Figure 12: Convergence rates of different methods for 4D Michalewitz test function (n = 2, m = 2). The EI-MC with fixed nummber of runs k = 3 (diamond), k = 5 (square), k = 10 (star), k = 20 (dotted), complete method by minimizing variance (EI-VAR) (solid) and simplified approach by sampling u (EI-Sample) (cross).

(a) Mean convergence for x 1 (b) Mean convergence for x 2 (c) Mean convergence for x 3 (Figure 13 :

 12313 Figure 13: Convergence rates of different methods for 6D Michalewitz test function (n = 3, m = 3). The EI-MC with fixed nummber of runs k = 3 (diamond), k = 5 (square), k = 10 (star), k = 20 (dotted), complete method by minimizing variance (EI-VAR) (solid) and simplified approach by sampling u (EI-Sample) (cross).