Mickaël Rivier

Pietro Marco Congedo

Surrogate-Assisted Bounding-Box Approach Applied to Constrained Multi-Objective Optimisation Under Uncertainty

Keywords: Multi-objective optimisation, Uncertainty-based optimisation, Error bounding boxes, Imprecise Pareto front, Surrogate-assisting strategy

This paper is devoted to tackling constrained multi-objective optimisation under uncertainty problems. A Surrogate-Assisted Bounding-Box approach (SABBa) is formulated here to deal with robustness and reliability measures, which can be computed with tunable and refinable fidelity. A Bounding-Box is defined as a multi-dimensional product of intervals centred on the estimated objectives and constraints that contains the true underlying values. The fidelity of these estimations can be tuned throughout the optimisation so as to reach high accuracy only on promising designs, which allows quick convergence toward the optimal area. In SABBa, this approach is supplemented with a Surrogate-Assisting (SA) strategy, which is very useful to reduce the overall computational cost. The adaptive refinement within the Bounding-Box approach is based on the computation of a Pareto Optimal Probability (POP) for each box. We first assess the proposed method on several analytical uncertainty-based optimisation testcases with respect to an a priori metamodel approach in terms of a probabilistic modified Hausdorff distance to the true Pareto optimal set. The method is then applied to two engineering applications: the design of two-bar truss in structural mechanics and the design of a thermal protection system for atmospheric reentry.

Le framework SABBa appliqué à l'optimisation contrainte multi-objectif sous incertitudes

Résumé : Ce travail se concentre sur les problèmes d'optimisation sous incertitudes multiobjectif. Le framework SABBa (Surrogate-Assisted Bounding-Box approach) est proposé pour traiter les mesures de robustesse et de fiabilité, qui peuvent être calculées avec une fidélité adaptable et raffinable. Une Boîte conservative est définie comme un produit multi-dimensionnel d'intervalles centrés sur les objectifs et contraintes estimés, de manière à contenir leur valeur réelle. La fidélité de ces estimations peut être ajustée au fur et à mesure du processus d'optimisation de manière à atteindre des estimations haute-fidélité uniquement sur les individus les plus prometteurs. Dans SABBa, cette approche est couplée à une stratégie d'assistance par modèle de substitution, qui permet de réduire le coût global du processus d'optimisation. Nous proposons de guider le raffinement adaptatif des Boîtes conservatives par le calcul de la probabilité pour chaque boîte d'être Pareto-optimale.

La méthode proposée est d'abord appliquée à plusieurs cas-tests analytiques et comparée à un approche par modèle de substitution a priori. Un indicateur par distance de Hausdorff modifiée probabiliste est proposée pour comparer les différentes stratégies. SABBa est ensuite appliqué sur deux cas ingénierie : une structure à deux poutres et un système de protection thermique ablatif pour réentrée atmosphérique.

Mots-clés : Optimisation multi-objectif, Optimisation sous incertitudes, Boîte d'erreur conservative, Front de Pareto imprécis, Assistance par modèle de substitution

Introduction

The topic of optimisation under uncertainty is of particular interest for companies nowadays as (i) the always increasing computational power now allows for multiple runs and reanalysis of simulations and (ii) robustness, reliability and cost optimality are critical factors for assessing the quality of a product and the efficiency of a company. In this context, two main areas are of primary interest: Robust Optimisation (RO) and Reliability-Based Optimisation (RBO). One of the central issues in conducting such analysis is the computational cost of attaining a prescribed level of accuracy. Most of the time, these optimisations are performed on numerical simulations, with run times ranging from minutes to weeks. This prohibits direct simulation-based methods such as nested Monte Carlo Simulations (MCS), performing a full Uncertainty Propagation (UP) at each optimisation iteration.

Contributions to reliability-based optimisation mostly rely on the computation of failure probability, such as the First-Order Reliability Method (FORM), and the Performance Measure Approach (PMA) [START_REF] Lee | A comparative study on reliabilityindex and target-performance-based probabilistic structural design optimization[END_REF][START_REF] Zhao | A comparison of deterministic, reliability-based topology optimization under uncertainties[END_REF]. Other decoupled approaches are based on local approximations of the failure probability and the associated gradient to perform line search optimisation, as in [START_REF] Jensen | Reliability-based optimization of stochastic systems using line search[END_REF], based on Subset Simulation (SS) computations, or in [START_REF] Valdebenito | Efficient strategies for reliability-based optimization involving non-linear, dynamical structures[END_REF], supplemented by a stabilisation strategy. The idea of using a surrogate model to tackle the computational cost limitation has been extensively used with Neural Networks [START_REF] Papadrakakis | Reliability-based structural optimization using neural networks and monte carlo simulation[END_REF], Polynomial Chaos Expansion (PCE) [START_REF] Keshavarzzadeh | Topology optimization under uncertainty via non-intrusive polynomial chaos expansion[END_REF] and PCE-Kriging metamodels [START_REF] Schöbi | Rare event estimation using polynomial-chaos kriging[END_REF]. Finally, in [START_REF] Bect | Sequential design of computer experiments for the estimation of a probability of failure[END_REF], a Stepwise Uncertainty Reduction (SUR) strategy for failure probability estimation is proposed, based on the volume of the excursion set. Other approaches to mention are Ref. [START_REF] Ching | Approximate optimization of systems with highdimensional uncertainties and multiple reliability constraints[END_REF] and [START_REF] Xu | Efficient reliability assessment of structural dynamic systems with unequal weighted quasi-monte carlo simulation[END_REF], based on "safety factor" constraints and Extreme Value Distribution (EVD) respectively.

In literature, there are also several contributions to the topic of Robust Optimisation, in particular applied to structural optimisation and Topological Optimisation (TO). Taguchi's robustness paradigm provides a typical RO formulation. It aims at maximising the performance (or mean performance) whilst minimising the associated variance. Refs. [START_REF] Amir | Efficient reanalysis techniques for robust topology optimization[END_REF][START_REF] Chen | Topology optimization of structures with interval random parameters[END_REF][START_REF] Doltsinis | Robust design of non-linear structures using optimization methods[END_REF][START_REF] Tang | Uncertainty based robust optimization method for drag minimization problems in aerodynamics[END_REF] are examples of this paradigm for RTO problems. To perform mean performance optimisation, Ref. [START_REF] Taflanidis | An efficient framework for optimal robust stochastic system design using stochastic simulation[END_REF] introduces the Stochastic Subset Optimisation (SSO), which is then supplemented by a local optimisation through Simultaneous Perturbation Stochastic Approximation (SPSA). SPSA was also used in [START_REF] Medina | Adaptive importance sampling for optimization under uncertainty problems[END_REF], coupled with an adaptive Importance Sampling (IS) strategy. Worst-case optimisation has been performed with monotonic assumptions in [START_REF] Elishakoff | Structural design under bounded uncertainty-optimization with anti-optimization[END_REF] and [START_REF] Mcwilliam | Anti-optimisation of uncertain structures using interval analysis[END_REF]. Surrogate-based methodologies are also massively exploited in the area of Robust Optimisation, for example with PCE [START_REF] Keshavarzzadeh | Topology optimization under uncertainty via non-intrusive polynomial chaos expansion[END_REF], and kriging [START_REF] Zhang | Sequential approximate optimization for design under uncertainty problems utilizing kriging metamodeling in augmented input space[END_REF]. A general discussion over the use of metamodels in the context of uncertainty-based optimisation is proposed in [START_REF] Jin | The use of metamodeling techniques for optimization under uncertainty[END_REF]. It shows that such methods yield promising results, notably the kriging or Gaussian Process (GP) surrogate models [START_REF] Lee | A global robust optimization using kriging based approximation model[END_REF][START_REF] Dellino | Robust optimization in simulation: Taguchi and krige combined[END_REF][START_REF] Eldred | Formulations for Surrogate-Based Optimization Under Uncertainty[END_REF]. Some works take full advantage of the nature of Gaussian processes. In [START_REF] Janusevskis | Simultaneous kriging-based estimation and optimization of mean response[END_REF], the projection of a GP surrogate model is conducted analytically to perform mean performance optimisation. Optimisation in the design space is then performed through Bayesian Optimisation (BO). In [START_REF] Le Riche | Gears design with shape uncertainties using controlled monte carlo simulations and kriging[END_REF], they exploit the ability of Gaussian process modelling to deal with heterogeneous noises in the observations. Finally, [START_REF] Binois | Quantifying uncertainty on pareto fronts with gaussian process conditional simulations[END_REF] uses conditional simulations in the context of multi-objective optimisation to sample possible Pareto optima. More generally, there have been extensions of GP-based optimisation (BO or Efficient Global Optimisation, EGO) to stochastic or noisy black-box evaluations, for example, with Sequential Kriging Optimisation (SKO) [START_REF] Huang | Global Optimization of Stochastic Black-Box Systems via Sequential Kriging Meta-Models[END_REF], and with the concept of Noisy Expected Improvement [START_REF] Picheny | Noisy Expected Improvement and on-line computation time allocation for the optimization of simulators with tunable fidelity[END_REF].

In several contributions [START_REF] Teich | Pareto-Front Exploration with Uncertain Objectives[END_REF][START_REF] Mlakar | Comparing solutions under uncertainty in multiobjective optimization[END_REF][START_REF] Fusi | An adaptive strategy on the error of the objective functions for uncertainty-based derivative-free optimization[END_REF], the authors introduced the concept of Bounding-Boxes and revisited the Pareto dominance rules. These concepts have been generalised in [START_REF] Rivier | Surrogate-Assisted Bounding-Box Approach for Optimization Problems with Approximated Objectives[END_REF] to tackle noisy estimations with tunable fidelity, and supplemented by a Surrogate-Assisting strategy to obtain the parsimonious Surrogate-Assisted Bounding-Box approach (SABBa) framework. The main contribution of the present paper is the extension of SABBa to the problem of optimisation under uncertainty, where both probabilistic objectives (RO) and probabilistic constraints (RBO) must be tackled. To this extent, we propose a Constrained Boxed Pareto Dominance rule to deal with estimated constraints. Parsimony is improved through the choice of GP surrogate models, allowing for natural computations of conservative errors and refinement criteria. In addition, the Pareto Optimal Probability (POP) is computed for each box to focus the refinements on boxes that are non-dominated with high probability. All numerical details for practical implementation and quantitative comparisons are given in the following. We show that the direct Double Loop approach, where a full Uncertainty Propagation (UP) is performed at each optimisation iteration, is extremely expensive and not competitive with respect to the so-called A Priori MetaModel (APMM) strategy. Hence, we systematically compare SABBa to the latest.

We depict a general uncertainty-based optimisation problem in Section 2, with some classical robustness and reliability measures. Then, we introduce the concept of Bounding-Boxes and Pareto Optimal Probability (POP) in Section 3. We extend the existing concepts from [START_REF] Teich | Pareto-Front Exploration with Uncertain Objectives[END_REF][START_REF] Mlakar | Comparing solutions under uncertainty in multiobjective optimization[END_REF][START_REF] Fusi | An adaptive strategy on the error of the objective functions for uncertainty-based derivative-free optimization[END_REF][START_REF] Rivier | Surrogate-Assisted Bounding-Box Approach for Optimization Problems with Approximated Objectives[END_REF] to constrained optimisation problems. The global algorithm is depicted in Section 4 alongside computational insights for measure estimation. Three variants of SABBa are quantitatively compared in Section 5 on analytical test-case and the most efficient approach is applied on two engineering applications in Section 6. Some conclusions and perspectives are drawn in Section 7.

Uncertainty-based optimisation: problem formulation

The objective of this section is to formulate a general expression for an Optimisation under Uncertainty problem including both Robust and Reliability-Based Optimisation (RO and RBO).

A deterministic constrained multi-objective optimisation problem can be described as follows:

minimise/maximise: f (x) , satisfying: g (x) ≤ 0, by changing: x ∈ X , where the m 1 objective functions are collected in a vector f ∈ R m1 , the m 2 constraint functions in g ∈ R m2 and x ∈ X are the n design variables included in the design domain X ⊂ R n .

In the presence of uncertain parameters ξ (e.g. environmental, material or geometrical parameters), the objective and constraint functions depend on both x and ξ. Hence, f (x, ξ) and g(x, ξ) must be replaced with robustness and reliability measures, ρ f (x) and ρ g (x) respectively. These allow the formulation of a general uncertainty-based constrained multi-objective optimisation problem, as follows: minimise/maximise: ρ f (x), satisfying:

ρ g (x) ≤ 0, by changing: x ∈ X , with ρ f ∈ R m1 and ρ g ∈ R m2 .
Many formulations are possible for ρ f and ρ g . The following are tackled in this paper:

Expectation: ρ f (x) = E ξ [f (x, ξ)],
Variance:

ρ f (x) = V ξ [f (x, ξ)],
Min/Max:

ρ f (x) = min ξ [f (x, ξ)] or max ξ [f (x, ξ)],
Quantile:

ρ f (x) = q p ξ [f (x, ξ)], p ∈ [0, 1]
, and the same for ρ g . Note that any other measure could also be treated similarly in the following.

Inria

Remark For conciseness, ρ f and ρ g are gathered in a vector ρ,

∀x ∈ X , ρ(x) = ρ f (x) ρ g (x) .
The main issue here is the accurate computation of these quantities. In this work, we approximate statistical measures through Monte Carlo Simulations (MCS) on surrogate models. To improve the accuracy of these approximations, one can increase the training set used for building the surrogate.

Hence, the uncertainty-based optimisation problem reduces to a constrained multi-objective optimisation problem with tunable fidelity computations for all objectives and constraints. The concept of tunable fidelity in the unconstrained case has been introduced in [START_REF] Rivier | Surrogate-Assisted Bounding-Box Approach for Optimization Problems with Approximated Objectives[END_REF] for multiobjective optimisation problems. Next section illustrates the concept of Bounding-Boxes and their applicability to constrained problems.

Bounding-Boxes

In this section, we introduce the concept of Bounding-Box. Notations and the Boxed Pareto dominance (introduced in [START_REF] Rivier | Surrogate-Assisted Bounding-Box Approach for Optimization Problems with Approximated Objectives[END_REF]) are recalled. The extension to constrained optimisation problems is then illustrated in Section 3.1.2. Finally, the concept of Pareto Optimal Probability is presented in Section 4.1 alongside computational insights.

Bounding-Box definition

A m-dimensional box is defined as follows, with a ∈ R m the center and r ∈ R m the positive half-width vector:

B(a, r) = b ∈ R m | b ∈ [a -r, a + r] ∈ ℘(R m ),
where ℘(R m ) is the power set of R m . The Bounding-Box approach consists in approximating an unknown value, here statistical measures ρ, by a conservative box containing the true values. More precisely, the unknown ρ is approximated by ρ l , l being the fidelity of the approximation. The exact approximation error ε l = ρρ l is conservatively approximated by ε l ≥ |ε l |. Hence it comes that ρ ∈ B( ρ l , ε l ) (see Fig. 1).

Boxed Pareto Dominance

To perform multi-objective optimisation, the Pareto dominance rules are usually exploited to assess dominance ( ), strict dominance ( ) or indifference (∼) between designs. In the Bounding-Box context, for comparing two boxes B(a, r), B(b, r ) ∈ ℘(R m ) 2 , the Boxed Pareto dominance is defined as follows:

• B(a, r) dominates B(b, r ) in the Boxed Pareto dominance sense: B(a, r) B B(b, r ) ⇐⇒ ∀j ∈ 1, m , ±a j + r j ≤ ±b j -r j and ∃j ∈ 1, m , ±a j + r j < ±b j -r j ;

• B(a, r) strictly dominates B(b, r ) in the Boxed Pareto dominance sense: ∀j ∈ 1, m , the symbol ± (implicitly ± j ) indicates the goal in the j th dimension:

B(a, r) B B(b, r ) ⇐⇒ ∀j ∈ 1, m , ±a j + r j < ±b j -r j ; RR n°9214 ρ ε l 1 ε l 2 ε l 1 ε l 2 ρ l B ρ l , ε l
± = + for minimisation -for maximisation.
Intuitively, B(a, r) dominates B(b, r ) if the worst outcome of B(a, r) dominates in the classical sense the best outcome of B(b, r ). An example is given hereafter in Fig. 2. In the case of bi-minimisation, the dominance are the following:

B 1 ∼ B B 2 , B 1 ∼ B B 3 , B 1 ∼ B B 4 , B 2 ∼ B B 3 , B 2 B B 4 , B 3 ∼ B B 4 .
Hence, only B 1 , B 2 and B 3 are non-dominated.

Constrained Boxed Pareto dominance

The Boxed Pareto dominance allows comparing boxes in a multi-objective optimisation context. For dealing with robust objectives as well as reliability-based constraints, the Pareto dominance should take into account constraint violation. To this extent, a Constrained Pareto dominance c is first proposed, where again refers to the classical Pareto dominance.

Let us assume that elements to compare are in R m and m = m 1 + m 2 with dimensions 1 to m 1 being the objective dimensions (subject to minimisation or maximisation) and the m 2 latest being the constraint dimensions (assumed as lower or upper bounds). We also introduce the vector a D , with D a set of indices, which is the projection of a on the associated dimensions. Formally, a D is defined as:

a D = d∈D a d e d ,
with e k the unit vector in the k-th direction.

In the multi-objective constrained problem, vector a 1,m1 represents the objective values and a m1+1,m the constraint values. Without loss of generality, these m 2 constraints can be considered of the form: The admissible and failure sets are then defined as follows:

∀j ∈ m 1 + 1, m , ±a j ≤ 0, Inria objective 2 (minimisation) objective 1 (minimisation) B 1 B 2 B 3 B 4
A = a ∈ R m | ∀j ∈ m 1 + 1, m , ±a j ≤ 0 , F = a ∈ R m | ∃j ∈ m 1 + 1, m , ±a j > 0 .
Contrarily to the classical Pareto dominance, the Constrained Pareto dominance compares boxes from the coupled objective/constraint spaces. The dominance rule when comparing two designs classically (without boxes) is defined as follows

a c b ⇐⇒ a 1,m1 b 1,m1 and a ∈ A or b ∈ F a c b ⇐⇒ a 1,m1 b 1,m1 and a ∈ A or b ∈ F a ∼ c b ⇐⇒ a c b and b c a
where the admissible set A contains the points satisfying all constraints while any point of the failure set F disrespects at least one constraint. Note that in the deterministic case, A = F , where the superscript • stands for complement. This relation would not hold in the context of boxed objectives and constraints, as shown in the following.

When comparing boxes of R m , the admissible A B and failure F B sets contain the boxes that lie entirely in A and F respectively:

A B = B ∈ ℘(R m ) | ∀a ∈ B, a ∈ A F B = B ∈ ℘(R m ) | ∀a ∈ B, a ∈ F RR n°9214

Pareto Optimal Probability

In the above, dominance rules allow to discriminate dominated boxes from non-dominated ones. This rule is exploited in [START_REF] Fusi | An adaptive strategy on the error of the objective functions for uncertainty-based derivative-free optimization[END_REF][START_REF] Rivier | Surrogate-Assisted Bounding-Box Approach for Optimization Problems with Approximated Objectives[END_REF] to refine only non-dominated boxes. To discriminate in a more rigorous way between non-dominated boxes, we propose to compute for each Bounding-Box its probability of being non-dominated. Such a computation is based on the assumption that the true robustness and reliability measures ρ in Figure 1 can be modelled with an aleatory variable following a uniform distribution within the Bounding Box B( ρ l , ε l ). The problem relapses to computing the Pareto Optimal Probability of aleatory variables with known uniform distributions.

With B all a given set of boxes B i i , the exact POP computation follows the formula below:

P OP true (B i ) = P {Z k } k ∩ Bj ∈B all j =i Z j c Z i , (1)
where for all k, Z k ∼ U B k . As assumed earlier, the realisations Z k are drawn uniformly within the set B k . Note however that this POP computation yields a combinatorial complexity that profoundly limits its calculation in closed form and would require Monte-Carlo approximation. Moreover, the returned score is sometimes not consistent in the presence of clusters of boxes in the objective and constraint space. For this reason, we propose two approximations of the POP, that will be quantitatively compared.

Ref. [START_REF] Teich | Pareto-Front Exploration with Uncertain Objectives[END_REF] proposes the use of a box probabilistic ranking as the fitness function for ESPEA (Estimated Strengh Pareto Evolutionary Algorithm). The score of a given box is computed by averaging the one-to-one domination probability with respect to all other boxes. Such an idea is used here, replacing the domination probability with the POP, as follows:

P OP av (B i ) = 1 N -1 Bj ∈B all j =i P Zj ,Zi Z j c Z i , (2)
with the same definition of Z k as above and where N is the total number of boxes in B all . The explicit details for computing the above expression are given in A. Since we compute each probability between only two boxes, the computational burden is very low. However, clustered boxes still yield inconsistent scores, as will be shown in the following.

A second metric, denoted as P OP min , is defined as follows, with B nd the set of non-dominated designs (using the Constrained Boxed Pareto dominance):

P OP min (B i ) = min Bj ∈B nd j =i P Zj ,Zi Z j c Z i . (3)
This approximation, yields a good relative ranking between boxes and is efficiently computable. It also shows robust behaviour when dealing with clustered boxes.

These three metrics i.e. P OP true , P OP av and P OP min are compared on three examples to illustrate the challenges associated with the POP computation. Each example is constituted of a given set of boxes and an optimisation problem: (1) a bi-minimization problem with highly clustered boxes (Figures 4), For its simplicity, interpretability and very low computational burden, P OP min will be used in the following.

SABBa framework

The objective of SABBa is to lower the computational cost of estimating the Pareto front by adaptively allocating the fidelity l(x) for computing the robustness and reliability measures ρ l Best and worst outcomes in green and red respectively. Associated POPs in Table 1 .

associated to a design x. The unconstrained version of this framework is developed in [START_REF] Rivier | Surrogate-Assisted Bounding-Box Approach for Optimization Problems with Approximated Objectives[END_REF] and relies on the coupling between the Bounding-Box approach presented above and a Surrogate-Assisting (SA) model built on the objectives in order to bypass function evaluations, notably toward the end of the optimisation. SABBa can be coupled to any optimisation algorithm and surrogate model for the SA strategy. This feature allows the approach to be very readily applicable and to benefit from new optimisation or metamodelling techniques. In this paper, the chosen optimisation algorithm is NOMAD [START_REF] Sébastien | Nomad: Nonlinear optimization with the mads algorithm[END_REF] for its reliable management of multiple objectives and constraints within a derivative-free framework (version 3.6.2). All surrogate models are handled with Gaussian Processes (using the python Gaussian Process package GPy [START_REF] Gpy | A gaussian process framework in python[END_REF] version 1.8.5).

In the following, we provide the global algorithm and the necessary numerical ingredients for the implementation. Specifically, Section 4.1 presents extensively the structure and algorithm of the framework, Section 4.2 deals with the computation and refinement of the Bounding-Boxes and Section 4.3 proposes a quality indicator for quantitative comparison on analytical test-cases.

Algorithm

The framework follows the structure depicted in Figure 6. The optimiser manages the design space exploration with the aim of computing an accurate approximation of the Pareto front. We illustrate now the typical sequence of operations within an iteration of the algorithm. Let us assume to have evaluated so far x i , ρ l (x i ) i for all the designs belonging to a set called X c (which includes each design evaluated by means of the Bounding-Box approach) and that the SA model ρ t SA and associate error ε t SA have been computed on these estimations. The algorithm features the following operation for each optimisation iteration, i.e. for a given design x:

• If the SA model yields accurate prediction on this design, meaning ε t SA (x) ≤ s 1 , the robustness and reliability measures returned by the surrogate ρ SA (x) are given to the optimiser;

• Else, the measures must be computed and the design x is added to the set X c . To this extent, f (x, ξ) and g(x, ξ) are computed for some values of ξ in order to build the surrogate models f (x, •) and g(x, •). This allows to estimate empirically, with a fidelity l mainly determined by the number of ξ samples, the robustness and reliability measures ρ l (x) = ρ l f (x), ρ l g (x) at very low cost. The Gaussian Process surrogate models also permit to compute an approximation ε l (x) of the associated error. The fidelity l at x can then be improved by increasing the size of the training data ξ i , f (x, ξ i ), g(x, ξ i ) on which the surrogates f (x, •) and g(x, •) are built.

In the case of designs estimated simultaneously, such as for the initial Design of Experiments (DoE) or for a generation in evolutionary algorithms, we propose to guide the refinements with the Pareto Optimal Probability (POP) of the estimated boxes, which is the probability of being not dominated. This refinement is performed until each box is either dominated, or has reached an estimated error ε l below a user-defined threshold s 2 . POP computation is presented in section 3.2.

Optimisation process

Bounding-Box approach

Surrogateassisting model Read new designs X t 7:

ρ t SA (x) ; ε t SA (x) Update X t ; Xc = x | ε t SA (x) > s1 ρ(x) = ρ l (x) ε(x) = ε l (x) ρ(x) = ρ t SA (x) ε(x) = ε t SA (x) ρ l (x) ; ε l (x) Dominated or ε l (x) ≤ s2 ∀x / ∈ Xc ∀x ∈ Xc Refinement: l(x * ) + + x * = argmaxP OP (x) t + +
for each x ∈ X t do 8:

if ε t SA (x) ≤ s 1 then 9:
Return ρ t SA (x) to the optimiser 10:

else 11: Add x in X c 12:
Compute a first approximation ρ 0 (x) of ρ 13:

end if

14:
end for 15:

while X r = x ∈ X P B ∩ X c | ε l (x) > s 2 is non-empty do 16: Find x * = argmax x∈Xr P OP (x) 17: l(x * ) + + 18:
Compute f (x * , ξ) and g(x * , ξ) at some ξ 19:

Update ρ l (x * ) and ε l (x * ) 20:

end while

21:

Update ρ t SA and ε t SA with the new ρ l (x) and ε l (x) values

22:

Return ρ l (x) to the optimiser 23: end while X P B contains all designs which associated boxes are non-dominated and X c refers to the designs which measures have not been returned by the SA model. Hence, the boxes in X r are all estimated without the SA model and both non-dominated and not refined up to s 2 .

Remark We consider to loop over finer and finer values of s 1 and s 2 . This allows more accurate intermediate results and slightly cheaper overall convergence. Practically, vectors of thresholds can be specified by the user.

The algorithm relies on a quite small set of parameters, namely:

• The predefined sequence of pair of thresholds (s 1 , s 2 );

• The number of function evaluations N f irst for the first approximations ρ 0 and N new the number of additional evaluations for refinement;

• The number of designs N at each optimisation iteration and N init at the first iteration.

The optimiser may impose these parameters.

Bounding-Boxes computation and refinement

To compute the boxes, robustness and reliability measures must be estimated, and the associated error must be quantified. In particular, we focus in this paper on the following measures:

Expectation, variance, minimum, maximum and quantile. The values returned to the optimiser are denoted by ρ(x) and ε(x). As shown in Fig. 6 and Alg. 1, they can be computed with tunable fidelity using the Bounding-Box approach or can be computed on the SA model. These two computations are presented in details in the following sections. Note that we make use of Gaussian Processes as surrogate models, giving both a predictive value and an associated predictive variance. Box widths are then computed based on this variance, relaxing the conservative assumption to a ±3σ paradigm. Practically, with σ 2 (z) the predictive variance of a GP surrogate model, an error can be computed as:

ε(z) = 3σ(z).
Although this does not imply ε(z) ≥ |ε(z)|, with ε the true error of the surrogate, the probability of dissatisfying the conservative assumption is P ε(z) < |ε(z)| < 0.3%.

Computed boxes

For each design x where the SA model is not exploited, the boxes must be directly computed from samples of f and g in the uncertain space. This corresponds to lines 9 to 15 in Algorithm 1. Gaussian Processes are used to approximate f (x, •) and g(x, •). Two approaches are considered:

• Separated Spaces: A GP is built in Ξ at a given x using a set of samples f (x, ξ i ), g(x, ξ i ) i .
This directly corresponds to f (x, •) and g(x, •).

• Coupled Space: A GP is built in X × Ξ using a set of samples f (x i , ξ i ), g(x i , ξ i ) i .
This gives approximations f and g of f and g. At a given x, the cuts f (x, •) and g(x, •) can be returned.

Remark While the Separated Spaces (SS) approach builds a GP in a space of lower dimension (Card(Ξ)) than the Coupled Space (CS) approach (Card(X × Ξ)), thus easing the surrogate modelling task, it does not use surrounding samples in the design space and restarts from scratch at each new design point x. These strategies will be compared on analytical test-cases in Section 5.

Both SS and CS approaches allow having predictive models f (x, •) and g(x, •), that will be written f x and g x , and their predictive variance σ f (x, •) and σ g (x, •). The fidelity l of these models, used to compute ρ l (x) and ε l (x), refers to the number of training data. This level of fidelity is implicit in the following equations. The ±3σ paradigm provides the box widths, denoted as ε fx = 3σ f (x, •) and ε gx = 3σ g (x, •). For each output measure, indexed by k, the classical empirical estimators are exploited for the computation of ρ k (x) = ρ l k (x) and the associated width ε k (x) is given in Eq. ( 4) for f x (the same can be written for g x ):

ε k (x) = ε l k (x) =                E ξ ε fx (ξ)
for expectation

E ξ ε µ (x) + ε fx (ξ) 2 + 2 f x (ξ) -µ(x) ε µ (x) + ε fx (ξ) for variance max m(x) -min ξ [ f - x (ξ)] , m(x) -min ξ [ f + x (ξ)] for minimum max M (x) -max ξ [ f - x (ξ)] , M (x) -max ξ [ f + x (ξ)] for maximum max q p (x) -q p ξ [ f - x (ξ)] , q p (x) -q p ξ [ f + x (ξ)] for quantile (4) where f + x (ξ) = f x (ξ) + ε fx (ξ) and f - x (ξ) = f x (ξ) -ε fx (ξ).
In practice, the expected values above are approximated by means of Monte Carlo Sampling (MCS) on surrogate model at very Inria low cost. A formal justification of these box sizes is given in B.

We also provide guidance for refining the surrogate model f x (ξ). While one could rely on a space-filling paradigm, geometrically filling the widest hole in the sampling, we choose here to use GP-based refinement criteria. A partial criterion is first computed for each statistical measure, with the following formulas:

c k (ξ) =                            ε fx (ξ)p(ξ)
for expectation

ε fx (ξ)p(ξ) for variance m -f - x (ξ) 2ε fx (ξ) + for minimum f + x (ξ) -M 2ε fx (ξ) + for maximum q p -f - x (ξ) 2ε fx (ξ) + f +
x (ξ) -q p 2ε fx (ξ) + p(ξ) for quantile [START_REF] Papadrakakis | Reliability-based structural optimization using neural networks and monte carlo simulation[END_REF] with the same definition of f + x (ξ) and f - x (ξ) as before and where • + = max(0, •). Note that we heuristically choose here to multiply the criteria for the mean, variance and quantile measures by the input Probability Density Function (PDF) p(ξ) in order to put more weight on the most likely area. Justifications for these formulas are provided in C. These partial criteria are then scalarised into the final refinement criterion through the following weighted sum:

c(ξ) = m k=1 w k ck (ξ) (6)
with ck (ξ) the normalised partial criteria, m the number of measures and w k the weights. The normalised partial criteria are computed as follows:

ck (ξ) = c k (ξ) -min ξ [c k (ξ)] max ξ [c k (ξ)] -min ξ [c k (ξ)] ∈ [0, 1].
We propose here to compute the weights as the ratio between the conservative error ε k and the target accuracy. In this manner, any partial criterion associated with high error compared to the target accuracy will profoundly influence the final criterion. Practically, to emphasise this behaviour, the dependence to the conservative error is chosen quadratic:

w k = ε k s 2 k 2 .
Note that in the case of multi-point refinement, we conduct a greedy sequential approach by assuming that previous refinements are performed, fixing the predictive values of the GP surrogate model and recomputing the predictive variance to obtain the updated refinement criterion. We can then perform Black-box evaluations in parallel on these points. This approach allows performing multi-point refinement efficiently without the need of any clustering heuristic. However, this strategy makes a lot of assumptions and may not yield optimal results when performing many refinements in parallel.

SA-based boxes

When the SA model error is low enough, below s 1 in all output dimension, the predictive value of the SA model is directly returned to the optimiser. As stated previously, the SA model is a Gaussian Process, which gives access to the predictive variance in order to compute box widths.

The SA model is built on the previously computed boxes to return an approximation ρ SA of ρ in X . Because this model is built on the previously computed boxes and evolves throughout the optimisation steps, the SA model is referred as ρ t SA , with t the optimisation step. In practice, the available data for constructing this model are boxes of different sizes. These boxes can be considered as noisy evaluations of ρ at different x, with heterogeneous noises. Using again the ±3σ paradigm, these different noises can be translated into heterogeneous variances

σ 2 i = ε l (x i ) 3
2 with x i the designs at which boxes were computed using the direct approach described in 4.2.1. This allows to make use of heteroscedastic Gaussian Processes, which naturally take into account heterogeneous noise variances, under the gaussianity assumption. The measure estimations ρ and associated error ε at a new design x are computed as follows:

ρ(x) = ρ SA (x) = k T * (K + ∆) -1 ρ l ε(x) = 3σ SA (x) = 3 k * * -k T * (K + ∆) -1 k * . (7)
The superscript t is not explicitly written as the above predictor is updated at each optimisation step. Here, K = K {x i } i , {x i } i represents the autocovariance matrix between training points, k * * = k x, x at the new design and k * = k x, {x i } i the covariance vector between x and training points. The diagonal matrix ∆ = diag σ 2 i i represents the heterogeneous gaussian noises presented previously. More details on GP surrogate models can be found in [START_REF] Edward | Gaussian processes in machine learning[END_REF].

Quality indicator

Assessing and comparing the performances of the method requires the computation of a quantitative quality indicator. To this extent, the Hausdorff distance is a classical choice for computing the closeness of the found optimal set to the true one. In practice, the modified Hausdorff distance proposed in [START_REF] Dubuisson | A modified hausdorff distance for object matching[END_REF] captures the similarities more efficiently. With Bounding-Boxes, there is no Pareto Optimal set but rather a set of non-dominated designs, each of them having a Pareto Optimal Probability. The approximated Pareto front P and Pareto optima X P are hence aleatory. We propose to compute the expected value of the modified Hausdorff distance with respect to the realisations of X P .

The deterministic modified Hausdorff distance d H is given as follows:

d H (A, B) = max d 1 (A, B), d 1 (B, A) , d 1 = 1 N a∈A d 2 (a, B), d 2 (a, B) = min b∈B [ a -b 2 ].
Under the box independence and uniformity assumptions, as in Section 3.2 we have ∀x,

ρ(x) ∼ U B ρ(x), ε(x) ,
where ρ(x) and ε(x) are the values returned to the optimiser in Figure 6. The proposed expected modified Hausdorff distance reads:

Q B = E X P d H (X P , X P ) . (8)
Inria where a realisation of X P corresponds to a realisation of ρ(x) for all x and the computation of the classical Pareto optima.

In the next section, the SABBa framework for optimisation under uncertainty problems is applied to analytical test-cases to assess its performance against classical approaches. The quality of the optimisation outputs is quantified using the quantitative indicator presented above.

Analytical comparisons

We apply here SABBa to two analytical test-cases. The first analytical test-case is low-dimensional and deals with a Taguchi multi-objective robustness formulation (mean optimisation and variance minimisation). The second test-case features a higher number of dimensions and consists in a mean optimisation under quantile constraint.

We systematically compare the performance of the proposed framework to an A Priori Meta-Model (APMM) strategy. In the latest, we build a surrogate model of the output in the Coupled Space (CS) before the optimisation process. A Double Loop, performing an Uncertainty Propagation at each optimisation iteration, can then be conducted on the surrogate at a very low cost. The overall computational cost of building this surrogate model is here chosen a priori and can be adapted to the available budget, with a substantial impact on the accuracy of the results.

We show in the first test-case the performance of the direct Double Loop approach with respect to the APMM strategy. The former performs N U P evaluations for measure computation at N opt designs, explored by the optimiser. It is shown to be extremely inefficient.

We compare here three variants of SABBa, denoted as follows:

• SA-SS: when using the Surrogate-Assisting strategy and Separated-Spaces models for measures computation;

• SA-CS: when using the Surrogate-Assisting strategy and a Coupled-Space model for measures computation.

• CS: when using the Coupled-Space model without the Surrogate-Assisting strategy;

Among the above variants, one can expect SA-CS to give better results than SA-SS and CS, taking advantage of both the low-dimensional Surrogate-Assisting model and the coupled space correlation to speed up box refinement.

Ten runs are performed for each strategy to capture both the mean convergence curve and the associated variability, represented as a translucent band around the mean. The log distance is assumed to show a Gaussian distribution over the repeated runs, implying a log-normal distribution over the actual distance. Note that a high variability of the convergence curves reveals a lack of reliability of the associated approaches.

Test-case 1: Unconstrained Taguchi optimisation

This problem is a bi-objective robust optimisation proposed in [START_REF] Fusi | An adaptive strategy on the error of the objective functions for uncertainty-based derivative-free optimization[END_REF]. There are two design variables x and one uncertain parameter ξ, and the problem reads:

minimise: ρ f (x) = µ(x) σ 2 (x)
where:

µ(x) = E ξ [f (x, ξ)] σ 2 (x) = V ξ [f (x, ξ)] with: f (x, ξ) = ξ -x 1 ξ 5 + cos(2πx 2 ξ) + 5 ξ ∼ U([0, 1]) by changing: (x 1 , x 2 ) ∈ [1, 2] 2 (9)
The Pareto front associated with this problem is discontinuous, and the optimal set in the design space is the union of a segment and a point (see Fig. 7). The Double Loop and A Priori MetaModel strategies are quantitatively compared in Figure 8. Each Uncertainty Propagation within the Double Loop is performed up to a relative accuracy of 5%. A very slow convergence of the Double Loop approach is observed.

For this first test-case, we consider two different surrogate-modelling capabilities. The aim is to illustrate SABBa performance both when the surrogate is able to represent the underlying functions accurately and when it shows poor efficiency. This highly impacts the quality of the APMM f , the Surrogate-Assisting model ρ SA (x) and the separated or Coupled-Space models f (x, •).

With high-quality surrogate model

The unconstrained Taguchi optimisation problem is first solved using high-quality metamodelling approaches. Practically, we build a heterogeneous GP surrogate model with a classical RBF kernel (also called squared-exponential, exponential quadratic or Gaussian), and Auto-

k(x, x ) = σ 2 exp - r 2 l 2 with r l = m i=1 x i -x i l i 2
We recall that the kernel function gives the covariance matrix in Equation ( 7). This model requires m + 1 hyperparameters {σ 2 , l 1 , . . . , l m } to be optimised but captures the characteristic lengthscale associated with each input dimension.

The convergence curve of the indicator defined in Equation ( 8) is plotted in Figure 9. Figures 9(a) and 9(b) show a very reduced variability of SABBa compared to the APMM strategy. However, only SA-CS (purple curve in Fig. 9(b)) shows a significant improvement for the mean convergence.

Note that here, the thresholds s 1 and s 2 are sequentially refined five times. To alleviate the tuning of these thresholds, SABBa can deal with normalised thresholds s1 and s2 . At each iteration, the range h i covered in the i th -dimension is updated, namely, ∀i,

h i = max x [ρ i (x)] -min x [ρ i (x)],
and the thresholds s1 and s2 are given in percentage of h. Here, s1 and s2 are both sequentially taken as 50%, 40%, 30%, 20%, 10% and 5% in all dimensions. The chosen parameters of SABBa are as follows: N init = 10, N = 1 (sequential optimiser), N f irst = 5 and N new = 1 (sequential refinement).

Figure 9 pictures the convergence curves of the APMM strategy and the three studied SABBa variants (SA-CS, SA-SS and CS). The mean curves are mostly comparable, with a slight advantage for SA-CS. However, the associated variability is much smaller when using any SABBa variant with respect to the APMM strategy. The indicator is assumed to follow a normal distribution in log scale, hence to be log-normal. The high variance associated to the APMM strategy coupled with the heavy tail of the log-normal distribution makes the approach very unreliable. We picture this risk in Figure 10, where the worst result out of ten repeated runs of SA-CS and the APMM strategy are plotted. These outputs correspond to N eval ≈ 100.

In these figures, optimal designs returned by SABBa are plotted in greyscale. This refers to the Pareto Optimal Probability (POP) of each design, ranging continuously from 0 (white) to 1 (black). Inria

With a low-quality surrogate model

In the following, the same optimisation problem is solved using a low-quality surrogate model in both SABBa and the APMM strategy. Contrarily to the previous formulation, we build a homogeneous GP surrogate model with only one lengthscale, that has to account for all dimensions. Practically,

k(x, x ) = σ 2 exp - r 2 l 2 with r l = 1 l m i=1 x i -x i 2
Only two hyperparameters {σ 2 , l} must be optimized. However, this model will notably fail when the characteristic lengths of the function in the different dimensions are very disparate. This test-case aims at simulating problems where the coupled space behaviour is hard to model. This would naturally arise when the number of dimensions is significant. Hence, one can expect the SA strategy to yield a significant cost improvement through low-dimensional measures surrogate modelling.

Such improvement is indeed revealed in Figure 11. The use of the SA strategy has two consequences: (i) non-Surrogate-Assisted approaches (APMM and SABBa CS) show poor performance compared to the other strategies and (ii) the gap between SA-CS and SA-SS is much narrower. Indeed, here, they both rely for the most part on the low-dimensional SA model. As previously, the worst optimisation outputs are plotted for SABBa SA-CS and the APMM approach. In Figure 12, SABBa shows much higher consistency and accuracy compared to APMM.

The use of coupled space surrogate models and Surrogate-Assisting strategy have shown to bring an significant cost reduction for the SABBa SA-CS framework. It performs better in average than the APMM strategy and shows far greater consistency and robustness.

Test-case 2: Quantile-constrained mean performance optimisation

We propose this second test-case to assess the performance of SABBa in a higher dimensional case and in the presence of a reliability-based constraint (derived from the Six-Hump Camel function). The objective is a robustness measure derived from a simplified Rosenbrock function.

We consider here four design variables and three uncertain parameters. The problem is stated as follows:

minimise:

ρ f (x) = µ(x) satisfying: ρ g (x) = q 95% (x) ≤ 1 where: µ(x) = E ξ [f (x, ξ)] q 95% (x) = q 95% ξ [g(x, ξ)] with: f (x, ξ) = 3 i=1 1 -x i + 3 1 + arctan 5(ξ i -0.5) 2 x i+1 -x 2 i 2 g(x, ξ) = 4 -2.1x 2 1 + x 4 1 3 x 2 1 + x 1 x 2 + -4 + 4x 2 2 x 2 2 + cos(2πξ 1 ) -sin( π 2 ξ 1 ) -ξ 1 -cos(2π0.05) -sin( π 2 0.05) -0.05 5 ξ ∼ U([0, 1] 3 ) by changing: x ∈ [-0.2, 1.2] 4 (10)
One can note that the robustness measure reduces to the classical formulation of the 4D

Inria

Rosenbrock function and the reliability measures to the Six-Hump Camel function. Analytically, it holds that:

µ(x) = 3 i=1 1 -x i + 3 x i+1 -x 2 i 2 , q 95% (x) = 4 -2.1x 2 1 + x 4 1 3 x 2 1 + x 1 x 2 + -4 + 4x 2 2 x 2 2 . ( 11
)
The optimum of this deterministic problem is found at x * ≈ (0.7033, 0.7035, 0.6212, 0.3859) with q 95% (x * ) = 1 and µ(x * ) ≈ 0.4981.

For this test case, both s1 and s2 are sequentially taken as 50%, 40%, 30%, 20%, 10%, 5%, 3%, 2%, 1% and 0.5%. As for the other parameters, N init = 10, N = 1 (sequential optimiser), N f irst = 5 and N new = 1 (sequential refinement).

Figure 13 shows a slight mean improvement when using SA-CS compared to the APMM strategy. The variability of the output is also halved with all SABBa variants, which is critical for real-world applications. Note that the plateau that is reached by the SA-CS curve is actually a plotting artefact. All runs have different final number of evaluations, and we considered the indicator Q B constant when the optimum is reached. Therefore, the average of these curves tends to flatten at the end. As done before, we plot the worst optimisation outputs in Figures 14. They are depicted in parallel coordinates plots, with a P OP min greyscale. Each efficient individual is represented by a grey curve and the true optimum is in red. This comparison reveals again the high unreliability of the APMM strategy with respect to the SABBa framework.

6 Physical applications

Two-bar truss

The two-bar truss optimization problem is notably illustrated in Refs. [START_REF] Baudoui | Optimisation robuste multiobjectifs par modèles de substitution[END_REF] and [START_REF] Jin | The use of metamodeling techniques for optimization under uncertainty[END_REF]. The uncertainty-based optimization problem is formulated as follows (with a schematic representation in Figure 15(a)): minimise:

ρ f (x) = V (x) satisfying: ρ g (x) = q 1 (x) q 2 (x) ≤ s max 0
where: q 1 (x) = q 0.999 [s(x, ξ)]

q 2 (x) = q 0.999 [s(x, ξ) -s crit (x, ξ)] with: ξ ∼ N 150000 210000 , 30000 2 0 0 21000 2 by changing: x ∈ [20, 80] × [800, 1200] × [700, 800] × [2, 3] (12)
In the above, the design variables are: (x 1 ) the diameter of the cross section d, (x 2 ) the bar length L, (x 3 ) the structure half-width B and (x 4 ) the thickness of the cross section T . The uncertain parameters refer to (ξ 1 ) the external force F and (ξ 2 ) the elastic modulus E. The objective is to minimize the total volume V (x) = 2πx 1 x 2 x 4 × 10 -6 while verifying that the probability of s(x, ξ) =

x 2 ξ 1 2πx 1 x 4 x 2 2 -x 2 3 exceeding s max = 400 N.mm -2 and s crit (x, ξ) = π 2 ξ 2 (x 2 1 + x 2 4 ) 8x 2 2 are both below 0.001.
The problem is represented in Fig. 15(b), where the objective value (V (x), to be minimized) Inria is plotted in the constraints space with a color scale. The abscissa and ordinate directions refer to the constraints q 1 (x) and q 2 (x) ; the associated thresholds are drawn with red lines. The admissible set is in the lower left quadrangle (q 1 < 400 N.mm -2 and q 2 < 0 N.mm -2 ) and one may note that the objective is nearly constant on the limit q 1 = 400 N.mm -2 . This optimisation has been solved with the SA-CS variant. The thresholds s 1 and s 2 are sequentially reduced up to 0.5% of the total ranges. A set of 6 non-dominated designs is returned, with different Pareto Optimal Probabilities (POP). This is reported in Figure 16, where the optimal designs are drawn in the output space with POP greyscale (black: 1; white: 0). Figure 16(a) shows the deterministic objective values associated to each non-dominated design ; the ordinate only avoids overlapping points for better readability. Figure 16(b) shows the boxes in the constraints space q 1 and q 2 . One can note that the boxes associated to the highest POP are located mostly in the admissible space (q 1 ≤ 400 N.mm -2 ). This optimisation was performed with 65 function evaluations, which is very inexpensive with respect to the accuracy of the output.

The design associated to the highest POP is x * = (63.292, 1021.4, 701.25, 2.0812). We can verify that P ξ [s(x * , ξ) ≥ s max ] 1.2 × 10 -3 and P ξ [s(x * , ξ) ≥ s crit (x * , ξ)] ≤ 1.0 × 10 -4 . The first quantile constraint is not perfectly satisfied. This is due to (i) the associated box not being entirely in the admissible space, (ii) the final accuracy not being very sharp and (iii) the quantile computation being performed with simple Monte Carlo, hence with high variance. A dedicated technique for accurate quantile computation would be required to avoid adding estimation noise, that is not taken into account here.

While the objective and q 1 values are well captured (optimal and close to constraint limit), the q 2 seems not to have a significant impact. To verify this, the optimisation is solved two more times. The results are plotted in Figure 17.

As expected, while similar optimal values for V and q 1 are found, the associated values of q 2 RR n°9214 vary dramatically. Specifically, the optimal V (in L) is in the ±10 -2 L range around the true optimum while the total spanned range is of 1.6 L. Similarly, optimal q 1 varies of less than 10 N.mm -2 over the total range of more than 2500 N.mm -2 . However, q 2 ranges approximately from -850 N.mm -2 to -150 N.mm -2 . This is far from negligible, as it can be seen in Figure 18 and validates the observation that V seemed nearly constant on the q 1 constraint limit. Note that in all three runs, designs associated to an objective value below the optimum have boxes that lay partially in the failure region q 1 > 400 N.mm -2 . Such an objective value cannot be reached while satisfying the constraints, which reveals that the true q 1 value are indeed over 400 N.mm -2 .

These two additional runs were performed with 72 and 65 function evaluations. Finer thresholds and additional evaluations would be required in order to find a unique global optimum on the q 1 = 400 N.mm -2 line. However, up to the chosen 0.5% threshold, the optimum has been well captured at a very low cost.

Design of a Thermal Protection System

The SABBa framework is then applied to the design of a thermal protection system for a re-entry vehicle. This test-case deals with 12 dimensions and aims at minimising the mean mass density under worst-case temperature constraint.

We study the re-entry of Stardust, that was the first mission using a low-density carbonphenolic ablator in 2006. Stardust was the fastest man-made object re-entering the earth atmosphere, at a velocity of 12.7 km/s. Surface total pressure and heat flux were computed from hypersonic computational fluid dynamics (CFD) simulations. In accordance with the state-ofthe-art design approach, we assume that the problem is locally mono-dimensional. The analysis is performed using the properties of the Theoretical Ablative Composite for Open Testing (TACOT). Its composition and properties are comparable to PICA ; nominal TACOT properties are available in the open literature. Volume-wise, TACOT is made of 10% of carbon fibers, 10% of phenolic resin, and is 80% porous. The thickness of the ablative material is two inches and adiabatic conditions are used at the bondline.

The physical model used here is a generic heat and mass transfer model for porous media that is presented in [START_REF] Lachaud | Porous material analysis toolbox based on openfoam and applications[END_REF]. The model is implemented in the Porous material Analysis Toolbox (PATO), distributed Open Source. First-order implicit finite-volume schemes in time and space were used for the simulations. The mono-dimensional problem was meshed with 300 finite-volume cells with a logarithmic refinement of parameter 0.2 towards the surface. In this study, we used an equilibrium chemistry model. We study the material response at the stagnation point, which reaches the highest temperature during the reentry. Figure 19 illustrates the typical evolution of the temperature inside the material and at the heated surface.

with x 1 the resin volume fraction (originally 10%) an x 2 the overall width of the system (originally 7.21 cm). We observe in Figure 19 that with initial parameters, the bottom temperature stays merely constant during the whole reentry. Hence, the search space is centred on lower x 1 and x 2 values.

In this test-case, Ξ is of dimension 12 and all uncertainties are assumed uniform. The choice of the uncertain parameters relies on a sensitivity analysis performed on this test-case in a previous paper [START_REF] Rivier | Ablative thermal protection system under uncertainties including pyrolysis gas composition[END_REF]. Here, retained uncertain parameters are: 1) Density and volume fraction of the fibrous preform (5% uncertainty each); 2) Density and volume fraction of the phenolic resin (5% uncertainty each); 3) Thermal properties of the charred material: heat capacity, conductivity and emissivity (5% uncertainty each); 4) Oxygen fraction in the pyrolysis gases (10% uncertainty); 5) Pyrolysis reaction activation energy (10% uncertainty); 6) Overall width of the system (0.1 cm uncertainty). Both design parameters are also affected by an uncertainty (±0.005 for x 1 and ±0.1 for x 2 ).

Inria

Each function evaluation requires the 1D simulation of heat transfer and shield ablation. This takes approximately 10 minutes to compute on a 2.90 GHz processor. An acceptable global cost should remain within a day, or equivalently below 100 to 150 evaluations. , where ρ f must be minimised and the admissible set is below the constraint threshold ρ g = 473.15 K. These optima were found using again the SABBa framework with Coupled-space model (SA-CS) for a computational cost of only 40 function evaluations. In Figure 21, we have reported the performances of all the designs evaluated during the optimization. They are plotted in blue when dominated, in red when entirely in failure zone and in green when non-dominated.

Among the three optimal designs depicted in Figure 20, only the right one lies entirely in the admissible set. Using the notations from Section 3.1.2, this box belongs to A B . Hence, contrarily to the two other optimal designs, it can dominate other boxes. In Figure 21, it can be observed that this design dominates all the blue boxes, that are not entirely in the failure zone but are strictly worse in the objective dimension. Contrarily, the red boxes belong to F B and lie entirely in the failure zone. Hence, they are all considered dominated. This behavior explains why only three designs are kept, as shown in Figure 20. One can note that the tunable fidelity has been well exploited by the Bounding-Box approach in Figure 21, where extreme boxes (on the very left or very right) are estimated with very wide uncertainty.

Figure 22 provides a visualisation of the final SA model. In particular, we plot the SA mean surface density ρ SA f (x), wherever the constraint ρ SAg (x) ≤ 473.15 k is satisfied. One can see that (i) the optima plotted in Figure 20 are coherent with this SA model, and (ii) the objective value ρ SA f (x) seems merely constant on the constraint limit. The final accuracy is imposed at 1%, which may not be enough to discriminate designs on the constraint limit. The above found x * ≈ (4.43 × 10 -2 , 4.25) is returned by SABBa but the whole frontier from (0.01, 5.1) to (0.07, 4) is of high interest.

Conclusions

In this work, we present an efficient framework, SABBa, for constrained optimisation under uncertainty problems. It handles most classical robustness and reliability measures, such as the Taguchi (mean and variance) robust optimisation or the quantile constrained optimisation problems. The parsimony of SABBa relies on several features presented in this paper. The Bounding-Box measures approximations is used to compare designs in a tunable fidelity context. It permits to compute the Pareto Optimal Probability (POP) associated to each design to perform a rigorous ranking among them. Gaussian Process (GP) surrogate models are exploited both for measure estimation and Surrogate-Assisting strategy, and allows to tract the estimation variability throughout the framework.

Inria

We have assessed the performance of the framework on both analytical and applicative test cases. We have systematically compared SABBa to the A Priori MetaModel strategy. To this extent, we have formulated a specific indicator based on a modified Hausdorff distance. Overall, SABBa shows both a faster convergence rate and a dramatic increase in the robustness of the optimisation process.

It has been successfully applied to a low-quantile reliability-based structural optimisation with six dimensions in the coupled space. The accuracy is very satisfactory, with an associated computational cost of less than 70 evaluations. A 12-dimensional Thermal Protection System (TPS) design with worst-case constraints has also been performed with only 40 evaluations.

This approach remains very general and broadly applicable to any optimisation process with statistics-based objectives and constraints. It shows very good parsimony, which is of primary importance for real world applications, where black-box evaluations can take several days to run. Most importantly, the optimiser is not imposed by the approach, allowing for easy coupling to any existing deterministic process.

Several steps in this work remain improvable. We used local refinement criteria, which are cheaper but less efficient than integral criteria. The use of GP surrogate models also limits the number of manageable input dimensions. Surrogate modelling on over few dozens of dimensions would require specific techniques such as feature selection or kernel adaptation. Likewise, in order to efficiently deal with very low probability constraints (e.g. 1×10 -7 ), one should couple SABBa to a dedicated tool for low quantile computation. Following actions could be the application of SABBa to other costly engineering cases, the inclusion of Value at Risk and Conditional Value at Risk (VaR and CVaR) measures, the definition of non-uniform boxes and the development of a specific Bayesian Optimisation technique for SABBa.

A POP computational details

The set of boxes B = B i i is given, with B i = B a i , r i . Every Z i ∼ U B i is of dimension m = m 1 + m 2 where m 1 is the number of objectives and m 2 the number of constraints. Hence, Z i 1,m 1 = Z i f are the objectives values and Z i m 1 +1,m = Z ig the constraints values. The probability in Equation ( 2) is computed as follows, using the independence assumptions between boxes and between dimensions:

P Zj ,Zi Z j c Z i = 1 -P Zj ,Zi Z j c Z i = 1 -P Zi g Z ig 0 (14)
+ P Zi g Z ig ≤ 0 P Zj g Z jg ≤ 0 P Zj f ,Zi f Z j f Z i f = P Zi g Z ig ≤ 0 1 -P Zj g Z jg ≤ 0 P Zj f ,Zi f Z j f Z i f where P Zi g Z ig ≤ 0 = m k=m1+1 max 0, min 1, -B - i k 2r i k = m k=m1+1 -B - i k 2r i k 1 0
, and • 1 0 means that values are taken between 0 and 1. The lower bound of the box is written as B - i = a i -r i . The second probability involved in the last line of Equation ( 14) is computed RR n°9214

as follows:

P Zj f ,Zi f Z j f Z i f = m1 k=1 P Zj k ,Zi k Z j k ≤ Z i k = m1 k=1 L 1 k 2r j k 1 0 + L 2 k 2r j k 1 0 1 2 L 2 k 2r j k 1 0 + L 3 k 2r j k 1 0 ,
where

L 1 k = B - i k -B - j k , L 2 k = min(B + i k , B + j k ) -max(B - i k , B - j k ), L 3 k = B + i k -B + j k .
In the above, B + i = a i + r i is the upper bound of the box. In each dimension, L 1 k can be interpreted as the portion of B j dominating B i , L 3 k the portion of B i dominated by B j and L 2 k the overlapping area. We depict these lengths in Fig. 23. Note that these values can be negative and that the computation is not symmetric between B j and B i .

B j k B i k L 1 k L 2 k L 3 k

B Justification for box sizes

Some developments are given here to explain the choice of the conservative error in Equation ( 4) associated with measure approximations. The expectation approximation conservative error ε µ is quite straightforward:

ε µ (x) = µ(x) -µ(x) = E ξ [f (x, ξ)] -E ξ [ f x (ξ)] = E ξ [ f x (ξ) + ε fx (ξ)] -E ξ [ f x (ξ)] = E ξ [ε fx (ξ)] ≤ E ξ [ε fx (ξ)] = ε µ (x)
with µ(x) the true statistical moment on f at a given x and µ(x) the approximated one, computed on the surrogate model f x . The true error ε fx is defined from f (x, ξ) = f x (ξ)+ε fx (ξ) and is conservatively approximated by ε fx (ξ) ≥ |ε fx (ξ)|.

Inria

The same can be conducted for the variance measure, with a bit more risks of overestimation:

ε σ 2 (x) = σ 2 (x) -σ 2 (x) = E ξ [(f (x, ξ) -µ(x)) 2 ] -E ξ [( f x (ξ) -µ(x)) 2 ] = E ξ [( f x (ξ) -µ(x) + ε fx (ξ) -ε µ (x)) 2 ] -E ξ [( f x (ξ) -µ(x)) 2 ] = E ξ [(ε fx (ξ) -ε µ (x)) 2 + 2( f x (ξ) -µ(x))(ε fx (ξ) -ε µ (x))] ≤ E ξ [(ε fx (ξ) + ε µ (x)) 2 + 2| f x (ξ) -µ(x)|(ε fx (ξ) + ε µ (x))] = ε σ 2 (x)
The conservative error associated with the expectation approximation could also be explained based on the monotonicity of the expectation operator. We write f ≥ g ⇐⇒ ∀x, f (x) ≥ g(x). By defining f + x (ξ) = f x (ξ) + ε f x (ξ) and f - x (ξ) = f x (ξ) -ε f x (ξ), the monotonicity of the minimum, maximum and quantile operators can also be exploited to obtain conservative error approximations. For the case of the minimum operator, with m = min . The same idea can be followed for the maximum and quantile approximations, resulting in the following errors:

ε max (x) = max M (x) -max ξ [ f - x (ξ)] , M (x) -max ξ [ f + x (ξ)]
ε q p (x) = max q p (x) -q p ξ [ f - x (ξ)] , q p (x) -q p ξ [ f + x (ξ)]

C Justification for refinement criteria

The chosen partial criteria for GP model refinement given in Equation ( 5) are also explained in the following. Both the expectation and variance are global measures over the whole domain. For this reason, it has been chosen to iteratively add a point at the maximum predictive conservative error, to converge the model on the entire space. This strategy is usually called Maximum Mean Square Predictive Error or MMSPE. Note that the partial criteria are multiplied by the input pdf to weight the predictive error according to the probability of occurrence. Thus, the final partial criteria are:

c µ (ξ) = ε fx (ξ)p(ξ) c σ 2 (ξ) = ε fx (ξ)p(ξ)
The minimum (resp maximum) measure partial criteria is simply the probability of exceeding the current minimal (resp. maximal) value m (resp. M ). This is performed with an assumption of uniform distribution with the conservative error box. Hence, the criteria can be written as such:

c min (ξ) = m -f - x (ξ) 2ε fx (ξ) + c max (ξ) = f + x (ξ) -M 2ε fx (ξ)
+ where f + x (ξ) = f x (ξ) + ε f x (ξ) and f - x (ξ) = f x (ξ) -ε f x (ξ) and with • + = max(0, •) referring to the value if positive, 0 either.

Finally, for the case of the quantile measure, it has been chosen to compute the product of the aforementioned probability of exceeding the quantile value, multiplied by the input density:

c q p (ξ) = q p -f - x (ξ) 2ε fx (ξ) + f + x (ξ) -q p 2ε fx (ξ) + p(ξ)
This product is maximised on the hyperplane f x = q p . Multiplying the criteria by p(ξ) puts more weight according to the probability of occurrence. The spread may not be optimal, but the optimisation of the criteria being performed on a fixed sampling, the tightening of the conservative error is very likely to decrease the value of the criterion in the surrounding area because samples are very unlikely to be localised exactly on the isoline of the function.

We can note that the chosen criteria are far from being optimal. However, these choices give a fast determination of the refinement point, and the number of samples is usually low enough so that the choice of the training point is not of significant importance compared to the quality of the metamodelling strategy.

Figure 1 :

 1 Figure 1: Bounding-Box approximation

Figure 2 :

 2 Figure 2: Boxed Pareto Dominance: Comparison of 4 boxes, best and worst outcomes in green and red, respectively.

Inria

4 Figure 3 :

 43 Figure 3: Constrained Boxed Pareto Dominance: Comparison of 4 boxes, best and worst outcomes in green and red, respectively.

  (2) a bi-minimization problem (Figure 5 (a)) and (3) a constrained mono-objective minimization problem (Figure 5 (b)

Figure 5 :

 5 Figure 5: (a) Example 2 with a bi-minimization problem ; (b) Example 3 with a constrained mono-objective minimization. Best and worst outcomes in green and red respectively. Associated POPs in Table1.

Figure 6 :Algorithm 1

 61 Figure 6: Structure of the SABBa framework

7 :

 7 Test-case 1: a) Discretisation of the design space in red and Pareto optimal sets in black. b) Image of the discretised points in the objective (µ,σ 2 ) space.

Figure 8 :

 8 Figure 8: Test-case 1: Cost comparison between the A Priori MetaModel and Double loop approach.

Figure 9 :

 9 Figure 9: Test-case 1-a: Cost comparison between APMM and three SABBa variants.

Figure 10 :

 10 Figure 10: Worst optimisation results: (a) SABBa SA-CS, 93 evaluations , Q B = 1.29 × 10 -2 . (b) APMM, 100 evaluations , Q B = 1.91 × 10 -1 .

Figure 11 :

 11 Figure 11: Test-case 1-b: Cost comparison between APMM and three SABBa variants.

Figure 12 :

 12 Figure 12: Worst optimisation results: (a) SABBa CS-REU-G, 196 evaluations , Q B = 1.23 × 10 -2 . (b) APMM, 200 evaluations , Q B = 7.94 × 10 -1 .

Figure 13 :

 13 Figure 13: Test-case 2: Cost comparison between APMM and three SABBa variants.

Figure 14 :

 14 Figure 14: Worst optimisation results: (a) APMM, 150 evaluations, Q B = 1.0 × 10 0 . (b) SABBa CS-REU-G, 161 evaluations, Q B = 1.12 × 10 -1 . (c) APMM, 200 evaluations, Q B = 3.24 × 10 -1 .

Figure 15 :

 15 Figure 15: (a) Schematic representation of the two bars, with bar section above (from [37]). (b) Representation in the objective/constraints space, constraints in abscissa and ordinate, objective in color.

Figure 16 :

 16 Figure 16: Optimisation output with POP greyscale in (a) objective space (optimum in red) and (b) constraints space (constraint limit in black).

Figure 17 :

 17 Figure 17: Optimisation output with POP greyscale in (a) objective space (optimum in red) and (b) constraints space (constraint limit in black).

Inria

Figure 18 :

 18 Figure 18: Found optima from the three runs (in black) in the constraints space, with objective value in color.

Figure 19 :

 19 Figure 19: Surface and in-depth temperatures obtained with nominal parameters.

Figure 20 :

 20 Figure 20: Outputs of the SABBa framework in (a) the input space, (b) the objective/constraint space.

Figure 21 :

 21 Figure 21: All constructed boxes, depicted in blue when dominated, in red when entirely in failure zone and in green when non-dominated. (a) Full view, (b) zoomed-in view.

Figure 22 :

 22 Figure 22: Final SA model, representing ρ SA f (x) where the constraint ρ SAg (x) ≤ 473.15 K is satisfied.

Figure 23 :

 23 Figure 23: Computational details

  By definition f (x, •) ∈ [ f x -ε fx , f x + ε fx ]. The monotonicity of the expectation operator then gives E ξ [f (x, ξ)] ∈ E ξ [ f x (ξ) -ε fx (ξ)], E ξ [ f x (ξ) + ε fx (ξ)] , or equivalently µ(x) ∈ µ(x) -E ξ [ε fx (ξ)], µ(x) + E ξ [ε fx (ξ)], thus giving the conservative error found in the above ε µ (x) = E ξ [ε fx (ξ)].

ξ[

  f x (ξ)], m(x) = min ξ [f (x, ξ)] ∈ min ξ [ f - x (ξ)], min ξ [ f + x (ξ)] = m(x) -m(x) -min ξ [ f - x (ξ)] , m(x) + m(x) -min ξ [ f + x (ξ)] ∈ m(x) -ε min (x), m(x) + ε min (x)whereε min (x) = max m(x) -min ξ [ f - x (ξ)] , m(x) -min ξ [ f + x (ξ)]

Table 1 .

 1 ). The computed values of POP are reported in Table 1. boxes are equally close to the Pareto front, whereas box B is sensibly less efficient than B 1 , hence further from the Pareto front. However, one can see in Table 1 that P OP true and P OP av yield higher values of POP for box B 2 than for boxes B 3 to B 6 . P OP av and P OP min are easily computable, but the former can lead to unintuitive probabilities for strictly dominated boxes. This is notably the case for box B 6 in Example 2 and box B 4 in Example 3, as shown again in

			objective 2 (minimisation)	3 4 5 6 7 8	B 1 B 2
				2	3	4	5	6	7	8
					objective 1 (minimisation)
			Example 1			Example 2			Example 3
		P OPtrue P OPav P OPmin P OPtrue P OPav P OPmin P OPtrue P OPav P OPmin
	B1	0.898	0.98	0.898		1.0		1.0	1.0	0.0	0.0	0.0
	B2	0.538	0.908	0.538		0.902		0.98	0.912	0.254	0.279	0.254
	B3	0.521	0.85	0.75		0.39		0.828	0.419	0.724	0.943	0.746
	B4	0.521	0.85	0.75		0.621		0.924	0.622	0.0	0.628	0.0
	B5	0.521	0.85	0.75		0.017		0.508	0.07	0.022	0.483	0.031
	B6	0.521	0.85	0.75		0.0		0.735	0.0	0.0	0.0	0.0

B 3 ,B 4 ,B 5 ,B 6

Figure 4: Example 1, with boxes B 3 to B 6 superimposed.

Best and worst outcomes in green and red respectively. Associated POPs in Table

1

. Figure

4

illustrates the example 1 with clustered boxes (B 3 to B 6 ). Geometrically, all these

Table 1 :

 1 POPs comparison on examples 1 to 3 (Figs. 6 and 7)

RR n°9214

rue Honoré d'Estienne d'Orves Bâtiment Alan Turing Campus de l'École Polytechnique 91120 Palaiseau Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399


Laurent Savy
email: laurent.savy@fr.thalesgroup.com

Mathieu Cattenoz
email: mathieu.cattenoz@gmail.com

Sylvie Marcos
email: sylvie.marcos@l2s.centralesupelec.fr

Robust Orthogonal Matching Pursuit Detector Applications to off-grid target issues in MIMO Radar

Keywords: Orthogonal Matching Pursuit, off-grid target, MIMO radar

The paper proposes to solve the off-grid issue arising in sparse processing involving Orthogonal Matching Pursuit (OMP) algorithm. It starts by establishing the Generalized Likelihood Ratio Test using a stochastic model to represent target location uncertainty within one resolution cell. In a second step, a linear approximate of the detector is built, requiring only a grid oversampling of 2 relative to the resolution cell size. Moreover, it involves only the formation of standard delta channels in each dimension, leading to low computational burden without loss of efficiency. The approach is here demonstrated on MIMO radar detection.

INTRODUCTION

In conventional radar processing, since the target parameters (for instance range, Doppler or angle) are unknown, different matched filter operations (range matched filter, Doppler processing, array beamforming) are applied for different parameter hypotheses in order to retrieve the target. This matched filter is mainly derived in the single target hypothesis and is however used in multi-targets environment with loss of processing. This usually consists in adding standard weighting windows onto the filter in order to lower the side lobes level of the mismatched filtering. In addition, the different single-source assumptions define a grid in the parameter space and the targets are never exactly on the resolution grid, causing a loss of processing due to a shift between the nearest grid point and the real target steering vector. However, for the specific kind of signals encountered in radar, the average losses over all target position within a resolution cell with conventional weighting functions are generally contained below 3dB. Theses losses are traditionally accepted to save computational load and remain in the most manageable case of the multiple single-source assumptions for detection with noise independence assumptions, which also allow the theoretical calculation of the threshold for the control of the false alarm probability.

However, in some situations, it appears that the mismatched filter in the multi-targets environment leads to poor performance and is not practical, due to the inability to find efficient apodization techniques. This is the case when the signal model deviates significantly form the linear phase complex exponential, found in particular in coherent MIMO radars where the range and angle dimensions are coupled due to the non-perfect orthogonality of the waveforms [START_REF] Rabaste | Signal waveforms and range/angle coupling in coherent colocated MIMO radar[END_REF].

In this situation, the main idea is to find a feasible approximate implementation of the global matched filter in the context of sparse targets hypothesis. In presence of spread nonpoint targets, a relevant solution is the Orthogonal Matching Pursuit (OMP) algorithm, see for instance [START_REF] Bosse | Matching pursuit via continuous resolution cell rejection in presence of unresolved radar targets[END_REF]. Because the projector onto the space orthogonal to one target steering vector is very sharp, the standard OMP solution requires to severe oversampling to mitigate the off-grid target issue. Anyway, in case of closely spaced targets (typically less than one resolution cell), estimation of target signal in the detection step becomes so biased that the rejection step can only perform poorly and spurious detections may occur. In [START_REF] Bosse | Matching pursuit via continuous resolution cell rejection in presence of unresolved radar targets[END_REF][START_REF] Rabaste | Off-Grid Target Detection With Normalized Matched Subspace Filter[END_REF], the authors proposed to mitigate the off-grid issue by using in the OMP rejection step the space spanned by the entire resolution cell. Unfortunately, the proposed methods are mainly dedicated to the case of steering vectors being complex exponential with linear phase evolution, otherwise computational load will become prohibitive.

In [START_REF] Cattenoz | Adaptive processing methods for MIMO radar experimental signals[END_REF] we proposed to mitigate the off-grid issue by using in the OMP rejection step the space spanned by adjacent steering vectors of a half cell. The basic idea is to enlarge the rejection region in the OMP rejection step. The method has been proven to work well on real data, to be computationally very efficient, and was justified using the development in Taylor series at second order of the steering vector in [START_REF] Cattenoz | Extended orthogonal matching pursuit for robust and fast target localisation in multiple-input multipleoutput radar[END_REF].

In this paper we go much deeper in the previous justification, showing the link with the Generalized Likelihood Ratio Test (GLRT) for target arbitrary localized within one resolution cell (or extended target) and deriving the discrete formulation of the detector in the Low Energy Coherence (LEC) case, already established for continuous signal in [START_REF] Harry | Detection, Estimation, and Modulation Theory, Part III: Radar-Sonar Signal Processing and Gaussian Signals in Noise[END_REF]. This is detailed in section II.

In section III, using the Taylor series development of the steering vectors in one resolution cell, we show that the space spanned by one resolution cell can be well defined by few steering vectors. It is furthermore shown that the computation of these steering vectors only requires an oversampling by two of the conventional one-cell grid, for the majority of practical radar waveforms and/or radar antenna configurations.

Section IV, illustrates in more details the application of the proposed method to the coherent MIMO radar detection processing, involving OMP algorithm.

II. LOW COHERENCE ENERGY CASE DETECTOR FOR OFF-

GRID TARGET -IMPACT ON OMP REJECTION STEP

A. Standard matched filter and OMP

As a clarification, we recall that the standard canonical hypothesis testing (i.e. on whitened data) for point target leading to the matched filter is the following situation:

) ( : 0 only noise H n x = n ~) , 0 ( Id CN ) ( ) ( : 0 1 noise signal H + + = n θ s x α (1)
α is the unknown target complex amplitude and ) ( 0 θ s is the steering vector, depending on a vector of parameters.

[

]

T P 0 02 01 0 θ θ θ L = θ
. 0 θ =(range, and/or Doppler frequency, and/or angles) in radar detection problems.

The GLRT for detection is then the matched filter:

) ( ) ) ) ) ( 0 1 0 0 2 0 0 threshold H H Y H H η η < > = s(θ (θ s x (θ s θ (2)
0 θ is usually sampled on a grid with mesh corresponding to a resolution cell. The OMP algorithm built on the previous model is then an iterative implementation of the multi-target matched filter that consists in detecting and rejecting each contribution of the strongest targets successively. It starts with x x = 1 and the k th iteration of the procedure can be summarized by the two following steps:

[ ] [ ] x A A A A I x A θ s A A s(θ (θ s x (θ s θ θ H k k H k k k k k k H k H k 1 1 0 0 1 0 0 2 0 0 ) ( ( , ) ( , : step Rejection ) ) ) Argmax find : step Detection 0 - + - - = = = = (3)

B. Low Coherency Energy detector and impact on OMP

In order to deal with the off-grid issue, we first change the target deterministic model in hypothesis H 1 to a stochastic one as follow:

n t x + = : 1 H t ~)) ( , 0 ( 0 2 θ S σ CN , n ~) , 0 ( Id CN ∫ ∫ ∆ + ∆ - ∆ + ∆ - θ θ θ ∆ ∆ ∆ = L 2 1 ∆θ around [ ] T P 0 02 01 0 θ θ θ L = θ .
2 σ is the unknown target power. It is worth noting that this model can also be used for extended target detection or for modeling steering vectors uncertainties (by adding extra parameters in steering vector model in the last case).

The likelihood ratio test with this stochastic model can be written now after some operations:

) ( ) )) ( ( ( 0 2 1 0 2 0 1 θ S Id x θ S Id Id x σ η σ + < > + - - Log H H H (5)
Denoting the singular value decomposition of ) ( 0 θ S by:

∑ ∑ ∑ + = + + - = + = = - n n n H n n n n n H n n n Log Log s NxN N 1 n ) 1 ( ) ( 1 )) ( ( ) ( of ize : : ) ( 2 0 2 2 2 1 0 2 0 0 λ σ σ λ σ λ σ σ λ θ S Id v v Id θ S Id θ S v v θ S (6)
After some more operations, the likelihood ratio test ( 5) can be finally written as :

∑ ∑ + < > + n n H H H n n n n Log ) 1 ( 1 2 2 2 2 0 1 λ σ η λ σ λ σ x v (7)
It is worth noting that ( 7) is a weighted and generalized version of the subspace-based matched detector [START_REF] Harry | Detection, Estimation, and Modulation Theory, Part III: Radar-Sonar Signal Processing and Gaussian Signals in Noise[END_REF], where projections on the eigenvectors of ) ( 0 θ S accounts for signal to noise ratio in each eigenvector related subspace. Only in the case of a strong target and a limited number of non-zero eigenvalues of ) ( 0 θ S we recover the matched subspace detector of [START_REF] Harry | Detection, Estimation, and Modulation Theory, Part III: Radar-Sonar Signal Processing and Gaussian Signals in Noise[END_REF]. One can also note that we recover the matched filter of (2) in case ) ( 0 θ S is of rank 1 and in that case (4) is simply written as

) ( ) ( ) ( 0 0 0 θ s θ s θ S H = H 0 0 v v = .
More interesting is the case of weak target (much below the noise), where we can get rid of the target power 2 σ in the test (7). Indeed, for weak target the test leads to:

{ } n n n n n H H H n n H n n n Log and Tr λ σ λ σ λ σ λ σ λ σ η λ λ 2 2 2 2 2 0 0 2 ) 1 ( 1 because ) ( ) ( 0 1 ≈ + ≈ + < > = ∑ ∑ θ S x θ S x x v (8)
Equation [START_REF] Chen | MIMO Radar Ambiguity Properties and Optimization Using Frequency-Hopping Waveforms[END_REF] shows that for weak target the detector simply turns to the weighted average of the projections onto the eigenvectors of ) ( 0 θ S by the eigenvalues. In our situation of interest, where

) ( 0 θ S
describes the signal in one resolution cell, the s n ' λ in (6) decrease very rapidly.

Equation [START_REF] Chen | MIMO Radar Ambiguity Properties and Optimization Using Frequency-Hopping Waveforms[END_REF] then can be reduced with little loss to the first term in the sum for detection purpose. However, for enough rejection of signal falling in one resolution, as required for OMP algorithm rejection step, a few more eigenvectors are needed. The OMP procedure with the stochastic target model then becomes:

[ ] x A A A A I x V A A s(θ (θ s x (θ s θ θ ) ) ( ( , : step Rejection ) ) ) Argmax find : step Detection 1 1 1 0 0 2 0 0 0 H k k H k k k k s k k H k H k - + - - = = = (9) [ ] k M k k s 1 0 , , - = v v V L : matrix

C. Illustration on mono dimensional complex exponential

signal with linear phase The previous results are illustrated for the particular but relevant case of a mono dimensional complex exponential with linear phase for signal ) (θ s , which can represent in radar Doppler analysis, range compression, or antenna beamforming on a uniform linear array. In that case: In that particular case, the eigenvalues of ) ( 0 θ S are the same than the ones of ) 0 ( S , and the n th eigenvector

T N N i N n N i N N i e e e θ       = - - - - - 2 ) 1 ( 2 ) 2 ) 1 ( ( 2 2
) ( 0 θ n v of ) ( 0 θ S is equal to ) ( 0 θ s ) 0 ( n v , with ) 0 ( n v
the n th eigenvector of ) 0 ( S . The eigenvalues of ) 0 ( S are plotted for illustration on Fig. 1. It is clear that most of the energy is captured by the first eigenvector with 80% (-1dB) of the total energy (see zoom on the first eigenvalues).

In next section, using 1 st order and 2 nd order expansion of the steering vector in one resolution cell, we show that the space spanned by one resolution cell can be well defined by only two or three steering vectors for each parameter. It leads to a very computationally efficient implementation of the projection of the signal in the space orthogonal to one resolution cell, needing only an oversampling by 2 of the standard one-cell stepped grid.

A. Taylor series expansion approximate

We start here by 1 st and 2 nd order linear approximation of the steering vector, quite valid in one resolution cell, around the value 0 θ (the center of the resolution cell): In order to compute very efficiently the delta channels x V H s (and hence the rejection step), we approximate them by finite differences derived from a twice oversampled grid to get half resolution cells. For a mono-dimensional signal, as the one in (10), this approximation is written classically as:

order 2 ) ( ) ( ) ( 2 1 order 1 ) ( ) ( ) ( ) (
4 / ) 2 ( ) ( 2 ) 2 ( ) ( ) 2 ( ) 2 ( ) ( 2 0 0 0 θ θ θ ∆ ∆ - + - ∆ + ≈ ∂ ∂ ∆ ∆ - - ∆ + ≈ ∂ ∂ = = s s s s s s s (14)
To compare the efficiency of the projections defined by ( 13) and its approximates defined by ( 14), we took the same signal example as in Fig. 1, and computed the projector response

) ( ) ) ( )( ( 1 θ θ s V V V V Id s H s s H s s H - -
, and added also the standard OMP projector using only

0 s V = s .
The results are plotted on Fig. 2. Clearly the delta channel projectors avoid the sharpness of the standard projector, while preserving a deep notch. It can be seen also that the approximate delta channels of (12) are closed to the exact delta channels of (11), which validates the proposed approach.

It is worth noting that at 2 nd order we get the same projector using the delta channels defined by (14) in s V or using directly

      ∆ - ∆ + = ) 2 ( ) 2 ( 0 0 0 θ θ θ θ s s s V s (15)
Indeed the delta channels are a linear combination of the three vectors involved in (15). This latest formulation is the half-cell projector already identified in an heuristic way in [START_REF] Cattenoz | Adaptive processing methods for MIMO radar experimental signals[END_REF].

IV. APPLICATION TO MIMO RADAR PROCESSING

In this section we will first briefly review the coherent MIMO radar with intrapulse coding and the associated sidelobes issues of the ambiguity function (to go further, detailed discussions can be found for instance in [START_REF] Rabaste | Signal waveforms and range/angle coupling in coherent colocated MIMO radar[END_REF]). In a second step, we will then give some applications of the technique described in III on simulated and real MIMO data.

A. Brief review of coherent MIMO radar processing and ambiguity function

For the sake of simplicity, let us consider a transmitting linear array composed of M colocated isotropic antennas. In a conventional phased array, all antennas transmit the same waveform, and only the phases of the waveforms are different to form the beam in a given direction. In coherent colocated intrapulse coding MIMO radar, each of the M antennas transmits a different pulsed waveform ) (t s m , with repetition interval T R so that the signal transmitted by the array in target located at direction c θ is given by:

∑ ∑ = - = l l M m R m i c E T t s e t s m E c T 1 ) ( ) ( ) , ( , x k θ θ (16)
where

m E,
x is a vector representing the position of the m th antenna, and

) ( c θ k is the wave vector.
Let us now assume without lack of generality that the receive array is composed of one isotropic antenna. Indeed it is shown in [START_REF] Chen | MIMO Radar Ambiguity Properties and Optimization Using Frequency-Hopping Waveforms[END_REF] that the angular reception processing can be completely decoupled. Then the received signal from a target of unknown amplitude α, with delay c r c / 2 = τ (r is target range) and Doppler frequency c υ is provided by:

) ( ) , , , ( ) ( t b t s t x c c c + = υ θ τ α with the noise ∑ ∑ = - - = l l M m t i c R m i c c c c m E c T e T t s e t s 1 2 ) ( ) ( ) , , , ( , πυ θ τ υ θ τ x k
Then, the response of the matched filter at the output of the receiver is:

∑ ∑ ∫ = - - - - = ' 1 ' 2 * ' ) ( ) ( ) ( ) , , ( ' , l l M m t i R m i e T t s t x e y m E T πυ θ τ υ θ τ x k
Moreover it is also shown in [START_REF] Chen | MIMO Radar Ambiguity Properties and Optimization Using Frequency-Hopping Waveforms[END_REF] that in case of no Doppler effect on the pulse, the Doppler processing can be completely decoupled from the range-angle processing, and that what matters in coherent MIMO processing is the range-angle matched filter for one elementary pulse given by:

dt t s t x e y M m m i m E T ∑ ∫ θ τ θ τ H w H s H y = = =
And the corresponding ambiguity function is: and then a decoupled range-angle ambiguity function (see [START_REF] Rabaste | Signal waveforms and range/angle coupling in coherent colocated MIMO radar[END_REF] for details):

( ) ) ( ) ( ) ( ) ( ) , , (
) ( ) ( ) ( ) , ( c ss c s H s c c decoupled C τ τ θ θ θ θ τ τ - ⋅ = - - Χ s s
However, some waveform families may approach this condition. Among these families, an interesting case is provided by the Code Division Multiple Access (CDMA) phase codes used in digital communications [START_REF] Rabaste | Signal waveforms and range/angle coupling in coherent colocated MIMO radar[END_REF]. For instance Gold codes which are optimized for having good periodic auto and cross correlation also provide good aperiodic correlations. Fig. 3 presents the range/angle 2D cut of the ambiguity function for τ c =0 and θ c =0 with Gold codes with N=127 chips, M=12 transmitters, pulse duration Tp=100µs, and bandwidth B=2 MHz. This set of parameters will be used further throughout the rest of the paper. We observe a thumbtack type MIMO ambiguity function instead of a 2D sinc function in the ideal decoupled case. There is no noticeable coupling effect (so we get the full resolution in range and angle) but the energy is spread over the entire parameter space with high sidelobe level. This spreading comes from the summation of the non-zero cross correlations of the codes. This spreading can moreover not be mitigated with the use of conventional techniques of apodization windows. There is therefore a strong need to decrease this level in multitarget environment to avoid the masking of weak targets by the sidelobes of strong targets. In the next paragraph, we present benefits provided by the robust OMP technique detailed in section III for that goal.

B. Extension of Robust OMP to MIMO radar

First we have to extend the approach of finite differences approximation defined in (14) for 1D on a 2D half-cell mesh grid. It was found, due to the finite differences approximation, that sufficient nulling of the space spanned by one cell requires not only to compute the regular 1 st and 2 nd order differences on each axes (4 parameters), but also the 1 st and 2 nd order differences obtained after a rotation of 45° of the parameter plane (4 additional parameters). This is equivalent to take for s V in (13) the cell under test and its eight neighbors on a half-cell mesh grid. This is the direct extension of (14) in 1D (2 neighbors) to 2D (eight neighbors). The rationale is described in detail in [START_REF] Cattenoz | Extended orthogonal matching pursuit for robust and fast target localisation in multiple-input multipleoutput radar[END_REF]. The extended projector response, ) , ( )

) ( )( , ( 1 θ τ θ τ s V V V V Id s H s s H s s H - - , and

C. Results on simulated data of Robust OMP for MIMO radar detection

In order to test the effectiveness of the robust OMP method, i.e. how the filter manages difficult situations, we use a particular scenario with two strong targets (same output SNR of 40 dB) with a decreasing distance d from 1.5 to 0.25 with steps of 0.25 range cell distances, and with a third weak target hidden in the sidelobes of the matched filter (output SNR of 15dB), located on a grid node in (-300m ; 0). The left-hand strong target is fixed on a grid node, while the right-hand strong target has a changing position, localized on a node only when its distance to the left target is a multiple of 0.5 cell width. We compute for each scenario the matched filter output, the standard OMP output and the proposed robust OMP output. These results are displayed on Fig. 5, Fig. 6, and Fig. 7 respectively. The proposed robust OMP avoids the spurious generation with standard OMP in case of closed or off-grid targets, and behaves as the matched filter for closed target within one resolution cell, giving only one detection. It also recovers well the weak target hidden in the sidelobes of the matched filter. In [START_REF] Cattenoz | Extended orthogonal matching pursuit for robust and fast target localisation in multiple-input multipleoutput radar[END_REF] the approach is also demonstrated on real data. In this paper we derived a robust version of OMP algorithm against off-grid/closed targets issues, justified with a stochastic model on targets relative to their location within one resolution cell in the parameter test grid. We have proposed a low computational load implementation involving only a half-cell mesh oversampled grid. The approach was demonstrated on simulated MIMO radar data. In [START_REF] Cattenoz | Extended orthogonal matching pursuit for robust and fast target localisation in multiple-input multipleoutput radar[END_REF] results on real data are given and it is also shown that the proposed method can be extended to model mismatch arising from hardware imperfect knowledge/calibration by incorporating small variants of the steering vector in the computation of (4).

  formed by the eigenvectors corresponding to the M dominant eigenvalues of ) ( 0 k θ S that are sufficient to get enough cancelation of all the signals within one resolution cell.

n

  can represent either the pulse index in Doppler analysis, the frequency bin index in pulse compression, or the antenna element index in a linear uniformly space array. θ represents then respectively the Doppler bin index, the range bin index, or the angle bin index. The resolution cell extent is defined by 1 : Hadamard product.

Fig. 1

 1 Fig.1Eigenvalues of S(0) for complex exponential signal with N=32

  and to use, for rejection step in OMP, the matrix:

Fig. 2

 2 Fig. 2 Projectors comparison

  MxM is the time cross-correlation matrix of the waveforms, and is not a diagonal matrix with practical waveforms, inducing rangeangle coupling in coherent MIMO radar. Note also that the ambiguity function is a 3D function, depending on tested angle θ and actual target angle θ c . Only hypothetic unrealizable orthogonal waveforms, with same autocorrelation function

Fig. 3 -

 3 Fig. 3 -Range/angle cut of the MIMO ambiguity function for CDMA Gold codes of length N=127, N=12 transmitters, Tp=100µs, B=2 MHz

4 .

 4 Clearly the notch covers one resolution cell with a sufficient depth of almost 60dB.

Fig. 4 -

 4 Fig. 4 -2D standard (left) and extended (right) projector responses

Fig. 5 -

 5 Fig. 5 -Matched Filter Output on simulated data

Fig. 7

 7 Fig. 7 Robust OMP Output on simulated data