Bijan Mohammadi
email: bijan.mohammadi@um2.fr

Ensemble Kalman Filters (EnKF) and geometric characterization of sensitivity spaces for Uncertainty Quantification in optimization

Keywords: EnKF, Uncertainty Quantification, principal angles, inverse problems, history matching, extreme scenarios

We present an original framework for uncertainty quantification (UQ) in optimization. It is based on a cascade of ingredients with growing computational complexity for both forward and reverse uncertainty propagation. The approach is merely geometric. It starts with a complexity-based splitting of the independent variables and the definition of a parametric optimization problem. Geometric characterization of global sensitivity spaces through their dimensions and relative positions by the principal angles between global search subspaces bring a first set of information on the impact of uncertainties on the functioning parameters on the optimal solution. Joining the multi-point descent direction and the quantiles on the optimization parameters permits to define the notion of Directional Extreme Scenarios (DES) without sampling of large dimension design spaces. One goes beyond DES with Ensemble Kalman Filters (EnKF) after the multi-point optimization algorithm is cast into an ensemble simulation environment. This formulation accounts for the variability in large dimension. The UQ cascade ends with the joint application of the EnKF and DES leading to the concept of Ensemble Directional Extreme Scenarios (EDES) which provides more exhaustive possible extreme scenarios knowing the Probability Density Function of our optimization parameters. A final interest of the approach is that it provides an indication of the size of the ensemble which must be considered in the EnKF. These ingredients are illustrated on an history matching problem.

INTRODUCTION

Forward and backward uncertainty propagation [START_REF]AIAA Guide for the Verification and Validation of Computational Fluid Dynamics Simulations[END_REF][START_REF] Caers | Modeling Uncertainty in the Earth Sciences[END_REF][START_REF] Kapur | Reliability in Engineering Design[END_REF] are obviously of great importance. For instance, taking into account manufacturing uncertainties is of utmost importance for a design procedure to be efficient as it is impossible to make sure that the final product will exactly correspond to the design specifications. In shape optimization, for instance, either in mono or multidisciplinary frameworks, this uncertainty is rarely accounted for [START_REF] Mohammadi | Shape Optimization in Fluid Mechanics[END_REF][START_REF] Martins | Multidisciplinary design optimization: A Survey of architectures[END_REF].

This question can be formulated through both forward and backward uncertainty propagation. Indeed, manufacturing uncertainties can be prescribed through probability density functions of the design parameters (often coming as characteristics of the product manufacturing process) [START_REF] Mohammadi | Value at Risk for confidence level quantifications in robust engineering optimization[END_REF][START_REF] Schillings | Efficient shape optimization for certain and uncertain aerodynamic design[END_REF]. They can also be introduced through uncertainties on the performances of a design (or other observations).

1.1. Context of the work. We consider a generic situation where the simulation aims at predicting a given quantity of interest j (e.g. the maximum value of a variable in a given area) and there are a few functioning u and several control parameters x involved. The ranges of the functioning parameters define the global operating/functioning conditions of a given design. This splitting of the independent variables in two sets is important for the rest of the paper.

The literature on uncertainty quantification (UQ) is huge. In short, in our situation, forward propagation aims at defining a probability density function for j knowing those of x and u [START_REF] Hoel | Introduction to Mathematical Statistics[END_REF][START_REF] Lindman | Analysis of Variance in Complex Experimental Designs[END_REF][START_REF] Spooner | Stable Adaptive Control and Estimation for Nonlinear Systems: Neural and Fuzzy Approximator Techniques[END_REF]. This can be done, for instance, through Monte Carlo simulations or a separation between deterministic and stochastic features using Karhunen-Loeve theory (polynomial chaos theory belongs to this class) [START_REF] Ghanem | On the construction and analysis of stochastic models: characterization and propagation of the errors associated with limited data[END_REF][START_REF] Iaccarino | Quantification of Uncertainty in Flow Simulations Using Probabilistic Methods[END_REF][START_REF] Xiu | Numerical Methods for Stochastic Computations: A Spectral Method Approach[END_REF][START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF][START_REF] Karniadakis | Multi-element generalized polynomial chaos for arbitrary probability measures[END_REF]. Examples of shape optimization with polynomial chaos and surrogate models during optimization to address the issue of functional evaluations are given in [START_REF] Congedoa | Shape Optimization of an Airfoil in a BZT Flow with Multiple-source Uncertainties[END_REF][START_REF] Cinnella | Robust optimization of dense gas flows under uncertain operating conditions[END_REF].

Backward propagation aims at reducing models bias or calibrating models parameters knowing the probability density function of j (or other constraints and observations) [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF][START_REF] Ide | Unified notation for data assimilation: operational, sequential and variational[END_REF][START_REF] Casella | Statistical Inference[END_REF]. This can be seen as a minimization problem and Kalman filters [START_REF] Kalman | A New Approach to Linear Filtering and Prediction Problems[END_REF] give, for instance, an elegant framework for this inversion assimilating the uncertainties on the observations.

Our aim is to propose a geometric framework to address, in our particular situation, the curse of dimensionality of existing approaches related to the explosion of their computational complexity due to the sampling necessary to access probabilistic information (momentum), even if this can be improved with intelligent sampling techniques [START_REF] Smolyak | Quadrature and interpolation formulas for tensor products of certain classes of functions[END_REF][START_REF] Bungartz | Sparse Grids[END_REF]. The different ingredients presented here can be applied with either high-fidelity or reduced order models, when available [START_REF] Schilders | Model order reduction: Theory, research aspects and applications[END_REF][START_REF] Obinata | Model reduction for control system design[END_REF][START_REF] Qu | Model Order Reduction Techniques with Applications in Finite Element Analysis[END_REF][START_REF] Veroy | Certified real-time solution of the parametrized steady incompressible Navier-Stokes equations: Rigorous reduced-basis a posteriori error bounds[END_REF]. Loworder models are often used instead of the full models to overcome the computational complexity of UQ.

After the splitting of the independent variables mentioned above, we define a multi-point formulation to account for the variability on u. This is feasible because the size of u is assumed small. We define a global sensitivity space using the sensitivities of j with respect to x for the multi-point problem. Once this space built, we analyze the dimension of its free generator subspace. Previous works have shown how to perform this task and how to use this information for adaptive sampling and robust optimization [START_REF] Mohammadi | Reduced sampling and incomplete sensitivity for low-complexity robust parametric optimization[END_REF][START_REF] Gallard | An adaptive multipoint formulation for robust parametric optimization[END_REF].

The next step is to analyze the impact of different modeling or discretizations on the results. Different models or solution procedures lead to different sensitivity spaces. Beyond their respective dimensions, principal angles [START_REF] Jordan | Essay on geometry in n dimensions[END_REF][START_REF] Gluck | Great circle fibrations of the three-sphere[END_REF][START_REF] Jiang | Angles between Euclidean subspaces[END_REF][START_REF] Shonkwiler | Poincare duality angles for Riemannian manifolds with boundary[END_REF] between the respective sensitivity spaces permit to measure the deviation due such changes. The dimensions of the spaces and the angles are interesting measures for both the epistemic and aleatory uncertainties. Indeed, suppose that, at given modeling procedure, the dimensions of the sensitivity spaces remain unchanged when enriching the sampling of the functioning parameter range, this stability would be a first indication of a low level of sensitivity of the simulations with respect to this parameter. Once this is established, principal angles between subspaces permit to analyze both the impact of a given evolution of the modeling on the sensitivity spaces or an enrichment of our sampling. Eventually, constant dimension and low angles will clearly indicate a situation of low uncertainty.

These ingredients can be used in a context of multi-point robust analysis of a system to define worst-case scenarios for its functioning. To this end we combine a multi-point sensitivity with the probabilistic features of the control parameters through their quantiles [START_REF] Mohammadi | Value at Risk for confidence level quantifications in robust engineering optimization[END_REF][START_REF] Ph | Value at Risk: The New Benchmark for Managing Financial Risk[END_REF]. These ingredients permit to define the concept of Directional Extreme Scenarios (DES) without a sampling of large dimension design spaces.

Ensemble Kalman filters (EnKF) [START_REF] Kalman | A New Approach to Linear Filtering and Prediction Problems[END_REF][START_REF] Anderson | Optimal Filtering[END_REF][START_REF] Jazwinski | Stochastic Processes and Filtering Theory[END_REF][START_REF] Gelb | Applied Optimal Estimation[END_REF][START_REF] Evensen | Advanced Data Assimilation for Strongly Nonlinear Dynamics[END_REF][START_REF] Evensen | Sequential Data Assimilation for Nonlinear Dynamics: The Ensemble Kalman Filter In Ocean Forecasting: Conceptual basis and applications[END_REF] permit to go beyond the directional uncertainty quantification concept when accounting for the uncertainties in large dimension. It also permits backward uncertainty propagation assimilating the uncertainty on the functional and constraints during the design. We propose to cast our multi-point optimization problem into the ensemble formulation. The geometric uncertainty quantification construction ends with the joint application of the EnKF and DES leading to the concept of Ensemble Directional Extreme Scenarios (EDES) which provides exhaustive possible extreme scenarios knowing the probability density function of the optimization parameters and this without any sampling of a large dimensional parameter space. Finally, the geometric information on the global sensitivity spaces provides lower bounds on the size of the ensemble which must be considered in the EnKF.

Summary of this work ingredients and calculation complexity.

To summarize this work will propose a cascade of ingredients to account for uncertainties avoiding any sampling of large dimensional spaces. We insist on the fact a sampling is only necessary for the functioning parameters u range (see section 2) leading to a multi-point optimization problem with a linear complexity with respect of the size of the sampling. One originality is the definition of the direction d in section 2.1. d is also used, together with quantiles from the PDF of x, to define quantile-based extreme scenarios in section 3. In sections 3 and 4 geometric analysis of the sensitivity space S m (space dimension, principal angles) gives some information on the impact of the variability of x on the results without any extra functional evaluation. To go beyond directional extreme scenarios and better account for the variability of x, the previous ingredients are cast into the ensemble Kalman filter framework in section 5. Finally, the geometric analysis ingredients are used for each of the ensemble run to introduce Ensemble Directional Extreme Scenarios in section 5.3. Concerning the computational cost of these analysis, one can say that, when using the same calculation ingredients that for a high-fidelity simulation (i.e. without calling for low-order models or cheaper discretizations), the best calculation complexity one might think of for a simulation under uncertainty is when its cost is comparable to the deterministic situation. This is clearly unreachable except if all the extra effort can be achieved in a fully parallel manner and parallel to the initial deterministic calculation. This permits for the time to solution to remain unchanged for a simulation when accounting for the presence of uncertainties. This is obviously the case with Monte Carlo approaches. But these are quite expensive and do not take advantage of available simulation environments. In particular, when an adjoint-based optimization environment exists. Hence, we propose an original way to upgrade existing platform without abandoning what has been built for the deterministic situations and with keeping the time to solution unchanged in the presence of uncertainties with two sources of parallelism coming from the multi-point formulation to account for the uncertainties on the functioning parameters and from the EnKF formulation for those on the optimization variables and observation data. In this work, m being the size of the sampling for the range of u introduced in section 2 and q the ensemble size introduced in section 5, the maximum calculation complexity (to build the Ensemble Directional Extreme Scenarios) at each iteration of optimization is: 2mq direct and adjoint solutions of the state equation. This is especially interesting if n, the size of the control parameter x, is large compared to q.

Parametric optimization

We are interested in a class of optimization problems where the cost function involves a functioning parameter u not considered as a design parameter: [START_REF]AIAA Guide for the Verification and Validation of Computational Fluid Dynamics Simulations[END_REF] min

x∈O ad j(x, u), u ∈ I ⊂ IR α , O ad ⊂ IR n .
where x is the design vector belonging to O ad the optimization admissible domain.

Usually, the functioning parameters u are just a few. On the other hand, the size n of x is usually large. Together, x and u fully describe our system. This splitting between functioning parameters (or operating conditions) and design variables is central to our discussion.

In [START_REF] Mohammadi | Reduced sampling and incomplete sensitivity for low-complexity robust parametric optimization[END_REF][START_REF] Mohammadi | Principal angles between subspaces and reduced order modeling accuracy in optimization[END_REF] we showed how to use multi-point optimization to address such optimization problem. The aim is to remove, during the optimization, the dependency in u. This is done minimizing a functional J(x) encapsulating this dependency expressed through A = {j(x, u l ), u l ∈ I m } over I m a given sampling of I:

(2) J = J(A), such that G(A) 0.

Several choices are possible for J and G to address the issue of robust design. For instance, following Taguchi's definition, one can look for minimal-variance design or only a given level for the variance. Indeed, a classical approach to extend single point design and improve off-design points is to control mean performance and variance of the functional [START_REF] Tang | Uncertainty based robust optimization method for drag minimization problems in aerodynamics[END_REF]. One can also look for information about the tails of the distributions which can be linked to the variance in the Gaussian framework (we will use this relationship in section 3 in quantile-based extreme scenarios).

To control the first two moments one can choose, for instance, J = µ and G = σ -σ 0 with

(3) µ = 1 m u l ∈Im j(x, u l ).
For the sake of simplicity, we consider a uniform weighting in the sum but more sophisticated choices are possible and in particular choices which enforce the performance we eventually want for the design: constant performance over the functioning parameters ranges for a robust design [START_REF] Mohammadi | Reduced sampling and incomplete sensitivity for low-complexity robust parametric optimization[END_REF]. σ 0 is the level of variability of the initial design or of a reference configuration. We proceed as in First-Order Second Moment (FOSM) methods [START_REF] Melchers | Structural Reliability Analysis and Prediction[END_REF]. σ measures the deviation between j(x, u l ):

(4) σ = 1 m -1 u l ∈Im (j(x, u l ) -µ) 2 .
We use this definition when the number of functioning parameters α is large. However, as we want the final performance to be as regular as possible with respect to u, we ask, when possible, for σ to also monitor the regularity of the performance defined for instance adding the following term to (4):

(5) 1 2

u l ∈Im < ∇ u j(x, u l ), ∇ u j(x, u l ) >,
where <, > is the Euclidean scalar product. When α is small, ∇ u j is evaluated by finite differences on the finite difference grid I m . This is reasonable because u is of small size and that we assume uniform PDF for u giving a Cartesian grid for a uniform sampling of the ranges of variation of all the components of u. Such minimization problems have brought new interest to descent methods and this not only because of their lower computational complexity than the gradient free methods [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF][START_REF] Jahn | Vector Optimization: Theory, Applications and Extensions[END_REF][START_REF] Nocedal | Numerical Optimization[END_REF]. Indeed, we will see that gradients are useful to see what should actually be the search space in a context of robust multi-point design. More precisely, the definition of the suitable search space is only possible after identification of the maximal free generator subspace in the space S m generated by the gradients of the functional at the sampling points I m (see section 3). This is the subspace generated by the orthonormal basis identified by the Gram-Schmidt procedure at given precision [START_REF] Mohammadi | Principal angles between subspaces and reduced order modeling accuracy in optimization[END_REF][START_REF] Mohammadi | Uncertainty Quantification by geometric characterization of sensitivity spaces[END_REF]. Indeed, without other information and considering vector spaces theory, the size of the sampling should be larger than the dimension of the control space (i.e. m = n + 1) . We noted, however, [START_REF] Gallard | An adaptive multipoint formulation for robust parametric optimization[END_REF][START_REF] Mohammadi | Reduced sampling and incomplete sensitivity for low-complexity robust parametric optimization[END_REF] that the dimension of the maximal free search subspace is usually much smaller than n.

2.1.

A multi-point descent algorithm. We consider an iterative descent algorithm for our constrained minimization problem involving a direct simulation chain linking the parameters (x, u) to the state y solution of a state equation y = F -1 (x, u) and to a functional j. We assume the functional is enough regular with continuous gradient in both u and x. When cast into the Kalman filters framework in section 5, y (or some part of it) will be called observation. More sophisticated minimization algorithms can be used instead of a simple descent method.

                                         Given x 0 , 0 < ρ, 0 < η 1, I m , k max , 0 < ε 1, optimization iterations k = 1, ..., k max                                  -m parallel state equation solutions y l = F -1 (x k , u l ), u l ∈ I m , -m parallel evaluations of j(x k , y(x k , u l )), u l ∈ I m , -m parallel solutions of the adjoint state V equation: V t F y (y l ) = j t y , u l ∈ I m , -m parallel evaluations of ∇ x j(x k , u l ) = j x + (V t F x ) t , u l ∈ I m , -define d the descent direction: d(µ(x k ), σ(x k )) = ∇ x µ-< ∇ x µ, ∇ x σ > ∇ x σ + η∇ x σ, -control parameter variation, x k+1 = x k -ρd, Stop if d ε,
where a = a/ a is the normalized vector a and ( 6)

∇ x µ = 1 m u l ∈Im ∇ x j(x k , u l ),
and following the definition chosen for σ we have (7)

∇ x σ = 2 m -1 u l ∈Im (j(x k , u l ) -µ)(∇ x j(x k , u l ) -∇ x µ),
with possibly a second term if the regularity issue is included in the definition of σ:

u l ∈Im < ∇ u j(x k , u l ), ∇ ux j(x k , u l ) >,
where ∇ xu j(x k , u l ) is obtained from ∇ x j(x k , u l ) by finite differences on I m , component by component, following what has been done for ∇ u j.

The definition of the descent direction permits to make sure that both µ and σ decrease for small descent steps. Indeed, a first order development in x gives:

σ(x k+1 ) -σ(x k ) = ∇ x σ ∇ x σ.(x k+1 -x k ) = -ρη ∇ x σ 0,
and we have

µ(x k+1 ) -µ(x k ) = ∇ x µ ∇ x µ.(x k+1 -x k ) = -ρ ∇ x µ (1 -ζ 2 + ηζ),
where ζ =< ∇ x µ, ∇ x σ >. Therefore, µ is also decreasing as 1 -ζ 2 + ηζ 0 for |ζ| 1 and 0 < η 1 as chosen in the algorithm. The analysis above can be extended to other functions in [START_REF] Caers | Modeling Uncertainty in the Earth Sciences[END_REF], other than µ and σ, with

d(J, G) = ∇ x J-< ∇ x J, ∇ x G > ∇ x G + η∇ x G
, as far as J and G are such that the quantities in this formula can be defined which implies that we need the functional and constraint to be differentiable.

Geometrical characterization of the global search space

During the iterations of the algorithm above a global sensitivity space can be defined using the gradients ∇ x j(x, u l ) ∈ IR n for u l ∈ I m (I m being a m-point sampling of I):

S m = Span{∇ x j(x, u l ), u l=1,...,m ∈ I m } ⊂ IR n×m .
A major quantity of interest is the dimensions p of this subspace as it can be used to identify if a design is robust or not with respect to u. This dimension depends obviously on the sampling I m of the range of variation of u. But it also depends on the modeling (the governing equations) and the way the state equations are solved (the numerical schemes). Each situation provides a different global sensitivity space. p is therefore also interesting to quantify the impact of a change in the modelling or discretization on a design, also called epistemic uncertainty. Indeed, p << m means that u will have low impact on the sensitivity space S m for the design x k . Of course, this is interesting if the dependency between the functional and the parameters involves a state equation. To illustrate this issue, consider a scalar u, a constant vector ψ and a functional j(x, u) = ψ t x + u. Here, p = 1 is independent of m since ∇ x j(x, u) is independent of u.

In the past we showed how to evaluate p using an incomplete Gram-Schmidt orthonormalization algorithm. This algorithm is parametric with a threshold to define the accuracy of the orthonormalization. One notice that new independent directions appear increasing this precision. This corresponds to the detection of smaller scales in the gradients which are filtered by a less accurate orthonormalization [START_REF] Mohammadi | Reduced sampling and incomplete sensitivity for low-complexity robust parametric optimization[END_REF][START_REF] Mohammadi | Principal angles between subspaces and reduced order modeling accuracy in optimization[END_REF][START_REF] Mohammadi | Uncertainty Quantification by geometric characterization of sensitivity spaces[END_REF].

The information on the dimension of S m for different modelling and/or discretizations can be completed by the respective positions of the different subspaces through the principal angles between them. These two sets of information permit to see uncertainty quantification from a geometric and quite deterministic point of view [START_REF] Mohammadi | Principal angles between subspaces and reduced order modeling accuracy in optimization[END_REF].

Let us briefly recall how to practically compute the principal angles [START_REF] Jiang | Angles between Euclidean subspaces[END_REF][START_REF] Shonkwiler | Poincare duality angles for Riemannian manifolds with boundary[END_REF][START_REF] Mohammadi | Principal angles between subspaces and reduced order modeling accuracy in optimization[END_REF].

θ i ∈ [0, π/2] [33,
For simplicity, suppose A and B are two subspaces of dimension k of IR n , n 2k, although this is not a prerequisite to define the principal angles. Now, let {a i , i = 1, ..., k} and {b i , i = 1, ..., k} be two arbitrary orthonormal bases for A and B. Consider M being the matrix of the projection operator P r A of B onto A defined by:

P r A (b i ) = k j=1 < b i , a j > a j , M = (< b i , a j >) i,j .
The principal angles can be linked to this operator through M = GΣH t , where G and H are orthogonal matrices and Σ = diag(cos(θ i )).

As G and H are orthogonal matrices, this is a Singular Vector Decomposition (SVD) of M . G and H are unknown, but we do not need them to get the principal angles θ i . We recall that the columns of G are the left-singular vectors of M and eigenvectors of M M t and the columns of H are the right-singular vectors of M and eigenvectors of M t M . And most important that cos 2 (θ i ) are the eigenvalues of P r t A P r A which writes in matrix form as:

M t M = (GΣH t ) t (GΣH t ) = HΣ 2 H t with Σ 2 = diag(cos 2 (θ i )).
Therefore, to find the principal angles between subspaces A and B, knowing an orthonormal basis in each subspace, one should calculate M and find the eigenvalues of M t M and take the square root of them. This last operation is valid as the angles are between 0 and π/2, and their cosine therefore always positive. In our case, an orthonormal basis in already available after the Gram-Schmidt orthonormalization has been used to identify the dimension of global search space S m .

Directional extreme scenarios

We showed how to introduce some geometric characterizations of the global sensitivity spaces built with the gradients of a functional ∇ x j(x, u l ) for u l ∈ I m a sampling of a functioning parameter variation range. This was feasible as we assumed there are only a few of such parameters in a system. On the other hand, it is not possible to account for the variability in the design parameter x ∈ IR n as 1 n. In [START_REF] Mohammadi | Value at Risk for confidence level quantifications in robust engineering optimization[END_REF] we showed how to introduce quantiles such as the Value at Risk (VaR) in optimization algorithms with the aim of quantifying our confidence on the optimal solution at low complexity without a sampling of the control space. This concept can be used in combination with the direction d of section 2.1 to define the notion of Directional Extreme Scenarios (DES).

We consider that a given control parameter x i=1,...,n represents, with a confidence level of a, the interval [x i + VaR - a , x i + VaR + a ] around x i and VaR - a 0 VaR + a . One expresses then possible deviations from x via two probability density functions (PDF) and the corresponding values at risk. One can make the hypothesis that the upper and lower bounds of the variations are symmetric, in which case VaR - a = -VaR + a . This is typically the case when uncertainties on a parameter follow a Gaussian distribution. Gaussian distributions are interesting as their Values at Risk are explicitly known: VaR 0.99 = 2.33 and VaR 0.95 = 1.65 for N (0, 1) and VaR a (N (0, σ)) = σVaR a (N (0, 1)).

To summarize x + VaR a (x) is a closed domain in IR n around x:

(8) B a (x) = {z ∈ IR n : x i + VaR - a (x i ) z i x i + VaR + a (x i ), i = 1, .., n}.
The knowledge of the level of uncertainties on x is by itself a problem which we do not discuss here. Manufacturing uncertainties are usually known through probability density functions coming as characteristics of the product manufacturing process: realization within margin errors. These PDFs can also be identified by assimilation observation data. Now, we would like to use this information together with the descent direction d given in algorithm 2.1 to define extreme values for the local variability of the functional along the direction d. We make the assumption that the functional is monotonic over B a (x), continuous and bounded. It is therefore uniformly continuous. This guarantees that a small perturbation in the design parameters will have small effect on the performance, something which is very reasonable to assume for any robust design. For the sake of simplicity, on can suppose the functional is convex in B a (x). If the functional does not behave as such it would be very hazardous to pretend providing any robust design.

These assumptions imply that extreme values of the functional over B a (x) are reached on ∂B a (x), the boundary of B a (x).

The intersection of ∂B a (x) and d permits to define two Directional Extreme Scenarios for x as worst-case scenarios: [START_REF] Ghanem | On the construction and analysis of stochastic models: characterization and propagation of the errors associated with limited data[END_REF] (DES) :

x ± = d ∩ ∂B a (x).
A sketch of this construction is given in figure 1. The uncertainty domain is of rectangular shape because we have assumed the components of x are independent (i.e. the correlation matrix is diagonal, see section 5).

In [START_REF] Mohammadi | Value at Risk for confidence level quantifications in robust engineering optimization[END_REF] we showed how this can be linked to the notion of over-solving where it becomes useless to solve accurately near an optimum when the variation in a control parameter falls below the level of the variability for this parameter. In other words, in the presence of uncertainties, one should not distinguish between the points in B a (x). This interpretation should also be used to define what me mean by uniqueness of the solution: one cannot distinguish between the solutions in B a (x). This construction meets the ingredients developed in the Unscent Kalman filter approaches [START_REF] Julier | A New Extension of the Kalman Filter to Nonlinear Systems[END_REF] and the extreme directional scenarios can be seen as two of the σ-points there.

The next step is to apply the geometric analysis of section 3 after the definition of two global sensitivity spaces S ± m at x ± using ∇ x j(x ± , u l ) for u l ∈ I m following the same procedure than for S m . The respective dimensions p ± of these spaces and their respective position through the principal angles between them give precious information on the impact of local variability on x on the global search space [START_REF] Mohammadi | Principal angles between subspaces and reduced order modeling accuracy in optimization[END_REF][START_REF] Mohammadi | Uncertainty Quantification by geometric characterization of sensitivity spaces[END_REF]. For instance, if the dimensions are the same and angles small, this would mean that variability in x will have very little impact on the design. Also, from the calculation complexity point of view, it is important that these have been built without a sampling of our parameter space of dimension n for x.

x d X+VaRa +/-along d Ba(x) Dq=span{d l ,l=1,…,q} Ba(x) d 1 d q X ) ( VaR - / x D q + ∩ Figure 1
. Sketch of Directional Extreme Scenarios (DES) given by x ± = d∩∂B a (x) and Ensemble Directional Extreme Scenarios (EDES) D q ∩ ∂B a (x) given [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF]. The grey zone is not necessary connected.

In some situations, we might be interested in going beyond just extreme scenarios and include the treatment of the uncertainties on the control parameters and observation data in the optimization solution. Kalman filters provide an elegant framework for such developments. Also, in section 5.3 we go beyond the assumptions above for j over B a (x) using the Ensemble Directional Extreme Scenarios (EDES) where the construction will be applied with the different directions {d l , l = 1, ..., q} built for the ensemble simulation. This propose an ensemble of worst case scenarios as a local sampling of ∂B a (x). The history matching example in section 6 will show a situation where the DES construction is not enough, indeed, and that the ensemble simulation is necessary for a more exhaustive sampling. Still, 2q (q being the ensemble size) the size of this sampling remains negligible compared to n the dimension of the design space. In addition, we will have some confidence on the choice of q through the dimensional analysis of the global sensitivity space D q .

Extended and Ensemble Kalman filters

The Kalman filter (KF) [START_REF] Kalman | A New Approach to Linear Filtering and Prediction Problems[END_REF] provides solution for state estimation of linear systems under Gaussian noise and permits to sequentially assimilate observation data. Extended Kalman filters (EKF) based on the use of the Jacobians of the forecast and observation operators with respect to the state permit to extend the Kalman filter analysis to the case of nonlinear systems [START_REF] Anderson | Optimal Filtering[END_REF][START_REF] Jazwinski | Stochastic Processes and Filtering Theory[END_REF][START_REF] Gelb | Applied Optimal Estimation[END_REF]. To avoid the introduction of these Jacobians, Ensemble Kalman filters (EnKF) [START_REF] Evensen | Advanced Data Assimilation for Strongly Nonlinear Dynamics[END_REF][START_REF] Evensen | Sequential Data Assimilation for Nonlinear Dynamics: The Ensemble Kalman Filter In Ocean Forecasting: Conceptual basis and applications[END_REF] introduce an approximation of the error matrices using statistics generated through scenarios run in parallel. Our aim is to cast our multi-point minimization into an EnKF formulation.

Let us recall the EKF and EnKF algorithms for the sequential solution x k of the following nonlinear dynamical system under uncertainties and assimilating observations y k . [START_REF] Iaccarino | Quantification of Uncertainty in Flow Simulations Using Probabilistic Methods[END_REF] x

k+1 = f (x k , U k ) + q k , y k+1 = g(x k ) + r k ,
where x k , q k ∈ IR n , U k ∈ IR ν and y k , r k ∈ IR π . We usually have n > ν and n > π as we would like to be able to control a system with less control than the size of the system and also there are usually less available observations than this size. Kalman filters are for Gaussian noises and we assume that q k and r k are stationary zero-mean Gaussian perturbations with covariance matrices Q and R. The initial condition x 0 and q k and r k are also assumed uncorrelated.

Defining the Jacobians

A k = ∇ x f (x k , U k ) and C k = ∇ x g(x k
), the EKF algorithms is made of an assimilation (upperscript a) and a forecast (upperscript f) step: [START_REF] Xiu | Numerical Methods for Stochastic Computations: A Spectral Method Approach[END_REF] EKF :

Assimilation:

x a k = x f k + K k (y k -g(x f k )), P a k = (I -K k C k )P f k , Forecast: x f k+1 = f (x a k , U k ), P f k+1 = A k P a k A t k + Q, where K k = P f xy k (P f yy k ) -1 , P f xy k = P f k C t k and P f yy k = C k P f k C t k + R.
We recover the original Kalman filter algorithm if f and g are linear

(f (x k , U k ) = A k x k + B k U k , g(x k ) = C k x k ).
The EnKF algorithm intends to avoid the progression steps for P a k and P f k requiring the Jacobians of f and g. Instead, taking advantage of q parallel forecast for an ensemble of states (x f 1 k , ..., x fq k ), the progression steps above are replaced by the following estimation of the forecast and assimilation states error covariance matrices:

P f k ∼ 1 q -1 E f k (E f k ) t , P a k ∼ 1 q -1 E a k (E a k ) t ,
where E f k and E a k are vectors of the deviation of each of the ensemble member from the mean previsions

x f k = (1/q) q i=1 x f i k and x f a = (1/q) q i=1 x a i k : E f k = (x f 1 k -x f k , ..., x fq k -x f k ), E a k = (x a 1 k -x a k , ..., x aq k -x a k ).
The same procedure is applied to generate from forecast observations (y

f i k = g(x f i k )) the error output vector E y f k = (y f 1 k -y k , ..., y fq k -y k )
which permits to define the following approximation for the forecast observation error covariance matrices:

P f xy k ∼ 1 q -1 E f k (E y f k ) t and P f yy k ∼ 1 q -1 E y f k (E y f k ) t .
The ensemble initialization is based on the probability density functions for x and y:

x f i 0 = x 0 + q i 0 , y i 0 = y 0 + r i 0 , i = 1, ..., q. The assimilation and the forecast steps are then the same for the EKF (and also the KF) and each of the ensemble members in EnKF except that the approximations above are used in the definition of the Kalman gain matrix K k : [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF] EnKF :

     Ensemble Assimilation: x a i k = x f i k + K k (y i k -g(x f i k )),
Ensemble Forecast:

x f i k+1 = f (x a i k , U k ) + q i k , y i k+1 = y k+1 + r i k+1 ,
The ensemble calculations are independent and can therefore be carried out in a fully parallel manner. The algorithm is in particular interesting when it succeeds with an ensemble of size q << n.

5.1. Optimization and Kalman filters. Our aim is to cast the optimization algorithm given in section 2.1 into the Kalman filters formulation. We will be using the ensemble Kalman filter because the systems are nonlinear in our context. We would like to use this approach for robust optimal design accounting for uncertainties in the optimization parameters x and observations y. These gather physical state variables and the constraints which will be seen as observations y. These can be, for instance, global or partial information on x with associated uncertainties.

The choice of the ensemble Kalman filter over the extended Kalman filter is because EKF would require the Hessian of j and the Jacobian of the governing equations and the constraints. These are not always available or easy to get. Our aim is to be able to account for the uncertainties using existing platforms.

The multi-point descent algorithm given in section 2.1 controlling the first and second momentum of a functional can be seen as iterations of:

x k+1 = f (x k , U k ) + q k , y k+1 = g(x k ) + r k ,
where the solution of the governing equations F (and other possible equality constraints) are gathered as observations in y k . U k ∈ IR m represents the sampling of the range of the functioning parameter u and is given by U k = (u 1 , ..., , u m ) ∈ IR m×α (i.e. ν = m × α). Recall that we have u ∈ I m ⊂ I ⊂ IR α with m the size of the sampling. q k and r k represent the uncertainties on x k and y k . These are zero-mean Gaussian with covariance matrices Q and R. In the sequel the approach is therefore limited to the Gaussian framework.

Illustration on quadratic minimization under equality constraints.

Let us start with a simple optimization problem under equality constraints and without any governing equation involved. We consider the minimization of a quadratic functional under two equality constraints:

j 0 (x, u) = 1 2 n i=1 (x i -u i n ) 2 , c 1 (x) = n/2 i=1 x i - 1 2 = 0, c 2 (x) = n i=n/2+1 x i -1 = 0.
The size of the optimization problem is n = 100. Let us start with u = 1 being a constant and not given over an interval I ⊂ IR. Without constraints the optimal solution is obviously x i = i/100. With constraints and in the absence of noise this can be solved considering a functional of the form:

(13) j(x, u) = j 0 (x, u) + 1 2 c 2 1 (x) + 1 2 c 2 2 (x).
First order optimality condition for j leads to a n × n linear system Ax = b(u) where A i,i = 2, A i =j = 1 when i&j n/2 or i&j n/2 and A i =j = 0 elsewhere and b i = u i/n -1/2 for i n/2 and b i (u) = u i/n -1 for i > n/2. The solution of this system is shown in figure 2.

In the presence of uncertainties on the optimization parameter x and the constraints c i , this can be cast into the Kalman filter formalism (classical here as we are in a linear dependency situation) with c i seen as two observations because these constraints can be considered as a priori global information on x k , and only partially satisfied (noise is the error one commits on the equality constraints):

(14) f (x k , u k ) = x k -ρ∇ x j(x k , u k ) y 1 k = c 1 (x k ), y 2 k = c 2 (x k ), where ∇ x j(x k , u k ) = Ax k -b(u k ) with U k = u k = u
and ρ = 0.01. The Gaussian uncertainties on the optimization parameter x k and constraints y 1,2 k are given through their correlation matrices Q = 0.02I 100 and R = 0.05I 2 . Figure 2 shows comparisons between the unperturbed solution and when uncertainties on the optimization parameters and the constraints realizations have been reduced by the EnKF algorithm. Ensembles of size 5 to 20 have been tested and the convergence histories for the deviations between the assimilation and forecast ensemble solutions and the deviation between ensemble solutions at given time k-index.

5.2.1.

Parametric quadratic minimization under equality constraints. Now consider the situation where u is no longer a constant but belongs to an interval I = [0.8, 1.8]. As we discussed in section 2, this is quite representative of what we often have in practical situation with u being a functioning parameter and we would like the optimal solution to be robust with respect to u: the functional to be as much as possible flat on its domain of variation, for instance. We consider I m a m = 20 points uniform sampling of I and J(x) = µ defined in (3) and with j(x, u l ) given by [START_REF] Karniadakis | Multi-element generalized polynomial chaos for arbitrary probability measures[END_REF]. With the notations above, the first order optimality condition leads this time to Ax = u∈Im b(u). As for the constant u case, this can be cast into the Kalman filter formalism in the presence of uncertainties on the optimization parameter x and the constraints c i with ∇ x J(x k , u k ) = (Ax k -u∈Im b(u))/m and the size of the control U k = (u k 1 , ..., u k m ) being now m. The multi-point descent direction d(µ, σ) can be built as in section 2.1 and the multi-point descent algorithm can be cast into the Kalman filter framework replacing ∇ x j(x k , u k ) by d in [START_REF] Congedoa | Shape Optimization of an Airfoil in a BZT Flow with Multiple-source Uncertainties[END_REF]. Figure 2 shows comparisons between the mono-point and this parametric minimization for ensembles of size 5 to 20. The parametric solution proposed by the multi-point minimization deviates quite from the deterministic solution. Even if this might appear being counter-intuitive, this means that it should be more difficult to fit the data in a parametric inverse problem. We will recover this behavior in section 6 for a history matching problem. 5.3. Geometry of the global EnKF search space and ensemble size. At this point one can take advantage of the descent directions ∇ x j(x l k , u) or d(µ(x l k ), σ(x l k )) built during our ensemble simulations and proceed as in section 3 to analyze the geometric characteristics of the related global search space D q in the context of single point optimization (for a single u value): D q = Span{∇ x j(x l k , u), l = 1, ..., q} ⊂ IR n×q , or multi-point (when u ∈ I) optimization:

D q = Span{d(µ(x l k ), σ(x l k
)), l = 1, ..., q} ⊂ IR n×q . An important interest of these spaces is to provide information on what should be the size of the ensemble in EnKF. This question is treated empirically in practice increasing the size q until the EnKF iterations succeed. An ensemble can be considered as enough rich if dim(D q ) < q. In addition to their respective dimensions, principal angles between two spaces D q and D q show how much deviation appears between subspaces for different ensemble sizes. This permits to define a new concept of convergence for the ensemble simulation.

This geometric construction can be used to define possible extreme scenarios introducing quantiles such as the value at risk (VaR) as presented in section 4. More precisely, ensemble directions d l k = d(µ(x l k ), σ(x l k )), l = 1, ..., q permits to define a sampling of what one might call an extreme scenarios region around a mean inversion solution x k : [START_REF] Cinnella | Robust optimization of dense gas flows under uncertain operating conditions[END_REF] D q ∩ ∂B a (x k ), where x k = (1/q) q l=1 x l k . Hence, applying the Directional Extreme Scenarios (DES) concept to the each of the ensemble directions provides a sampling of the extreme scenarios region through the following Ensemble Directional Extreme Scenarios: [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF] (EDES) : {(x l k ) ± , l = 1, .., q} = {d l k ∩ ∂B a (x k ), l = 1, ..., q}. A sketch of this construction is given in figure 1 and compared to the situation presented in section 4. These extreme scenarios will be illustrated in the examples of section 6 and compared to the results by simple DES when the Ensemble Kalman filter is not used and the uncertainty of x only accounted for through two directional VaR-based extreme scenarios (in mono and multi-point situations for u) permitting for low-complexity analysis of the impact of variability in x on the inversion.

Application to a history matching problem

We consider an inverse problem appearing in several science and engineering domains such as electroencephalography (EEG) inverse source localization and history matching in reservoir modeling [START_REF] Caers | Modeling Uncertainty in the Earth Sciences[END_REF]. We discuss this latter problem.

Reservoir simulation is a major activity for many oil companies. It starts with building the initial geological reservoir and fluid models and their discrete versions. These can be either PDE-based or relying on reduced order models, for instance involving response surfaces. Our reservoir model is PDE based but other solutions can be envisaged and this choice is not central to our discussion. In all cases, and as for many models in science and engineering, one should be aware that reservoir models are not reality and have inevitable errors also called reducible or epistemic uncertainties. After the reservoir model is established, it must be improved through history matching which must include a quantification of uncertainties. This step consists of adjusting the model of a reservoir until it closely reproduces its past behavior available through the historical production and pressures data which are both noisy. These are not reducible and are called aleatory uncertainties.

In addition to these data, observations and constraints (which can be considered as a priori information and somehow as observations) are also available. Examples of such are well logs and information there on the core, seismic data, geological and production scenarios, etc.

Previous works exist on the application of ensemble Kalman filter methods to history matching problems (see [START_REF] Seiler | Using the EnKF for history matching and uncertainty quantification of complex reservoir models , Computational Methods for Large-Scale Inverse Problems and Quantification of Uncertainty[END_REF][START_REF] Liu | Ensemble Kalman filter for automatic history matching of geologic facies[END_REF][START_REF] Heidaria | History matching of petroleum reservoir models by the Ensemble Kalman Filter and parameterization methods[END_REF] and references therein). EnKF being in the class of Monte Carlo simulation methods, it is usually stated that EnKF is suitable as it avoids gradient based minimization algorithms which can be trapped by local minima. We think that one should take advantage of the gradients information when available and rather than dropping gradient-based minimization algorithms cast them into the EnKF framework. Gradient information also reduces the risk for a pure EnKF estimates to deviate strongly from petrophysical properties. In addition, thanks to the knowledge of the dynamic global sensitivity spaces built during EnKF iterations, one have a good estimation of what the ensemble size should be.

To illustrate our purposes we consider a synthetic case with two injection and two productions wells in a 100m × 100m domain. We consider a simple time dependent Darcy-like model [START_REF] Darcy | Les Fontaines Publiques de la Ville de Dijon[END_REF] over a domain Ω linking a single phase pressure p to the flow velocity v (gravity is neglected):

(17) ∂ t p + ∇.v = i=1,2 u i δ s i , v = -x∇p,
which is solved on a 100 × 100 cartesian mesh following [START_REF] Chavent | Discontinuous upwinding and mixed finite element for twophase flows in reservoir simulation[END_REF] and using a second order time integration. The literature on the solution of this equation is very rich.

A particular choice of a numerical scheme is not central to our discussion. The algorithm given in section 2.1 requires the adjoint of the governing equation in order to make the cost of sensitivity evaluation independent of the size n of the optimization space for x. To make the code development as automatic as possible, the adjoint solver has been generated using automatic differentiation by Tapenade [START_REF] Hascoet | Tapenade 2.1 user's guide[END_REF][START_REF] Griewank | Computational derivatives[END_REF].

In equation [START_REF] Ide | Unified notation for data assimilation: operational, sequential and variational[END_REF] x describes the soil characteristics and our aim is its identification. More precisely, x = k(µφ) -1 with k(m 2 ) the permeability, µ(bar . hour) the fluid viscosity and 0 < φ 1 is the porosity of the soil. An a priori available information on x is that it must be piecewise constant over a number of plateau ω i . This is accounted for through a projection step on x. We have here an example of information provided by geophysicists from seismic observations. Hence, the optimization parameter K over the computation mesh is projected over this representation defining a set of couples {(K i , ω i ), i = 1, ..., n p } where

K i = 1 |ω i | ω i x(s)ds, ∪ n i=1 ω i = Ω.
In the examples presented here, ω i describes transverse strips as shown in figure 3 and we have a number of plateau n p n with n = 100 the size of the optimization space. This a priori information is very useful for an efficient ensemble simulation permitting to keep the ensemble size low as shown in figure 5.

In addition to the information above we have well logs and x is known at wells and this information defines the observation equations: y i=1,...,4 k = x k (s i ) -x * (s i ) = 0, at s i the locations of the 4 wells.

The control variables u ∈ IR 2 + represents the injection pressure at the two injection wells s 1,2 . This is an example where u ∈ IR 2 + in j(x, u) of section 2. We aim at recovering x minimizing:

j(x, u) = i=3,4 T 0 (p i (t, x, u) -p * i (t)) 2 dt,
from histories p * i=3,4 (t) recorded at production wells s 3,4 . p * i (t) is uncertain but this uncertainty is assumed small in comparison to what expected on the permeability and the observations. The physical time (100 hours) over which these histories have been recorded and used for history matching is different from the k-index time in the EnKF iterations. These correspond to optimization iterations.

The uncertainties on the optimization parameters x k and constraints y 1,...,4 k are both assumed Gaussian with correlation matrices Q = 0.02I 100 and R = 0.05I 4 .

Our aim is to identify x by our gradient-based+EnKF inversion approach accounting for the different uncertainties in a mono-point situation and when the control u is defined through its range of variation. . Target x after projection and definition of the set {(K i , ω i ), i = 1, ..., 10} where K(= k(µφ) -1 ), field pressure at t = 100h and pressure histories p * i=3,4 (t) recorded at production wells s 3,4 over 100 hours and used as target for inversion.

6.1. Mono-point inversion under uncertainties. Before considering the geometric analysis of the global sensitivity spaces in a multi-point situation, let us show an example of the results one gets with the EnKF-based inversion in a mono-point configuration defined with constant values over time for the two components of the control u at 1 bar and 1.3 bar at injection wells. These are not realistic values and are just illustrative. The governing equation being linear in p, there is no difficulty in considering larger values for the injection pressure and the inversion will succeed as far as the target pressure histories are compatible.

To be sure the problem is admissible, an unperturbed target solution has been generated with a piecewise constant permeability distribution with n p = 10 plateau (shown in figure 3).

We consider three ensembles of size q = 10, 20 and 40. Figure 4 shows the difference between the target x and p and the mean ensemble estimates for the ensemble of size 40. It also shows the deviation between ensemble solutions. These results permit to see the regions of maximum deviation between the EnKF estimates and also between these EnKF estimates and the solution by a deterministic inversion (i.e. when the inversion does not account for the uncertainties on the data and observations). It is also interesting to see that the target pressure histories at production wells s 3,4 are well recovered by the ensemble mean solution. But that the uncertainty on the production at the two wells are different.

Figure 5 shows the histories over time k-index of the deviation between the ensemble analysis and forecast EnKF estimates and the root-mean-square deviation of the forecast estimates for three ensembles of size q = 10, 20 and 40. This problem is more complex than the quadratic minimization of section 5.2. We have the same level of uncertainties in both problems and the same size of the problem (n = 100). Still an ensemble of size 20 seems to produce convergent statistics. This can be explained with our geometric analysis. Figure 5 also shows the outcome of the Gram-Schmidt algorithm detecting independent directions in the global sensitivity space D q = Span{∇ x j(x l k , u), l = 1, ..., q} for q = 20 and 40 as defined in section 5.3 and eventually their dimensions which remain below 10. This gives an indication on how to choose the ensemble size in EnKF.

The figure also illustrates the application of the concept of principal angles between subspaces to estimate the relative position of the global sensitivity spaces D 20 and D 40 . In addition to the dimension of the spaces, these angles provide additional information on the impacts on the search spaces due to an enrichment of the ensemble.

6.1.1. Var-based Directional Extreme Scenarios in the mono-point case. As described in section 5.3, the EnKF inversion permits to define a sampling of the extreme scenarios region [START_REF] Cinnella | Robust optimization of dense gas flows under uncertain operating conditions[END_REF]. We also saw in section 4 how to define two VaR-based directional extreme scenarios during a gradient based minimization. We would like to compare the outcome of the two approaches with the EnKF analysis considered as the reference solution.

Knowing the Gaussian probability density function for x k , we define two VaRbased scenarios along ∇ x j(x k , u):

(x k ) ± = ∇ x j(x k , u) ∩ ∂B 0.95 (x k ), where x k = x k-1 -ρ∇ x j(x k-1 , u) with ∂B 0.
95 is the boundary of the VaR 0.95 -based variability ball defined by (8) around x k . Knowing the correlation matrix Q of Gaussian variables x k , the extreme scenarios have their i th component defined by:

(18) (x k ) ± | i = x k | i (1 ± Q ∇ x j(x k , u) ∇ x j | i ).
In our case Q is diagonal (Q = 0.02I 100 ) and from section 4 we have VaR 0.95 (N (0, 0.02)) = 0.02 VaR 0.95 (N (0, 1)) = 0.02 × 1.65 = 0.033, which defines the value of the i th component of x ± k following the sign of the i th component of ∇ x j(x k , u):

(x k ) ± | i = x k | i (1 ± 0.033 sgn(∇ x j) i ).
In the same way, 2q (x l k ) ± VaR-based Ensemble Directional Extreme Scenarios can be defined applying [START_REF] Casella | Statistical Inference[END_REF] to each of the ensemble members x l k , l = 1, ..., q. Figure 6 shows pressure histories for these ensemble extreme scenarios after an inversion with an ensemble of size q = 40 giving 80 scenarios {(x l 30 ) ± , l = 1, ..., 40} after 30 iterations of EnKF. As reported in figure 4, the production at well 4 is more sensitive to the variability in x than at well 3. The figure also compares these with the pressure histories for the two VaR-based scenarios (x 30 ) ± by [START_REF] Casella | Statistical Inference[END_REF]. This lowcomplexity estimation of the impact of the variability in x gives satisfactory bounds in this single-point situation. The situation is less satisfactory in the multi-point situation justifying ensemble simulations. ) ± , l = 1, ..., 40 described in section 6.1.1. Production at well 4 is more sensitive to the variability in x. Lower: pressure histories for the two VaR-based scenarios given by ( 18) with a single-point gradient based minimization (no ensemble simulation involved).

[1.1, 1.5]. We would like to have the inversion to be robust also with respect to this variability. The multi-point minimization algorithm given in section 2.1 is cast into the Kalman filters framework as described in section 5.1 to control the first two moments µ and σ of j defined in section 2 accounting for the variability on both x and u.

We consider again three ensembles of size q = 10, 20 and 40. The multi-point optimization of 2.1 is defined for I m being a 10 × 10 sampling of I.

Figure 7 shows the histories over time k-index of the deviation between the ensemble analysis and forecast EnKF estimates and the root-mean-square deviation of the forecast estimates for three ensembles of size 10, 20 and 40. Figure 7 also shows the outcome of the Gram-Schmidt algorithm detecting independent directions in the global sensitivity space D q = Span{d(µ(x l k ), σ(x l k )), l = 1, ..., q} defined in section 5.3 and eventually its dimension which remains below 15. This is often larger than what we had in the mono-point situation in figure 5, but sometimes smaller (at given iteration) and would then mean that in these situations the multi-point direction constructed to account for the variability on u actually reduces the effect of the variability on x. In all cases, this explains why an ensemble of size 10 is clearly not enough rich. The figure also illustrates the application of the concept of principal angles between subspaces to estimate the relative orientation of the global sensitivity spaces D 20 and D 40 for the ensemble sizes of 20 and 40. These spaces are generated by the 20 (resp. 40) directions d(µ(x l k ), σ(x l k )) as described in section 5.3. In addition to the dimension of the spaces, these angles provide additional information on possible impacts on the search spaces due to an enrichment of the ensemble. In the present case, the angles increase with the EnKF iterations which mean that one needs more accuracy with convergence and that the additional information by the ensemble of 40 is not present in what available with an ensemble simulation of size 20.

Figure 8 shows the difference between the target x and p and the mean ensemble estimates for the ensemble of size 40, to be compared with figure 4 in the monopoint situation for u. The inversion gives a very different estimation for x when the control u is not anymore single valued. Also field deviations between the EnKF estimates and between these EnKF estimates and the solution by a deterministic inversion are different from what obtained in the mono-point situation showing the utility of the multi-point formulation. The pressure histories at production wells s 3,4 also deviate more from the targets which were recovered by both the inversions in the deterministic and uncertain mono-point cases. The deviation is not the same for the two wells. One should have in mind that there is no guarantee an inversion finding a distribution of x maximizing the likelihood with respect to the uncertain parameters and observation data would produce a design giving the deterministic target pressure output. Indeed and surprisingly, despite one has lower deviation between ensemble pressure simulations for well 3 (Figure 8, lower-right), the inversion under uncertainty produces a solution farther from the deterministic target than for well 4. This is a kind of information which help discriminating between different production points in term of reliability and necessary monitoring.

As in the mono-point situation, these can be compared to the two scenarios obtained by the Directional Extreme Scenarios (DES) concept presented in section 5.3 which permits a low-complexity quantification of the uncertainties due to the variability on x.

6.2.1.

Var-based Directional Extreme Scenarios in the multi-point case. We proceed with the same analysis than in section 6.1.1 in the single-point situation but using the directions d(µ(x k ), σ(x k )) instead of ∇ x j(x k , u) and the multi-point control algorithm of section 2.1 instead of a simple descent for j as presented earlier.

Figure 9 (to be compared to figure 6 in the mono-point case) shows pressure histories for the ensemble extreme scenarios after an inversion with the multi-point optimization algorithm of 2.1 and an ensemble of size q = 40 after 30 iterations of EnKF. The figure also compares these with the pressure histories for the two directional VaR-based scenarios (x 30 ) ± defined by:

(19) (x 30 ) ± = x 30 (1 ± Q d(µ(x 30 ), σ(x 30 )) d(µ, σ) ).
This low-complexity estimation does not give satisfactory bounds for the extreme scenarios in this multi-point situation. This is unfortunate as otherwise there would be no need for an ensemble analysis to access the extreme scenarios. This situation shows that we need to be cautious, especially when uncertainties are present in both

x and u. In these cases, the variability in the search direction d(µ(x l ), σ(x l )) is larger than what we have in the mono-point situation with ∇ x j(x l , u) for l = 1, ..., q. This can also be seen from increasing principal angles comparing figures 5 and 7. Indeed, when accounting for the variability in u (i.e. going from mono to multi-point) the angles between subspaces D 20 and D 40 increase as an indication of larger variability. As a consequence, when accounting for the uncertainty on u we need to a posteriori define the worst case scenarios from a more exhaustive ensemble analysis.

CONCLUDING REMARKS

In order to be easily integrated in engineering environments and quantify our confidence on the optimal solution without intensive sampling of large dimensional parameter spaces a cascade of geometric uncertainty quantification concepts, with growing computational complexity, has been introduced.

The analysis starts with the geometric characterization of global sensitivity spaces through their dimensions and relative positions by the principal angles between global search subspaces. Then, joining a multi-point descent direction and extreme values information from the probability density functions of design variables the concept of Directional Extreme Scenarios (DES) has been introduced.

The construction goes beyond DES with Ensemble Kalman Filters (EnKF) after the multi-point optimization algorithm is cast into an ensemble simulation environment. This permits to account for the variability on the functioning parameters through the multi-point formulation and for the variability on the optimization parameters and observation data through the ensemble Kalman filter formulation.

The UQ cascade ends with the joint application of the EnKF and DES leading to the concept of Ensemble Directional Extreme Scenarios (EDES) which provides exhaustive possible extreme scenarios knowing the Probability Density Function of the optimization parameters and this without a sampling of the admissible space.

These ingredients have been illustrated on an history matching application. The interest of casting a multi-point minimization algorithm into an EnKF framework has been demonstrated when simple directional extreme scenarios do not provide satisfactory bounds for the impact of the uncertainties on the optimization parameters on the results.

The method (both the DES and EDES constructions) can theoretically be applied to any level of uncertainty on x. But this will increase the regularity constraints on the functional over B a (x) as discussed in section 3 and 5.3. Also, the Kalman ensemble size will need to be much larger as the dimension of the ensemble sensitivity space D q will certainly increase due to this variability.

The Gaussian framework is necessary for the Kalman filter framework. But, the other ingredients of the paper can be applied to non-Gaussian situations: dimensional analysis, principal angles, quantile-based extreme scenarios as far as one can (even by numerical integration of tabulated laws) find the a-quantiles quantifying the tail of the PDFs.

The framework presented features several sources of natural parallelism making the time to solution comparable to a mono-point deterministic inversion.

Looking at the bigger picture, this analysis shows the importance of considering all sources of variability in the problem to make sure the global search space is neither too large nor too small. The latter should obviously be avoided for the design to be robust. But, the former tells us that a more efficient design is maybe possible better monitoring the design space.

  -Square deviation of EnKF estimates

Figure 2 .

 2 Figure 2. Mono-point (left column) and parametric (right column) quadratic minimization under equality constraints: comparison between the unperturbed solution and when uncertainties on the optimization parameters and the constraints realizations have been reduced by the EnKF algorithm.

Figure 3

 3 Figure 3. Target x after projection and definition of the set {(K i , ω i ), i = 1, ..., 10} where K(= k(µφ) -1 ), field pressure at t = 100h and pressure histories p * i=3,4 (t) recorded at production wells s 3,4 over 100 hours and used as target for inversion.

Figure 4 .

 4 Figure 4. Inversion with EnKF (ensemble of size 40) in a monopoint u case. Differences between the ensemble mean solutions from the target x and the field pressure p(t = 100h) and the deviation between ensemble solutions. Target pressure histories at production wells s 3,4 are well recovered by the ensemble mean but uncertainties on the production histories do not follow the same pattern.

Figure 5 .Figure 6 .

 56 Figure 5. Inversion with EnKF (ensemble of size 10, 20 and 40) in a mono-point u case. Deviation between EnKF analysis and forecast estimates during EnKF-gradient iterations and the RMS of EnKF x f k estimates. The histories of the Gram-Schmidt orthonormalization during ensemble iterations (middle-left) and the dimension of the free subspace generator for spaces D 20 and D 40 (middle-right) explains why an ensemble of size 10 seems to be efficient. Lower: principal angles between subspaces D 20 and D 40 defined in section 5.3.

Figure 7 .

 7 Figure 7. Inversion with EnKF (ensemble of size 10, 20 and 40) in a multi-point u configuration. Deviation between EnKF analysis and forecast estimates during EnKF-gradient iterations and the RMS of EnKF x f k estimates. The histories of the Gram-Schmidt orthonormalization during ensemble iterations for the ensemble of size 40 (middleleft) and the dimension of the free subspace generator for space D q (middle-right) to be compared to figure 5. Lower: principal angles between subspaces D 20 and D 40 defined in section 5.3.

Figure 8 .

 8 Figure 8. Inversion with EnKF (ensemble of size 40) in a multipoint u configuration. Differences between the ensemble mean solutions from the target x and the field pressure p(t = 100h) and the deviation between ensemble solutions. Target pressure histories at production wells s 3,4 are not well recovered by the ensemble mean and the uncertainties on the production histories do not follow the same pattern than in the mono-point case. The results are quite different from those with a mono-point inversion in figure4.

Acknowledgements

The adjoint code for the history matching problem has been obtained using Tapenade AD tool developed at INRIA-Sophia Antipolis by L. Hascoet and his team Tropics.

) ± , l = 1, ..., 40 described in section 6.1.1. Lower: pressure histories for the two directional VaR-based scenarios given by [START_REF] Hoel | Introduction to Mathematical Statistics[END_REF] with the muti-point minimization algorithm 2.1 (no ensemble simulation involved) inefficient in this multi-point case.