F Zhang

B Cosson

S Comas-Cardona

C Binetruy

Efficient stochastic simulation approach for RTM process with random fibrous permeability

Keywords: B. Permeability, C. Stochastic simulation, C. Random field, C. Level set method (LSM), E. Mold filling process

An efficient numerical approach is presented for the simulation of mold filling process with random fibrous permeability as input. Suitable dimension-reduction techniques are employed to represent all the random fields in the physical system. The required accuracy in tracking the moving resin flow front is guaranteed by the Level Set Method (LSM). The Probabilistic Collocation Method (PCM) is applied in the efficient solution of the non-linear stochastic system. Numerical examples for different injection schemes are presented to demonstrate the cost-effectiveness of the current approach in predicting the variability in mold filling results in comparison to the traditional Monte Carlo Simulation (MCS).

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Introduction

The resin transfer molding (RTM) technique, consisting of the impregnation of a fibrous preform by a reactive resin and the curing of the part, is widely applied in advanced composites manufacturing. Irregular flow patterns or defects in the final products often occur due to the high variability in the transport properties of the fibrous media.

Therefore, RTM robustness has to be studied to avoid trial and error procedures.

Experimental studies [START_REF] Lundström | In-plane permeability measurements: a nordic round-robin study[END_REF][START_REF] Hoes | New set-up for measurement of permeability properties of fibrous reinforcements for RTM[END_REF][START_REF] Desplentere | Micro-CT characterization of variability in 3D textile architecture[END_REF][START_REF] Endruweit | Influence of the fabric architecture on the variations in experimentally determined in-plane permeability values[END_REF] show the scatter in measured global permeability data can be up to 0.5 in terms of coefficient of variation from batch to batch of the same material.

Besides, local permeability exhibits spatially-dependent fluctuation due to its inherently heterogeneous architecture. Although studies have been performed on the probabilistic characterization of the random transport properties of fibrous media [START_REF] Hoes | New set-up for measurement of permeability properties of fibrous reinforcements for RTM[END_REF][START_REF] Desplentere | Micro-CT characterization of variability in 3D textile architecture[END_REF][START_REF] Endruweit | Influence of the fabric architecture on the variations in experimentally determined in-plane permeability values[END_REF], the random effects are usually neglected in the available optimization or on-line control schemes [START_REF] Hsiao | Flow sensing and control strategies to address race-tracking disturbances in RTM. Part I: design and algorithm development[END_REF][START_REF] Nielsen | Closed-loop flow control in RTM using real-time numerical process simulations[END_REF][START_REF] Lee | Analysis and minimization of void formation during resin transfer molding process[END_REF] which employ common simulation softwares for mold filling. Available stochastic studies in the RTM context [START_REF] Markicevic | Stochastic modeling of preform heterogeneity to address dry spots formation in the VARTM Process[END_REF][START_REF] Padmanabhan | Stochastic modeling of nonisothermal flow during resin transfer molding[END_REF][START_REF] Verleye | Permeability prediction for the meso-macro coupling in the simulation of the impregnation stage of RTM[END_REF][START_REF] Li | Statistical characterization and robust design of RTM processes[END_REF][START_REF] Skordos | Stochastic simulation of woven composites forming[END_REF] mostly depend on the Monte Carlo simulation (MCS), which is computationally expensive for a fast uncertainty analysis or to be combined with other algorithms for optimization or sensitivity analysis.

In this paper, an efficient stochastic simulation tool for mold filling process is developed, which is able to accurately predict the statistics of responses at a fraction of the cost of the MCS. The approach is implemented based on available numerical techniques in the stochastic and flow simulation domain, which are explained in details in Section 2.

Numerical examples are given in Section 3 to illustrate the cost-effectiveness of the current method in uncertainty quantification, using the MCS as reference. Main conclusions and future outlook are discussed in Section 4.

Stochastic simulation approach for mold filling process

Dimension reduction for the stochastic system

As fiber reinforcement permeability is function of local architectural parameters such as fiber volume fraction and fiber tow orientation [START_REF] Desplentere | Micro-CT characterization of variability in 3D textile architecture[END_REF][START_REF] Endruweit | Influence of the fabric architecture on the variations in experimentally determined in-plane permeability values[END_REF][START_REF] Verleye | Permeability prediction for the meso-macro coupling in the simulation of the impregnation stage of RTM[END_REF][START_REF] Skordos | Stochastic simulation of woven composites forming[END_REF], it can be described by a random field ( ) , K x ω with as a variable in the random space. For this application, ( ) , K x ω is simplified as a weak stationary random field with mean ( K µ ) and autocovariance independent of the absolute location, i.e.

( ) [

]

K x K E µ ω = , (1.a) [ ] ( ) ( ) [ ] ( ) ( ) ( ) ( ) [ ] K K KK x x K x K E x x K x K Cov x C µ ω µ ω ω ω - ∆ + - = ∆ + = ∆ , , , , (1.b)
The reliable data on the marginal Probability Density Function (PDF) of fiber permeability is very limited. Some studies [START_REF] Hoes | New set-up for measurement of permeability properties of fibrous reinforcements for RTM[END_REF] presumed a Gaussian distribution for the measured global permeability data, while the log-normal distribution is adopted for ( )

, K x ω in this study, in order to ensure the nonnegativeness of permeability data even in case of a high variance.

In analogy with the discretization of a deterministic property field for numerical simulation, the random field ( )

, K x ω has to be discretized in both spatial and probabilistic spaces to reduce the dimensions of the corresponding stochastic system and to characterize it by the Stochastic Partial Differential Equations (SPDEs). Focusing on the underlying Gaussian field ( )

( ) ( ) ω ω , log , x K x Y =
, its discretization can be performed using the Karhunen-Loève expansion (KLE) [START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF], by projecting ( ) , Y x ω onto a set of uncorrelated random basis functions with corresponding deterministic coefficients, i.e.

( ) ( ) ( )

1 , M Y i i i i Y x x ω µ λ φ ξ ω = = + (2)
In Eq. ( 2), M is the order of the KLE which controls the accuracy of the expansion. It has been proved that the random basis

( ) { } i ξ ω (i = 1, 2, …, M) is a collection of 4
independent standard Gaussian random variables for a Gaussian random field [START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF]. The corresponding deterministic coefficients ( ) ( )

x x Y i i i φ λ =
are determined by solving Eq.

(3) for the eigenvalues { } i λ and eigenfunctions

( ) { } x i φ of the covariance of ( ) , Y x ω (denoted by ( ) 1 2
,

Y C x x ). ( ) ( ) ( ) 1 1 2 2 2 2 2 2 1 1 , , , ; , y x dy dx y x y x y x C i i D i Y φ λ φ = (3)
In order to reduce the computational effort in solving Eq. ( 3), the eigenvalues and eigenfunctions are assumed in the form of

( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 2 1 2
, ,

i i i k i i x y x y λ λ λ φ φ φ = =
, and the covariance is chosen as a parametric model separable in two directions:

( )

1 2 1 2 2 1 1 2 2
; exp

Y Y x y x x y y C x , y x , y l l σ - - = - - (4)
with the correlation length ,

x y l l (in x, y direction) depending on the fabric architecture.

Using the architectural parameters of reinforcement (e.g. inter-tow gap width for directional fabrics, porosity for random mat, etc.) acquired from image analysis or other measurements, realizations of permeability field can be obtained by means of numerical or analytical methods based on representative unit cell of the fabric structure. Then, optimization algorithms [START_REF] Li | Statistical characterization and robust design of RTM processes[END_REF][START_REF] Skordos | Stochastic simulation of woven composites forming[END_REF] can be applied to fit the experimental realizations into Eq. (4) to estimate ,

x y l l . While this procedure is under progress, a trial value for the correlation length is used in the numerical examples.

Once Eq. ( 3) is reduced to a set of 1-D equations for each direction, it can be transformed by means of the Galerkin method into generalized eigenvalue problems, from which the set { }

( ) ( ) ( ) 1 , exp M Y i i i i K x x ω µ λ φ ξ ω = = + (5)
As the random dimension of the input of the stochastic system is reduced to M, any random output, as non-linear function of the input, should be represented by the same random basis { } i ξ . Without knowing the marginal PDF, the Polynomial Chaos Expansion (PCE) [START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF] is employed to represent various outputs. For example, the pressure field is expressed in the PCE as

( ) ( ) ( ) ( ) 0 , P j j j p x p x ω ψ ξ ω = = (6)
where …,M) following the joint Gaussian PDF [START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF], and ( ) j p x denotes unknown deterministic coefficients. The order of the PCE (p) is taken as the highest order of the Hermite polynomials, related to the number of terms by P+1=(p+M)!/p!M! [START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF]. Once the deterministic coefficients are solved, the statistics of the target random response can be easily obtained by performing the MCS on its PCE.

( ) { } ξ ψ j (j = 0, 1, …, P) are M-dimensional Hermite polynomials in terms of { } i ξ ξ = (i = 1, 2,

Stochastic model for mold filling process

As a first step, only permeability uncertainty and the isothermal processes are considered (other random inputs or process models are ready to be incorporated in the presented stochastic framework). Following the well-established deterministic simulation models for RTM process [START_REF] Padmanabhan | Stochastic modeling of nonisothermal flow during resin transfer molding[END_REF][START_REF] Lam | Numerical simulation of the mould-filling process in resin-transfer moulding[END_REF], resin flow through fibrous media is approximated by a quasi-steady process. Each intermediate state is governed by Darcy's law and the continuity equation of an incompressible fluid with relevant random quantities represented by the suitable expansions:

( ) ( ) ( ) ( ) 1 0 exp M P Y i i i j j i j x p x Q x µ λ φ ξ ψ ξ = = ∇ ⋅ + ⋅ ∇ = (7)
where ( )

x Q
is the source term (which can also be random if needed). In order to simulate the moving resin flow front, most commercial codes employ the Control Volume / Finite Element Method (CV/FEM) [START_REF] Lam | Numerical simulation of the mould-filling process in resin-transfer moulding[END_REF], in which the flow front position is represented by the filling fraction in the control volumes. Due to the dependence of the flow front shape on the discretization of the computational domain, it is not efficient in simulating distinct and accurate flow fronts, which are required for the current study.

As an alternative, the Level Set Method (LSM) [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF] motivated by the multi-phase flow problems [START_REF] Tornberg | A finite element based level set method for multiphase flow[END_REF] begins to be adopted for LCM simulation [START_REF] Soukane | Application of the level set method to the simulation of resin transfer molding[END_REF][START_REF] Gantois | Numerical simulation of resin transfer molding using BEM and level set method[END_REF]. The basic idea of this technique is to represent the moving interface by a continuous function ( )

( ) ( ) ( ) 0 , , , = ∇ ⋅ + ∂ ∂ ω φ ω ω φ x x u t x ( 8
)
with the Darcy velocity derived from the expansions of permeability and pressure

( ) ( ) ( ) ( ) ( ) ( ) 1 0 1 , exp M P Y i i i j j i i j u x x p x ω µ λ φ ξ ω ψ ξ ω µ = = = - + ⋅ ∇ (9)
Besides, the signed distance feature of the level set function can only be preserved by performing reinitialization (Eq. 10) at certain time intervals.

( ) ( )( ) ( ) ( ) = = - ∇ + ∂ ∂ ω φ ϕ ϕ ω φ τ ϕ , 0 , 0 1 , x x x sign ( 10
)
where τ denotes the artificial time variable. ( )

ω φ ,
x is then replaced by the solution ( )

x ϕ (also random actually) which is a signed distance function with the same zero-contour, i.e. flow front. The details of numerical solution for Eq. ( 8) and ( 10) can be referred to [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF][START_REF] Tornberg | A finite element based level set method for multiphase flow[END_REF].

Solution of the stochastic model

Before choosing a suitable solution technique for Eqs. ( 7)~ [START_REF] Verleye | Permeability prediction for the meso-macro coupling in the simulation of the impregnation stage of RTM[END_REF], it is noted that, for the stochastic modeling of RTM process, a particular difficulty consists in the characterization of the random resin/air interface in the computation domain, represented by the zero-level set. Due to the signed distance requirement for any valid realization of a random level set function ( )

, x
φ ω , it is impossible to assume ( )

, x
φ ω as a stationary random field. For this reason, application of the random level set function has only been tempted for the random boundaries with simple geometries [START_REF] Nouy | An extended stochastic finite element method for solving stochastic partial differential equations on random domains[END_REF], which is far from the current case of an irregular interface. However, the uncertainty propagation in the model necessitates the randomness in the flow front or the level set function. Therefore, a stochastic solution technique able to circumvent the direct representation of ( ) , x φ ω , i.e. a so-called non-intrusive stochastic method [20] which deals with the whole non-linear system of equations as a black-box, is needed.

The Probabilistic Collocation Method (PCM) [START_REF] Stefanou | The stochastic finite element method: Past, present and future[END_REF]21] is found to be an optimal choice.

Focusing on an interesting random output of the model, e.g. the fill time ( ) , T x ξ , the stochastic system (Eqs. ( 7)~( 10)) can be represented by an implicit function

( ) ( ) 1 , exp M Y i i i i T x f x ξ µ λ φ ξ = = + (11)
As mentioned in Section 2.1, the random output can be approximated by the PCE:

( ) ( ) ( ) 0 , P j j j T x T x ξ ψ ξ = ≈ (12)
Thus, the solution is interpreted as seeking the set of unknown coefficients

( ) { } j T x to
minimize the residual of Eq. ( 12) in the statistical sense. Following the notion of weighted residual method, the mean residual can be expressed as

( ) ( ) ( ) ( ) ( ) 0 , 0 P j j k j T x T x w p d ξ ξ ψ ξ ξ ξ ξ = - = (13
k k w ξ ξ δ ξ - =
).

The integration is then reduced to a linear system of equations of the unknown coefficients

( ) { } j T x ( ) ( ) ( ) 0 , P j j k k j T x T x ψ ξ ξ = = (14)
with the coefficient matrix Therefore, the computational effort for solving the stochastic system is reduced to that of solving the deterministic system for N times, which means the PCM has great advantage in terms of the efficiency and convenience in operation.

However, the accuracy in predicting the response statistics is the primary factor in assessing the performance of the PCM. Since the zero residual is only ensured at the set of collocation points, they must be able to capture the high probability region in the sample space. In analogous to the Gaussian quadrature, the collocation points are usually chosen from the roots of the polynomials of one order higher than the PCE, and kept as close to the origin as possible in sample space [START_REF] Stefanou | The stochastic finite element method: Past, present and future[END_REF]21]. For example, if the 2 nd -order PCM is employed, the optimal collocation points should have the coordinates (in sample space) 0 or 3 ± (as roots of ( )

3 3 3 H ξ ξ ξ = -).
Up to now, the overall framework (implemented in MATLAB programming language) for uncertainty quantification in the RTM process model is completed based on the PCM and the LSM, as illustrated in the flowchart in figure 1.

Numerical examples

The effectiveness of the approach developed above for the uncertainty analysis is demonstrated by applying it to a simple two-dimensional mold filling model. The mold has a thickness of 5mm, a single inlet/outlet and an internal insert (figure 2(a)), containing an isotropic fabric of 30% fiber volume fraction. The fibrous permeability, assumed as a stationary log-normal random field with the mean value K 0 = 1×10 -9 m 2 , a coefficient of deviation CV=1.0 (to illustrate the robustness of the method) and the covariance function described by Eq. ( 4) with a correlation length equaling half the mold size, is taken as the only source of input uncertainty of the system. The and also the MCS with 10 3 realizations for validation.

Constant flow rate driven injection

First, the injection case using a constant flow rate (5×10 -6 m 3 /s) is considered. In order to design the appropriate stiffness and closing force of the mold for the targeted thickness or fiber volume fraction [START_REF] Comas-Cardona | A generic mixed FE-SPH method to address hydro-mechanical coupling in LCM processes[END_REF][START_REF] Walbran | Measurements of normal stress distributions experienced by rigid liquid composite moulding tools[END_REF], the range of resin pressure magnitude needs to be known. However, without taking into account the uncertainty in pressure field, the maximal pressure magnitude may be underestimated by traditional simulation codes, which increases the risk of preform deformation and mold deflection.

The capability of the current method is first validated by the converged solution of the MCS, taking the peak value of the inlet pressure as the random response of interest. As shown in figure 3(a), the converged mean value can be taken as obtained by the MCS with more than 300 realizations, while the PCM only needs a cost of 5% of the former to get the same result with negligible error. To get more detailed probabilistic information, the Cumulative Density Function (CDF) of the random response can be employed.

Constant pressure driven injection

As an alternative, the constant pressure injection is widely adopted in practice because of its low cost of injection device and the steady pressure level. However, the total fill time has to be controlled respecting the resin reactivity and process automation. Some techniques aiming at the prediction of accurate fill time are tempted [START_REF] Deleglise | Solution to filling time prediction issues for constant pressure driven injection in RTM[END_REF], but none focus on the random effects. Using a test example with inlet pressure given as 0.2MPa, the influence of the input uncertainty on the statistics of the fill time is analyzed by means of the current approach.

In analogy to the precedent injection case, the convergence for the typical statistics of the mold fill time is first studied. For the mean value, the MCS requires at least 600 realizations to obtain the steady value, which can be predicted by the 2 nd -order PCM with a relative error less than 2% using about 2.5% computational effort (figure 6(a)).

Comparing the CDF curves in figure 6(b), the PCM result can match the overall shape and capture the tail region of the MCS estimation using 10 3 realizations (approximately taken as converged).

In addition to the overall fill time for the mold, the randomness in the flow pattern can be easily investigated using the current method. Still using the converged MCS solution

  . The log-normal random field is accordingly represented in terms of discrete random variables:

  e. level set function) taking opposite signs in different media (impregnated or dry) and zero on their interface (resin-air flow front). In the current context, the level set function is a random variable ( ) ω φ ,xrepresenting the signed distance from location x to the flow front with a random location. Following the classical LSM[START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF], the movement of resin flow front during a certain time step is characterized using the advection equation[START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF]

  stochastic solution technique. For the instance of the PCM, the zero residual is enforced at a set of collocation points { } k ξ (k = 1, 2, …, N) in sample space (equivalent to using the Dirac-functions at the collocation points as weight functions ( ) ( )

ξ

  by means of the PCE and the deterministic system of equations, respectively.

  4 th -order K-L expansion is performed on the permeability field, from which the generated realizations can reproduce the input statistics (e.g. the local PDF and the autocorrelation) with acceptable accuracy, as shown in figures 2 (a)~(b). Resin with viscosity of 0.1 Pa sis injected into the mold by means of two different schemes respectively. In each case, the deterministic solution resulting from a homogeneous fabric with permeability K 0 is used to normalize the mean of the stochastic solution to highlight their discrepancy. The stochastic simulation is performed (on a CPU of Intel Xeon X5550/2.67GHz) using the 2 nd -order PCM (with 15 collocation points selected following the rules in Section 2.3)

Figure 3 (

 3 b) compares the CDF estimated by the 2 nd -order PCM to those by the MCS with 10 2 or 10 3 realizations, showing that the convergence of the overall CDF can be accurately reached by the PCM (except the tail region of CDF>0.87) at a cost less than 3% of that of the MCS.Using the convergent stochastic solution, the evolution of the mean and the CV of the inlet pressure during the injection process are displayed in figures 4(a)~(b), respectively, showing the consistency of the PCM estimation with the MCS standard through the time history. Results show that the long-term response (i.e. the average response of a large sample) is underestimated by the deterministic simulation, and the discrepancy may be significant for an input with high variance, as figure4(b)shows that the CV of the response (around 0.65~0.70) is comparable to that of the input. Besides, the PCM prediction of the CV of inlet pressure increases as flow front advances, demonstrating the uncertainty propagation in the model. In addition, the peak in the mean value between 300 and 350 s with a magnitude comparable to the maximal value reflects the disturbance of the internal insert (referring to the isochrones in figure2(a)). Due to the simple configuration of the mold, the spatial distribution of the mean and CV of the random pressure field are both analogous to that of the deterministic pressure field and not displayed here. Results show that the CV of pressure at any location (except those close to outlet) is largely in the range of 0.6~0.7, which confirms the importance of the pressure variability for mold design.The irregular flow pattern resulting from the random permeability field is the main obstacle for mold design. As a basis of the optimization schemes, the simulation tool should be able to assess the impact of the random input on the flow pattern. The resin-arriving time (referred to as fill time below) at any location is taken as the random variable quantifying the uncertainty in flow front. First, the PCM prediction for the mean and CV of fill time is compared to the respective MCS solution and validated (both have a slight dependence on discretization but only at a local scale). As can be seen in figures 5(a)(b), the mean value has small fluctuation around the deterministic one, implying a minor fluctuation in the flow pattern. Figures5(c)(d)show that the CV of the fill time has a maximum of around 0.15, with relatively higher values in the regions close to the corners or affected by the insert. However, the uncertainty in permeability field has no major influence on the flow pattern for this injection scheme.

  for validation, the mean value of the local fill time (normalized by the deterministic value) is displayed in figures 7(a)(b) for both methods. Simulation results show that at any location, the mean of random fill time is longer than the deterministic prediction by about 15%~30%. The two methods give similar prediction in terms of magnitude and spatial distribution. The discrepancy between the stochastic and the deterministic solutions demonstrates the importance of the random effects for choosing suitable process parameters to account for the variation in the material properties.The uncertainty propagation can be observed from the CV of the fill time (figure7(c)(d))which increases as the flow front moving forward, reaching about 0.8 (comparable to the CV of the input random field) near the outlet. Comparing to the results of constant flow rate injection (figure5), the constant pressure injection appears more sensitive to the input uncertainty, in terms of the variability in the flow pattern and in mold fill time.By comparison to the MCS (10 3 realizations) estimation, the PCM is able to predict accurately the statistical moments of the fill time which is an implicit function of the permeability, and its capability of randomness propagation is proved.4 ConclusionsA stochastic simulation code has been developed by combining the Level Set technique with the Probabilistic Collocation Method, taking the advantage of efficiency in both the moving interface simulation and uncertainty propagation, and thus suitable to be embedded as a module into optimization or sensitivity analysis algorithms[START_REF] Sudret | Comparison of finite element reliability methods[END_REF], instead of the time-consuming MCS. Numerical examples illustrate the capability of the developed simulation code in uncertainty quantification in the mold filling problem, as well as its great advantage over the MCS in terms of computational efficiency (e.g. for the given numerical tests, the CPU cost of the PCM is on average less than 5% of that of the MCS) for predicting the response statistics with comparable accuracy. Stochastic simulation results of simple mold filling processes reveal distinct discrepancies between the stochastic and deterministic solutions for various responses, which indicates the necessity to take into account the intrinsic randomness in transport properties.In spite of the simple numerical tests shown in this paper, potential application in many aspects can be developed based on the current framework, such as the simulation of RTM process with multiple sources of input uncertainty (e.g. the race-tracking permeability, defects distributed in fabric, randomness in resin viscosity and curing), or the optimization of process parameters to reduce the sensitivity of outputs to the input uncertainty.

Figures

  Figures

Fig. 1

 1 Fig. 1 Flowchart of simulation procedure

  Fig. 1 Flowchart of simulation procedure

Fig. 4 Fig. 5

 45 Time history of statistical moments of inlet pressure: (a) Mean (normalized by deterministic value); (b) CV. Statistical moments of fill time (constant flow rate injection) comparing between MCS (10 3 realizations) and PCM (2 nd -order) estimation: (a) MCS estimation of mean value (normalized by deterministic value); (b) PCM estimation of mean value (normalized by deterministic value); (c) MCS estimation of CV; (d) PCM estimation of CV.

Fig. 6 Fig. 7

 67 Convergence of statistics of mold fill time (constant pressure injection): (a) Mean (normalized by deterministic value); (b) CDF. Statistical moments of fill time (constant pressure injection) comparing between MCS (10 3 realizations) and PCM (2 nd -order) estimation: (a) MCS estimation of mean value (normalized by deterministic value); (b) PCM estimation of mean value (normalized by deterministic value); (c) MCS estimation of CV; (d) PCM estimation of CV.

Acknowledgements:

The authors would like to acknowledge the Conseil Régional du Nord -Pas de Calais (France) for the financial support for this study.


Olivier Teytaud Tao
email: olivier.teytaud@inria.fr

Comparison-Based Complexity of Multiobjective Optimization

Several comparison-based complexity results have been published recently, including multi-objective optimization. However, these results are, in the multiobjective case, quite pessimistic, due to the huge family of fitness functions considered. Combining assumptions on fitness functions and traditional comparison-based assumptions, we get more realistic bounds emphasizing the importance of reducing the number of conflicting objectives for reducing the runtime of multiobjective optimization. The approach can in particular predict lower bounds on the computation time, depending on the type of requested convergence: pointwise, or to the whole Pareto set. Also, a new (untested yet) algorithm is proposed for approximating the whole Pareto set.

Introduction

The comparison-based model [START_REF] Teytaud | General lower bounds for evolutionary algorithms[END_REF][START_REF] Fournier | Lower bounds for comparison based evolution strategies using VC-dimension and sign patterns[END_REF] is a recent tool for modeling evolutionary computation. Most evolutionary algorithms are comparison-based, and it is known that comparison-based black-box complexity is equivalent to black-box complexity for some invariant families of fitness functions [START_REF] Gelly | Comparison-based algorithms are robust and randomized algorithms are anytime[END_REF].

[20] has shown lower complexity bounds for multiobjective optimization. Unfortunately, the results therein, based on a worst case analysis on very wide families of objective functions, are quite pessimistic and exponential in the number of objectives. [START_REF] Rudolph | Convergence rates of sms-emoa on continuous bi-objective problem classes[END_REF] has shown a runtime analysis for multiobjective optimization, using the results in [START_REF] Witt | Rigorous runtime analysis of a (mu+1) es for the sphere function[END_REF] on the mono-objective case, with much more realistic rates, thanks to assumptions on fitness functions, supposed to be quadratic. We will here keep their assumptions on families of functions; basically, their model is a natural extension of the classical sphere function to the multiobjective case. This paper extends results from [START_REF] Rudolph | Convergence rates of sms-emoa on continuous bi-objective problem classes[END_REF], limited to small values of the offspring size λ and the number d of objective functions, by using comparison-based complexity results from [START_REF] Teytaud | General lower bounds for evolutionary algorithms[END_REF]. Importantly, the results emphasize the importance of the number of conflicting objectives in the runtime analysis.

Section 2 introduces comparison-based complexity. Section 3 recalls fundamental lower bounds in this setting (basically, the convergence is at best linear in log-scale, and the constant decays linearly as a function of the dimension). Section 4 provides lower bounds for the multiobjective case: essentially, there is, for convergence to the whole Pareto set, an additional factor linear as a function of the number of conflicting objectives. Section 5 shows corresponding upper bounds on the comparison-based complexity. Section 6 summarizes the paper and discusses limitations and further work.

2 Comparison-based algorithms of type (µ + , λ)

This section is devoted to a formal definition of comparison-based multiobjective optimization algorithms. The maths are strongly based on [START_REF] Teytaud | General lower bounds for evolutionary algorithms[END_REF][START_REF] Teytaud | Lower bounds for evolution strategies using vc-dimension[END_REF][START_REF] Fournier | Lower bounds for comparison based evolution strategies using VC-dimension and sign patterns[END_REF], but the application to multi-objective optimization is new and provides interesting specific hints.

We here follow the lines of [START_REF] Fournier | Lower bounds for comparison based evolution strategies using VC-dimension and sign patterns[END_REF] for defining four classes of evolutionary algorithms: Selection-Based Non-Elitist Strategies, Selection-Based Elitist Strategies, Full-ranking Non-Elitist Strategies, Full-Ranking Elitist Strategies. Let λ and µ be two integers. The first case is the case of Selection Based evolution strategies (SB-µ, λ-ES). In this case, there is a set I of internal states. The algorithm starts in the initial state I 0 returned by the function initial state. At each iteration, the algorithm follows these three successive steps. First generate a set of λ points, called the offspring. Then select only µ of these points, by any procedure, as a function of the fitness values. Finally, the internal state is updated. In the mono-objective case, this includes for example cumulative step-length adaptation (CSA) by [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF], covariance matrix adaptation (CMA) by [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF], mutative self-adaptation (SA) by [START_REF] Rechenberg | Evolutionstrategie: Optimierung Technischer Systeme nach Prinzipien des Biologischen Evolution[END_REF][START_REF] Schwefel | Adaptive Mechanismen in der biologischen Evolution und ihr Einfluss auf die Evolutionsgeschwindigkeit[END_REF], covariance matrix adaptation by SA (CMSA) by [START_REF] Beyer | Covariance matrix adaptation revisitedthe CMSA evolution strategy[END_REF], the one-fifth rule for step-size adaptation [START_REF] Rechenberg | Evolutionstrategie: Optimierung Technischer Systeme nach Prinzipien des Biologischen Evolution[END_REF][START_REF] Beyer | The Theory of Evolution Strategies[END_REF]. In the multi-objective case, it includes NSGA-II [START_REF] Deb | A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II[END_REF], SPEA2 [START_REF] Zitzler | Spea2: Improving the strength pareto evolutionary algorithm[END_REF]. Techniques which use gradients or surrogate models are not concerned.

SB-(µ + λ)-ES (the second case, termed elitist) selects µ points among the offspring and the points selected at the previous step.

Please notice that the framework is very general: it includes adaptation by cross-entropy methods (CEM by [START_REF] De Boer | A tutorial on the cross-entropy method[END_REF]) or Estimation of Distribution Algorithms like e.g. UMDA [START_REF] Mühlenbein | Evolutionary computation and Wright's equation[END_REF], Compact Genetic Algorithm [START_REF] Harik | The compact genetic algorithm[END_REF], Population-Based Incremental Learning [START_REF] Baluja | Population-based incremental learning: A method for integrating genetic search based function optimization and competitive learning[END_REF], Relative Entropy [START_REF] Mühlenbein | The estimation of distributions and the minimum relative entropy principle[END_REF] and Estimation of Multivariate Normal Algorithms (EMNA by [START_REF] Larranaga | Estimation of Distribution Algorithms[END_REF]), Direct Search Methods [START_REF] Conn | Recent progress in unconstrained nonlinear optimization without derivatives[END_REF], the Nelder-Mead algorithm [START_REF] Nelder | A simplex method for function minimization[END_REF], the Hooke&Jeeves algorithm [START_REF] Hooke | Direct search" solution of numerical and statistical problems[END_REF] -and MOO techniques based on these algorithms, without surrogate model learning. This points out the links between all these algorithms.

We present SB-(µ, λ)-algorithms (resp. SB-(µ + λ)-algorithms) in Algo- ω,n ← proposal(In) end while rithm 1. In this algorithm (also in Algorithm 2), the concatenation of the two vectors x = (x 1 , . . . , x k ) and x ′ = (x ′ 1 , . . . , x ′ ℓ ) is denoted by x ⊕ x ′ = (x 1 , . . . , x k , x ′ 1 , . . . , x ′ ℓ ); we also use the shortcut v ∈ (x 1 , . . . , x k ) to express that there exists i ∈ [ [1, k]] such that x i = v.

Notice that the length of the vector v n is equal to µ except maybe during the first iterations of the algorithm in the case λ < µ: this explains the use of the auxiliary variable ℓ in Algorithm 1 for the (µ + λ)-case.

Finally, we would like to explain a generalization of SB-(µ + , λ)-ES, called Full Ranking (µ + , λ)-ES (FR-(µ + , λ)-ES). Instead of just giving the best µ points (i.e., the µ points with the lowest fitness values), we can consider a selection procedure which returns the best µ points ordered with respect to their fitness.

The outline of these algorithms is given in Algorithm 2. More precisely, the selection step described in this algorithm works as follows. Given the vector of points E n = (z 1 , . . . , z p ) considered at step n, the function select returns a vector of µ distinct integers v n = (i 1 , . . . , i ℓ ), but as the full ranking matters, there's no constraint on v n ; it's just a vector of ℓ distinct points in [ [1, p]] (Once again, the length of the vector v n may not be equal to µ at the beginning of the algorithm in the "+" case).

Note that both Algorithms 1 and 2 define a class of algorithms: in order to obtain an algorithm, one has to specify the distribution of ω, how the offspring is generated (function generate), the space of states I as well as the functions initial state and update, and finally the function proposal. Throughout the chapter, Algorithm 2 Full Ranking (µ, λ)-ES (resp. Full Ranking (µ + λ)-ES). Framework for evolution strategies based on full ranking, working on a fitness function f . The random variable ω is a random seed. Compared to Algorithm 1, the vector of integers v n obtained at the selection step is now ordered with respect to the fitness values of points from E n ; this framework is thus more general.

Initialization: I0 ← initial state(ω); S0 ← ∅; n ← 0 while true do n ← n + 1 Generation step (generate an offspring On of λ distinct points): On ← generate(In-1) En ← On (resp. En ← On ⊕ Sn-1) ℓ ← min(µ, |En|) Selection step: vn ← select(En, f ), i.e.:
The vector vn = (i1, . . . , i ℓ ) is the vector of indices of selected points (i.e. an increasing sequence of length ℓ in [ [START_REF] Baluja | Population-based incremental learning: A method for integrating genetic search based function optimization and competitive learning[END_REF]min(µ,|En|)]] Update the internal state:

In ← update(In-1, En, vn) Sn ← (En,i 1 , . . . , En,i ℓ ) x (f ) ω,n ← proposal(In) end while
we assume that all functions involved in these algorithms are measurable.

Importantly, in spite of the formulation in Algorithms 1 and 2, the formalization is not restricted to deterministic algorithms. In Algorithms 1 and 2, all steps are deterministic, but there is an initial source of randomization, ω, which can be as large as required, e.g. it might be an infinite sequence of realizations of a random Gaussian variables or random uniform variables. So this is not a restriction, the algorithm must just be rephrased accordingly, so that ω is reported in the internal state and used in e.g. the "proposal" function.

Complexity lower bounds

We consider a bounded domain E and a distance dist : E × E → R on E or on a superset of E. For ε > 0, we define N E (ε) to be the maximum integer n such that there exist n distinct points x 1 , . . . , x n ∈ D with dist(x i , x j ) 2ε for all i = j. N E (ε) is termed the packing number, and quantifies the "size" and "dimensionality" of the domain as explained below. In particular, N E (ε) = |E| when ε is small enough in the case of a finite domain D, and log N E (ε) ∼ z log(1/ε) when ε → 0 if the domain E is bounded with non-empty interior [START_REF] Kolmogorov | ε-entropy and ε-capacity of sets in functional spaces[END_REF] of dimension z. Please notice that these definitions are consistent both for continuous domain E ⊂ R z and discrete domain E ⊂ R z .

We denote by x * (f ) ⊂ D the set of optima. It can be a point, in particular in the mono-objective case, or a multi-objective algorithm in which the Pareto set is restricted to one point; but, in the general case, we get a set of optima, a Pareto set of dimension typically equal to the number of objectives minus one. For keeping generality in our results, we will need a metric on subsets of D.

Typically, x (f ) ω,n is the proposal by the multi-objective algorithm -it can be an approximation of the target Pareto set, or it can be just a point, if the algorithm is intended to propose a pointwise convergence. Then, dist(x (f ) ω,n , x * (f )) can be:

• inf x∈x * (f ) ||x -x (f ) ω,n || for a pointwise convergence (i.e. x (f ) ω,n is a point). • The Hausdorff distance between x (f ) ω,n and x * (f ) (i.e. x (f ) ω,n is a set).
The precise definition of the metric does not matter for the moment, and we'll discuss these examples later.

For any given optimization algorithm as in Algorithm 2, and for ε > 0 and δ > 0, we define n ε,δ be the minimum number n of iterations such that with probability at least 1 -δ, an optimum is found at the n-th iteration within distance ε; i.e., n ε,δ is minimal such that for all n n ε,δ and for all f ∈ F, [START_REF] Teytaud | General lower bounds for evolutionary computation[END_REF] is the branching factor K of algorithms as defined above.

P r ω [dist(x (f ) ω,n , x * (f )) < ε] 1 -δ. A crucial notion from
K = sup E |{select(E, f ) : f ∈ F}|,
where the supremum holds for:

• E any set of size at most λ in the case of SB-(µ, λ)-ES or Full Ranking (µ, λ)-ES;

• E any set of size at most µ+λ in the case of SB-(µ+λ)-ES or Full Ranking (µ + λ)-ES.

Various bounds on K can be derived [START_REF] Teytaud | General lower bounds for evolutionary algorithms[END_REF]:

K µ λ 2 λ / √ 2πλ for SB -(µ, λ) -ES; (1)
K µ (λ + µ) for SB -(µ + λ) -ES; (2)
K λ!/(λ -µ)! for F R -(µ, λ) -ES; (3)
K (λ + µ)!/λ! for F R -(µ + λ) -ES. (4)
Theorem 1 (Lower bound hitting time for (µ + , λ)-ES.) Consider a (µ, λ)-ES or (µ + λ)-ES, as defined in Algorithm 2, and a set F of fitness functions on domain D, i.e. F is a set of functions from D to R d , i.e. F ⊂ (R d ) D , and define x * (f ) the optimum of f (this is a set, possibly restricted to a point), and such that {x * (f ) : f ∈ F} = D. Let ε > 0 and δ ∈]0, 1[. In particular, if K denotes the branching factor of the algorithm, then

n ε,δ log(1 -δ) log K + log N (ε) log K .
In the next section, we will see the consequences of this theorem in multiobjective optimization.

Results in multiobjective optimization

The section above, based on earlier results in comparison-based optimization and inspired by communication complexity, provided general results which can be used in many frameworks around comparison-based algorithms, and also around iterative algorithms with bounded branching factor. We here specialize the analysis in the case of multi-objective optimization.

Section 4.1 will provide necessary tools in terms of packing numbers for later sections. Section 4.2 will provide general results in terms of computation time before reaching a given precision, and section 4.3 will provide the main results, namely the number of fitness evaluations before reaching a given precision.

Packing numbers

Packing numbers go back to [START_REF] Kolmogorov | ε-entropy and ε-capacity of sets in functional spaces[END_REF]. A basic case is N E (ε) = Θ(1/ε z ) for E an open bounded subset of R z ; many cases are derived from this one. We assume in the rest of this paper that N d, and we work on a search space D = [0, 1] N .

Standard case. If E is equal to [0, 1] N and the distance is the standard euclidean metric, then we get N E (ε) = Θ(1/ε N ). This is widely used in monoobjective optimization [START_REF] Teytaud | General lower bounds for evolutionary algorithms[END_REF][START_REF] Teytaud | On the hardness of offline multi-objective optimization[END_REF]. We now switch to more important cases (more important for the multiobjective case), namely cases involving sets (which will be, for us, Pareto sets).

Sphere functions. If E is the set of possible Pareto sets in [0, 1] N defined by d spherical objective functions with optima in [0, 1] N , i.e.

E = {convexHull(x 1 , . . . , x d ); (x 1 , . . . , x d ) ∈ [0, 1] N d } (5)
and if E is equipped with the Hausdorff metric, then

N E (ε) = Θ(1/ε N d ). (6)
i.e. we have the same set of possible Pareto sets as above, but not the same distance), then the situation is more complicated: the algorithm is only allowed to provide a point, whereas the target is a set (the Pareto set); this dissymetry implies, for a branching factor theorem, that we need a slightly adapted definition of N E (ε), namely N E (ε) is the minimum number of points x 1 , . . . , x N (ε) in D such that for all e ∈ E, there is i such that d(x i , e) ε. This is a covering number; covering numbers can be used as well as packing numbers in branching factor theorems.

Furthermore, we have to distinguish different cases. Converging to just one point of the Pareto set can be made by optimizing just one objective function (provided that there is an objective function different from all others), and we get the same convergence as for mono-objective functions. A more interesting question is what happens when the objective functions are constrained to be different. We therefore define the set E as follows

E = {intersections of af f ine subspaces of R N with [0, 1] N } (8)
which means that the optima of the d sphererical objective functions are sufficiently far away (outside the domain possibly -but obviously we request convergence only to the Pareto set in the domain), so that the Pareto set has its frontiers on the frontiers of the domain D = [0, 1] d . We point out that considering smaller affine subspaces (e.g. included in [a, b] N ) would be sufficient for our results (both upper bounds and lower bounds, in Table 1); we use E as above for simplicity. Then, we claim: Lemma 2.Consider E as defined in Eq. 8. Then, for the pointwise metric (Eq. [START_REF] Fournier | Lower bounds for comparison based evolution strategies using VC-dimension and sign patterns[END_REF],

N E (ε) = Ω(1/ε N +1-d ). (9)
Proof: Consider the set

E ′ = {(x 1 , . . . , x d-1 , 0, . . . , 0); (x 1 , . . . , x d-1 ) ∈ R d-1 } e ′ is isomorphic to [0, 1] N +1-d and therefore N E ′ (ε) = Ω( 1 ε N +1-d ). As E ′ ⊂ E, N E ′ (ε) N E (ε) and therefore N E (ε) = Ω(1/ε N +1-d ).

Lower bound on the number of comparisons

The number of comparisons performed during an optimization run is always a lower bound on the computational cost; if the cost of the fitness function is low and if the internal cost of the optimization algorithm is low, then this becomes the main computational cost. There are other cases in which the number of comparisons is a very natural criterion: when it is the only available information, as e.g. when parametrizing a strategy for two-players games. The following result is an immediate consequence of Theorem 1; note that we have a 3 in the equations because the branching factor, for each comparison, is 3, as two points might lead to the dominance of the first by the second, or a dominance of the second by the first, or they might be uncomparable; there are 3 cases.

Theorem 2 (W.r.t nb of comparisons) Consider assumptions of Theorem 1 with K = 3 corresponding to the result of a comparison between the fitnesses of two visited points (the three outcomes are uncomparable, better than, worst than; we are in a multiobjective setting). Then, with log 3 (x) = log(x)/ log(3), the number of comparisons n c required for ensuring with probability 1 -δ a precision ε is n c log 3 (1 -δ) + log 3 (N (ε)). I.e., formally, P (

x nc -x * (f ) < ε) 1 -δ ⇒ n c log 3 (1 -δ) + log 3 (N (ε)).
This very general result is relevant when the cost of the comparisons is not negligible. In many cases, the cost of the comparisons is negligible in front of the fitness evaluations -in next section, we will switch to complexity measures based on the number of fitness evaluations.

Lower bound on the number of function evaluations

We can now derive the main bounds. Table 1 applies Theorem 1 using Eq. 9 for the first row, and Eq. 6 for the second row, and Eqs 1-4 for bounding K.

Category SB- SB- FR- FR- of ES (µ, λ) (µ + λ) (µ, λ) (µ + λ) Pointwise convergence N +1-d log( µ λ ) N +1-d log( µ λ+µ ) N +1-d log( λ! (λ-µ)! ) N +1-d log( (λ+µ)! λ! ) to the Pareto set Hausdorff convergence N d log( µ λ ) N d log( µ λ+µ ) N d log( λ! (λ-µ)! ) N d log( (λ+µ)! λ!
) to the Pareto set Table 1: Each equation is the coefficient of log(ε), within constant multiplicative factor, in the lower bound on the runtime; see Th. 1 for the complete lower bound (e.g. for pointwise convergence in the SB-(µ, λ)-ES case, we get a lower bound n ε,δ = Ω log (1-δ) log K + N +1-d log( µ λ ) log ε ). N is the dimension of the search space and d is the number of objective functions; we assume N d.

Upper bounds on the comparison-based complexity

This section is devoted to the tightness of some dependencies in the bounds above.

We have shown (Theorem 1 and Eq. 6) that the runtime is lower bounded by Θ(N d log(1/ε)) for convergence to the whole Pareto set for sphere functions. We have the following upper bound counterpart:

Proposition 1: For any set of d spherical objective functions as defined in Eq. 5 and Hausdorff metric, algorithm 3 outputs an approximate Pareto set with precision ε (for the Hausdorff metric) within runtime O( N d log( 1 ε ) of the Pareto set.

Remark: The constant in the O(.) does not depend on the chosen spherical functions.

Algorithm 3 Algorithm for approximating the whole Pareto set in the case of d spherical objective functions. We assume that each e i is an instance of a comparison-based algorithm with runtime N log(1/ε) in the monoobjective case (e.g. the Hooke& Jeeves pattern search algorithm), i.e. the number of fitness evaluations is O(N log(1/ε)) for reaching precision ε. The convergence properties of variants of the Hooke&Jeeves algorithm can be found in Proof: There are monoobjective optimization algorithms with runtime bounded by O(N log(1/ε)); see e.g. the proof in [START_REF] Fournier | Lower bounds for comparison based evolution strategies using VC-dimension and sign patterns[END_REF] for a variant of the Hoock&Jeeves algorithm.

Algorithm 3 applies such an algorithm to each of the d objective functions. After O(dN log(1/ε)) evaluations, we get d distinct optima within precision ε. The convex hull of these d points is, within precision O(ε), the Pareto set.

Proposition 2: For any set of d spherical objective functions as defined in Eq. 8 and the metric defined in Eq. 8, algorithm 4 outputs a point at Euclidean distance at most ε of the Pareto set within runtime O (N + 1 -d) log( 1 ε ) . Proof: Algorithm 4 proceeds as follows:

• It spends N + 1 function evaluations for finding a linear subspace of R N generated by N + 1 -d linearly independent vectors which are not parallel to the Pareto set.

• Then, it performs an optimization in this affine subspace by any comparison-based algorithm with runtime O dim(B) log

1 ε = O (N + 1 -d) log 1 ε .

Conclusions and discussion

Our results show that for convergence to the whole Pareto set, the number of conflicting objectives has a (multiplicative) linear impact on the runtime, as well as the dimension of the search space (Proposition 1 and complexity bounds in Table 1, second row). On the other hand, for pointwise convergence, the runtime is linear in N + 1 -d; a high number of conflicting objectives makes the Evaluate the fitness values at all points of the form (0, 0, . . . , 0, 1, 0, . . . , 0). If f (0, 0, . . . , 0, 1, 0, . . . , 0) = f (0, . . . , 0), then (0, 0, . . . , 0, 1, 0, . . . , 0) is not parallel to the Pareto set. Pick N + 1 -d such vectors (which are not parallel to the Pareto set). Let B be the vector space generated by these N + 1 -d vectors. while (true) do Perform one iteration of (e.g.) the Hoock&Jeeves algorithm in the vector space B on fitness function f 1 . Let x be the current iterate of this algorithm. end while problem easier. This can be used for choosing between pointwise convergence and global convergence.

These intuitively satisfactory results are essentially a consequence of the realistic model used in [START_REF] Rudolph | Convergence rates of sms-emoa on continuous bi-objective problem classes[END_REF].

There are several possible extensions:

• Further analysis could include an iterative model: monoobjective optimization for each optimum separately, and then successive choices by the user of weights of each fitness function.

• Extend the results to quadratic fitness functions and not only spherical fitness functions.

• Apply the VC-type results as in [START_REF] Fournier | Lower bounds for comparison based evolution strategies using VC-dimension and sign patterns[END_REF] for improving the results in the case of λ large.

Also, whenever Algorithm 3 has been proposed for theoretical purposes (proving bounds), the idea of finding good points and then extrapolating them into an approximate Pareto set might be a good idea in the real world as well, probably with more than d points however -experiment this is the most immediate further work.

  This is proved by considering the convex hulls of sets of size at most d (which are exactly the Pareto sets corresponding to d sphere functions). Sphere functions, pointwise convergence. If E is the set of possible Pareto sets in [0, 1] N defined by d spherical objective functions and equipped with d(a, b) = inf y∈b ||a -y|| (pointwise metric);

Algorithm 4

 4 Algorithm for finding one point close to the Pareto set. Inputs: f 1 , . . . , f d , which are d spherical objective functions on [0, 1] N . Output: at each iteration, a point x hopefully converging to the Pareto set. Perform a design of experiments as follows:

  Fournier et al. Inputs: f 1 , . . . , f d , which are d spherical objective functions on [0, 1] N . Initialize e 1 , . . . , e d , some linearly converging comparison-based algorithms for mono-objective optimization. while (true) do for i ∈ [[1, d]] do Perform one iteration of algorithm e i on objective function f i . Let xi be the approximation of arg min f i provided by e i . end for The approximate Pareto set is the convex hull of {x 1 , . . . , xd }. end while

Acknowledgments

O. Teytaud is grateful to the LISC laboratory, where he wrote a part of this document. We are also grateful to Gunter Rudolph, Dimo Brockhoff, Marc Schoenauer and Dagstuhl's seminar 10361 for fruitful discussions.