# Low-complexity moment-based robust design and uncertainty back-propagation

Radia BOUABDALLAH¹, Bijan MOHAMMADI¹, Robert RAPADAMNABA¹

¹IMAG, Univ. Montpellier, CNRS, Montpellier, France

*E-mail: radia.bouabdallah@gmail.com

DOI: 10.46298/arima.7160

Submitted on February 8, 2021 - Published on September 3, 2021

Volume: 32 - Year: 2021

Editors: Eric Badouel, Nabil Gmati, Maurice Tchuenté, Bruce Watson

**Abstract**

The paper shows how to take advantage of a possible existing linear relationship in an optimization problem to address the issue of robust design and backward uncertainty propagation lowering as much as possible the computational effort.

**Keywords**

uncertainty quantification, robust optimization, backward propagation, low-complexity.

# I INTRODUCTION

The performance of a system designed for given functioning conditions often seriously degrades when these conditions are modified. Today's industrial robust design mainly relies on reduced-order modeling and intelligent sampling [10, 13, 16, 20] which either does not use high-fidelity simulations during design or uses lower accuracy than what would be affordable in a single-point optimization.

By robust design we mean a proposition which ensures similar performances over a given operation range. We will discuss the implication of this requirement on the moments of the performance functional. Our aim is to propose a methodology which permits to design a system having similar performances over a given range of its operating condition or functioning parameters. From a practical point of view, we would like this to be achieved modifying as less as possible an existing single-point optimization platform. In particular, we would like to avoid replacing the high-fidelity ingredients of the platform by low-complexity solvers. Finally, we would like the parallel time-to-solution to remain comparable to the single-point situation. We will show that our direct optimization algorithm permits to achieve such a task.

We address the solution of robust moment-based optimization problems after a multi-point reformation. The first four moments are considered (i.e. mean, variance, skewness and kurtosis).
going beyond classical engineering optimization based on the control of the mean and variance. In particular, the impact on the design of a control of the third and fourth moments is discussed. The multi-point formulation leads to discrete expressions for the moments. linking moment-based and multi-point optimizations. The linearity of the sums in the discrete moments permits an easy evaluation of their gradients with respect to the design variables. Optimal sampling issues are analyzed and a procedure is proposed to quantify the confidence level on the robustness of the design. The proposed formulation is fully parallel and the parallel time-to-solution is comparable to single-point situations.

The literature on multi-point optimization is vast and exhaustive referencing is out of our scope. This formulation has been used, in particular, to extend single-point optimization to account for additional operating conditions.

Forward and backward uncertainty propagation are obviously of great importance with a huge literature dedicated to uncertainty quantification (UQ) [7, 8, 21]. Forward propagation aims at defining, for instance, a probability density function (PDF) for a functional $J$ knowing those of optimization variables [2, 3, 11]. 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) [5, 14, 16, 19]. Backward propagation aims at reducing models bias or calibrating models parameters knowing the PDF of $J$ [4, 6, 9]. This can be seen as a minimization problem and Kalman filters [1] give, for instance, an elegant framework for this inversion assimilating the uncertainties on the observations.

We particularly discuss one aspect of UQ where the target state $u^*$ used in an inverse problem is uncertain. This is the case, for instance, minimizing $J(u(\mathbf{x}), u^*) = \|u(\mathbf{x}) - u^*\|$ to reduce the distance between a model state $u(\mathbf{x})$ and observations or target $u^*$.

Targeting uncertain data is a realistic situation as the acquired data are usually uncertain. It is therefore interesting to be able to quantify the impact of this uncertainty on the inversion results. An important information will be the sensitivity of the design to a given level of uncertainty on the data at some location. Indeed, if this sensitivity is low, this would be an indication that a more accurate acquisition there is unnecessary.

Considering the target as uncertain is also interesting because we do not always have existence of a solution for an inversion problem as $u^*$ is not necessary solution of the state equation making an exact or deterministic inversion pointless. Also, the approach permits to go beyond inversions based on least-squares minimization involving a mean state target.

The uncertainty in measurements is also an interesting way to account for the presence of variability in the state. More generally, as a model or numerical procedure is by nature imperfect and partial, we can consider this uncertainty as a representation and estimation of these imperfections. These imperfections are even more present in inverse problems where one cannot afford the same level of resolution than for a single simulation. We therefore need to be able to quantify the impact of these weaknesses. Being able to handle uncertain targets in inverse design is therefore also useful to account for epistemic uncertainties related to possible model or solution procedures deficiencies.

The paper also addresses the issue of backward uncertainty propagation through the estimation of the covariance matrix of the optimization variables. Two procedures for the estimation of this matrix are presented taking advantage of the existence of local linear relations between the target state $u^*$ and the optimization variables $\mathbf{x}$.
We present the different ingredients of the paper in the context of a linear problem for load distribution optimization with the objective of achieving a target deformation or displacement of a shape. The state equations we consider are those of linear elasticity and we would like to take advantage, as much as possible, of the linear relationship between the optimization variables (e.g loads) and the state variables (e.g the displacements) in the solution of the problem and the different uncertainty quantifications.

The paper starts with the description of our single-point optimization problem to illustrate the situation where the solution of the inverse problem can be explicitly expressed thanks to the mentioned linearity. Then we address robust parametric optimization and its solution through moment-based and multi-point formulations. The paper ends with two low-complexity approaches for backward uncertainty propagation and estimation of the covariance matrix of the optimization variables.

## II STATE EQUATIONS AND SINGLE-POINT OPTIMIZATION PROBLEM

Let us start with the single-point optimization problem of interest. We are interested by inverse problems where the functional typically measures a distance to a target solution.

We consider a generic state equation with a linear state equation and also linear dependency between the optimization variables x (controlling in our case the loads distribution) and the state variable u: $\mathcal{F}(u(\mathbf{x})) = f(\mathbf{x})$, where u is therefore a function of x.

Our single-point optimization problem reads:

$$ \min_{\mathbf{x} \in O_{ad}} J(u(\mathbf{x}), \alpha) \quad \text{under the constraint that} \quad \mathcal{F}(u) = f(\mathbf{x}). \qquad (1) $$

Where $\alpha$ represents the other independent variables of the problem such as geometry entities or operating conditions [4, 23]. Here, for simplicity we do not mention all dependencies in u, such as the spatial dependency and so forth. In what follows, $\alpha$ is the target solution $u^*$ which we consider as uncertain and the functional measures the least-square misfit between the state and this target.

In single-point optimization $\alpha$ is fixed in opposition to multi-point or robust optimizations where $\alpha$ components will no longer be fixed but given either through intervals or more sophisticated probability density functions.

The source term $f(\mathbf{x})$ contains $N$ loads expressed as $x_i\delta(s_i)$ where $\delta(s_i)$ is a spatial localization function in space variable $s$. One would like to optimize:

$$ f(\mathbf{x}) = \sum_{i=1}^{N} x_i \delta(s_i). \qquad (2) $$

The optimization variable is then $\mathbf{x} = (x_1, ..., x_N) \in (\mathbb{R}^+)^N$.

The inverse problem of interest aims at making the solution $u$ as close as a target $\alpha$ over the calculation domain $\Omega$, or some subset of it. Using a least-squares formulation over the domain with spatial integration, this reads:

$$ J(u(\mathbf{x}), \alpha = u^*) = \frac{1}{2} \int_{\Omega} (u(\mathbf{x}) - u^*)^2 ds. \qquad (3) $$
Taking advantage of the linear dependency between x and the solution u and using the super-position principle, the solution can be expressed through projections over $N + 1$ elementary solutions:

$$u = u_0 + \sum_{i=1}^{N} x_i u_i,$$

In this decomposition, $u_0$ is the homogeneous solution (the zero load situation together with non-homogeneous boundary conditions) and the other elementary solutions $u_{i=1,...,N}$ are obtained respectively with $x = (0,..., 0, 1, 0,..., 0)$ with 1 at the $i^{th}$ position. The elementary solutions are obtained with homogeneous Dirichlet boundary conditions.

Once the homogeneous and elementary solutions known, the optimization variable x is solution of a linear system:

$$Ax = b(\alpha = u^*), \quad (4)$$

of size $N$ derived from first order optimality condition for $J$, $\nabla_x J(u(x), \alpha) = 0$, with:

$$A_{i,j} = \langle u_i, u_j \rangle \quad \text{and} \quad b_i(\alpha = u^*) = \langle (u_0 - u^*), u_i \rangle, \quad \text{for } i, j = 1, \dots, N,$$

$$\text{where } \langle v, w \rangle = \int_{\Omega} v(s) w(s) ds.$$

A fundamental remark now is that the cost of this inversion is negligible compared to the cost of calculating the elementary solutions. Also, it is important to notice that the matrix $A$ does not depend on $\alpha$. Therefore, any numerical transformation necessary for its inversion is made only once. Uncertainty quantification is situation where the target solution $\alpha$ is uncertain becomes therefore feasible at low cost as it will not require any extra solution of the state equations. In the sequel we will discuss different ways to handle this problem.

## 2.1 3D elasticity

We will illustrate the ingredients of the paper on a 3D problem where the state equation $F(u) = f(x)$ is the elasticity system whose solution is given using the Cast3M [22] industrial solver on tetrahedral meshes which we consider as a black-box. Something which is made possible thanks to the formulation of the optimization problem and this is one strength of the approach.

Briefly speaking, Cast3M is a finite element software for structural and fluid mechanics developed by the French Alternative Energies and Atomic Energy Commission (CEA). Cast3M provides a mesh generation tool, partial differential equations simulation modules coupled with visualization tools.

The 3D problem we consider consists of a geometry of an elastic piece represented in a $(10m \times 10m \times 1m)$ concrete slab with hollow bodies [12], supported by point supports at the four corners modeled by homogeneous Dirichlet boundary conditions for the displacement. The mechanical characteristics of the concrete considered are $E = 30GPa$ for Young's modulus and $\nu = 0.2$ for Poisson's ratio. This is a complex problem and a full 3D calculation is quite time consuming. The present solution based on the use of elementary 3D solution is therefore quite adapted for the optimization of load distributions for such a cases as the state equations calculations can be carried out independently from the optimization problem and in a fully parallel way. The hollow bodies involve complex mesh generation issues as shown in Figure 1. The 3D mesh has about 14400 grid points.
<figure><img src="image_1.png" /><figcaption>Figure 1: Example of a complex mesh generated by Gmsh [18] taking into account the presence of hollow bodies.</figcaption></figure>

We assume negligible volume forces and consider a situation with a relatively small number of loads to ease the comparison between the different approaches. Considering higher number of loads does not pose any supplementary difficulty and only requires parallel calculations of additional independent elementary solutions.

## 2.2 Illustration using $N = 5$ optimization variables

We apply an arbitrary five point loads distribution (i.e. $N = 5$). Any other combination can be considered, five is suitable as enough large to illustrate the approach and enough low to show images of all elementary contributions.

We calculate the elementary displacements $u_{i=1,...,5}$ obtained after application of unit forces $x_i$.

The total displacement is the sum of these elementary displacements $u = u_0 + \sum_{i=1}^{N=5} x_i u_i$ where $u_0$ is the displacement without any forces applied but in the presence of Dirichlet boundary conditions. In our problem we consider $u_0 = 0$.

Figure 2 shows the elementary and total structural displacements which will be considered as target deterministic displacement $u^*$. Uncertainty analysis will take place around this target displacement.
<figure><img src="image_2.png" /><figcaption>Figure 2: The 5 first pictures represent elementary displacements $u_{i=1,...,5}$ and the bottom right picture shows the total structural displacement $u = u_0 + \sum_{i=1}^{N=5} x_i u_i$, where $u_0 = 0$ in this case. The elementary displacements cannot be retrieved using symmetry considerations from each other.</figcaption></figure>

### III ROBUST PARAMETRIC OPTIMIZATION

Consider a cost function involving the optimization variables $x$ in an admissible ensemble $O_{ad}$ and another parameter $\alpha$, not considered as a design variable:

$$ \min_{x \in O_{ad}} J(x, \alpha), \quad \alpha \in I \subset \mathbb{R}^n, \ O_{ad} \subset \mathbb{R}^N. \qquad (5) $$

As mentioned, $\alpha$ can be an operating condition or the target solution in an inverse problem. Both the operating conditions and the target solution can be uncertain and given only through
probability density functions. Handling the uncertainty on the operating conditions is necessary to achieve robust design and usually requires a multi-point formulation. The number of operating condition parameters is usually small compared to the size of the optimization problem. This is a very general context and we previously visited it to address robustness issues in optimization with respect to x and α in general frameworks [25–27, 30].

## 3.1 Moment-based optimization

Multi-point optimization can be used to address such optimization problem [24]. The aim is to remove the dependency in α during optimization. This can be done, for instance, minimizing a functional $J(x)$ encapsulating this dependency under a constraint on the higher moments of $J$:

$$J = \mu = \mathbb{E}[J(x, \alpha)], \quad \text{such that} \quad \sigma = \mathbb{E}[(J(x, \alpha) - \mu)^2] \le \sigma_0, \qquad (6)$$

where $\mathbb{E}[\cdot]$ is the expectation operator.

We can go beyond the first two moments and use also the third and the fourth moments. Going beyond the first two moments is important when the PDF of $J$ deviates from a pure Gaussian distribution. Indeed, even with interval-based (with uniform PDF) or Gaussian entries there is no reason for the PDF of the solution of a simulation to remain uniform or Gaussian.

The third and fourth moments, the skewness $\gamma$ and the kurtosis $\kappa$, are defined as:

$$\gamma = \mathbb{E}[(\frac{J(x, \alpha) - \mu}{\sigma})^3], \quad \kappa = \mathbb{E}[(\frac{J(x, \alpha) - \mu}{\sigma})^4]. \qquad (7)$$

The skewness is an indication of the deviation from symmetry of the PDF and the kurtosis permits to quantify if it is tall and skinny or short and squat.

One can consider that a robust design should favor symmetry in the distribution which means a lower absolute value of skewness. This is an assumption and other non-zero values can be targeted for the absolute value of the skewness. However, the assumption of symmetry appears reasonable. Indeed, when driving a car on a straight line, one expects the car to have the same behavior for small and symmetrical disruptions in holding the steering wheel by the driver. In a Gaussian distribution we have $\gamma = 0$. Also, in a normal distribution the mean and median coincide and if a PDF is not too far from a normal distribution, the median will be near $\mu - \gamma\sigma/6$. Therefore, if $|\gamma| \to 0$ the PDF tends toward a normal distribution. This provides an inequality constraint on $|\gamma|$ as $\gamma$ can be either positive or negative. For an uni-modal PDF a reduction of the skewness comes when the mean and the mode of the distribution converge toward each other at given standard deviation.

Robust design means higher density near the mean which means higher kurtosis, but this is more subtle. Indeed, despite higher kurtosis means concentration of the probability mass around the mean, it could also imply thicker tails in the PDF. This means that more of the variance is the result of infrequent extreme deviations. We need therefore to define what we mean by robust design: acceptance of frequent modest deviations or acceptance of infrequent extreme ones. If operational security is a major concern the latter should be definitely avoided. We therefore consider that a reasonable requirement would be to have a design reducing the initial kurtosis value: $\kappa \le \kappa_0$ together with a constraint on the variance $\sigma$.
To summarize, robust moment-based optimization can be seen through the following constrained minimization problem:

$$ \begin{gathered} \min_{\mathbf{x} \in O_{ad}} J(\mathbf{x}) = \mu(\mathbf{x}) \quad \text{such that} \\ \sigma(\mathbf{x}) - \sigma_0 \le 0, \quad |\gamma(\mathbf{x}) - |\gamma_0| \le 0, \quad \kappa(\mathbf{x}) - \kappa_0 \le 0. \end{gathered} \qquad (8) $$

Equality constraints on the moments are particular cases of these. In cases, higher kurtosis are targeted (see discussion above), the last constraint becomes $-\kappa(\mathbf{x}) + \kappa_0 \le 0$.

## 3.2 Discrete expressions

Monte Carlo simulations permit to recover these moments with an error decreasing as $\sigma/\sqrt{M}$ with M the number of functional evaluations and this with a convergence rate independent of N. But, for small N, classical numerical integration over-performs Monte Carlo simulations in term of complexity based on the number of functional evaluations to recover the moments at a given accuracy. As we are interested by small values of N, this latter may therefore be preferred.

Both Monte Carlo trials and numerical integration lead to the introduction of weighted sums [24] over a M-point sampling $I_M$ as estimators of the previous moments (denoted with the same notation):

$$ \mu = \sum_{\alpha_k \in I_M} \omega_k J(\mathbf{x}, \alpha_k), \qquad \sigma^2 = \sum_{\alpha_k \in I_M} \omega_k (J(\mathbf{x}, \alpha_k) - \mu)^2, \qquad (9) $$

$$ \gamma = \sum_{\alpha_k \in I_M} \omega_k (J(\mathbf{x}, \alpha_k) - \mu)^3, \qquad \kappa = \sum_{\alpha_k \in I_M} \omega_k (J(\mathbf{x}, \alpha_k) - \mu)^4. \qquad (10) $$

A major difference between Monte Carlo and numerical integration is that in the former the sampling $I_M$ and the weights are chosen according to the PDF of $\alpha$. In this paper we consider uniform PDF. We have therefore uniform sampling in each of the dimension of $\alpha$ and $\omega_k \sim 1/M$. Unbiased estimates use slightly different coefficients and also introduce corrections. For instance, one should consider $\omega_k = 1/(M - 1)$ for $\sigma$, $\omega_k = M/((M - 1)(M - 2))$ for $\gamma$ and $\kappa$. However, because we are interested by the gradients of these quantities with respect to $\mathbf{x}$ and gradient-based minimization algorithms, these changes will only have slight impacts on the amplitude of the gradient with no real incidence on the optimization history, especially if optimal descent step sizes are used. This discussion also shows that moment-based and multipoint optimization can be seen as a whole. If different PDFs are targeted than uniform, we need to redefine the sampling $I_M$ and the weights $\omega_k$ accordingly. But, this is beyond the scope of this paper.

## 3.3 Gradients of the moments

The linearity in the sums permits to access to the gradients of the moments with respect to the optimization variables $\mathbf{x}$ from the gradient of the functional at the sampling point $\alpha_k$. We recall
that we took $\sigma = 1$ in Equation (7) to ease this linearization. The different gradients read:

$$
\begin{align}
\nabla_{\mathbf{x}} \mu &= \sum_{\alpha_k \in \mathbf{I}_M} \omega_k \nabla_{\mathbf{x}} J(\mathbf{x}, \alpha_k), \\
\nabla_{\mathbf{x}} \sigma^2 &= \sum_{\alpha_k \in \mathbf{I}_M} 2\omega_k (J(\mathbf{x}, \alpha_k) - \mu) \nabla_{\mathbf{x}} J(\mathbf{x}, \alpha_k) - \mu, \tag{11} \\
\nabla_{\mathbf{x}} \gamma &= \sum_{\alpha_k \in \mathbf{I}_M} 3\omega_k (J(\mathbf{x}, \alpha_k) - \mu)^2 \nabla_{\mathbf{x}} J(\mathbf{x}, \alpha_k) - \mu, \\
\nabla_{\mathbf{x}} \kappa &= \sum_{\alpha_k \in \mathbf{I}_M} 4\omega_k (J(\mathbf{x}, \alpha_k) - \mu)^3 \nabla_{\mathbf{x}} J(\mathbf{x}, \alpha_k) - \mu.
\end{align}
$$

Knowing the elementary gradients $\nabla_x J(\mathbf{x}, \alpha_k)$ at sampling points $\alpha_k$ brings interesting information on the geometry of the problem. Indeed, they can be used to build a vector space and the analysis of this space provides important information on the complexity of the problem [24–26].

The natural presence of parallelism due to the *M* independent evaluations of the state, functional and its gradient is interesting to monitor computational complexity. In addition, as mentioned earlier for our problem, we notice that only the right-hand-sides $b_{i=1,...,N}$ need to be assembled after each perturbation of $\alpha = u^*$. Therefore, the optimality condition mentioned in section II for *J* can be extended to $\mu$ writing:

$$ \nabla_x \mu = \sum_{\alpha_k \in \mathbf{I}_M} \omega_k \nabla_x J(\mathbf{x}, \alpha_k) = 0, $$

which gives the following linear system to solve:

$$ A\mathbf{x} - \sum_{\alpha_k \in \mathbf{I}_M} \omega_k b_{\alpha_k} = A\mathbf{x} - \sum_{\alpha_k \in \mathbf{I}_M} \omega_k < (u_0 - u^*), u > = 0, $$

with the matrix $A_{i,j} = \langle u_i, u_j \rangle$ unchanged.

In the presence of constraints on moments, we can proceed in the same way to get the gradients of the moments, again using the nullity of the gradient of the first moment at the optimum:

$$
\begin{align*}
\nabla_{\mathbf{x}} \sigma^2 &= \sum_{\alpha_k \in \mathbf{I}_M} 2\omega_k (J(\mathbf{x}, \alpha_k) - \mu) \nabla_{\mathbf{x}} J(\mathbf{x}, \alpha_k) \\
&= C_1(\mathbf{x}) A \mathbf{x} - \sum_{\alpha_k \in \mathbf{I}_M} 2\omega_k (J(\mathbf{x}, \alpha_k) - \mu) b_{\alpha_k},
\end{align*}
$$

$$
\begin{align*}
\nabla_{\mathbf{x}} \gamma &= \sum_{\alpha_k \in \mathbf{I}_M} 3\omega_k (J(\mathbf{x}, \alpha_k) - \mu)^2 \nabla_{\mathbf{x}} J(\mathbf{x}, \alpha_k) \\
&= C_2(\mathbf{x}) A \mathbf{x} - \sum_{\alpha_k \in \mathbf{I}_M} 3\omega_k (J(\mathbf{x}, \alpha_k) - \mu)^2 b_{\alpha_k},
\end{align*}
$$
$$ 
\begin{aligned}
\nabla_{\mathbf{x}} \kappa &= \sum_{\alpha_k \in \mathbf{I}_M} 4\omega_k (J(\mathbf{x}, \alpha_k) - \mu)^3 \nabla_{\mathbf{x}} J(\mathbf{x}, \alpha_k) \\
&= C_3(\mathbf{x}) A \mathbf{x} - \sum_{\alpha_k \in \mathbf{I}_M} 4\omega_k (J(\mathbf{x}, \alpha_k) - \mu)^3 b_{\alpha_k}.
\end{aligned}
$$

We see that these are four vectors involving a same matrix A and elementary contributions $b_{\alpha_k}$ with different weightings and scaling $C_1, C_2$ and $C_3$ which are defined by:

$$ C_1(\mathbf{x}) = \sum_{\alpha_k \in \mathbf{I}_M} 2\omega_k (J(\mathbf{x}, \alpha_k) - \mu), $$

$$ C_2(\mathbf{x}) = \sum_{\alpha_k \in \mathbf{I}_M} 3\omega_k (J(\mathbf{x}, \alpha_k) - \mu)^2, $$

$$ C_3(\mathbf{x}) = \sum_{\alpha_k \in \mathbf{I}_M} 4\omega_k (J(\mathbf{x}, \alpha_k) - \mu)^3. $$

## 3.4 Controlling the moments

A classical approach to improve robust design in engineering is to control the first moment under the constraint of leaving unchanged or even reduce the second one. An easy way to adapt this concept to our situation is to consider the following descent direction orthogonal to the gradient of the variance:

$$ d = - \left( \nabla_{\mathbf{x}} \mu - \langle \nabla_{\mathbf{x}} \mu, \frac{\nabla_{\mathbf{x}} \sigma}{\| \nabla_{\mathbf{x}} \sigma \|} \rangle > \frac{\nabla_{\mathbf{x}} \sigma}{\| \nabla_{\mathbf{x}} \sigma \|} \right), $$

where $\nabla_{\mathbf{x}} \sigma / \| \nabla_{\mathbf{x}} \sigma \|$ is the unit vector along $\nabla_{\mathbf{x}} \sigma$. Unfortunately, this formulation does not permit anymore to get the optimal solution directly inverting a linear system as it is the case in the single-point or even when minimizing solely the first moment. An iterative procedure is then necessary using some descent algorithm. We will see examples of these below.

One can generalize the previous construction to the case of more than two moments using a Gram-Schmidt orthonormalization procedure. Let us organize the four vectors in the following order:

$$ (\nabla_{\mathbf{x}} \kappa, \nabla_{\mathbf{x}} \gamma, \nabla_{\mathbf{x}} \sigma^2, \nabla_{\mathbf{x}} \mu). \qquad (12) $$

To reduce the impact of the solution of the first order optimality condition for the first moment on the higher moments, we require the descent direction to be orthogonal to the subspace generated by the gradients of these moments. This can be easily achieved using a Gram-Schmidt orthonormalization which extracts from the set of the vectors above a new set of orthonormal vectors $v_{i=1,...,4}$:

$$ (v_1, v_2, v_3, v_4). $$

The descent direction we use in minimization is defined by:

$$ d = - \| \nabla_{\mathbf{x}} \mu \| v_4. $$
In this way, we make sure the descent direction is orthogonal to the subspace generated by moments 2, 3 and 4. The organization of the vectors as given in (12) is important to make sure the direction monitoring optimization is the one given by the first moment. Once this is achieved, this descent direction can be used in an iterative procedure which is again necessary as the constrained optimality condition cannot simply be cast into a single linear system.

The examples below show that this choice of search direction leads to a reduction of the second and the fourth moments as well. This is positively surprising from a robust optimization view point as one would have actually expected the higher moments to be rather conserved.

## 3.5 Multi-point optimization for a quadratic functional

Let us illustrate our ingredients on a simple example. The functional involves a least-squares minimization:

$$J(\mathbf{x}, \alpha) = \frac{1}{2} \sum_{i=1}^{N} (x_i - \alpha)^2, \quad -0.5 \le \alpha \le 0.5, \quad -5 \le x_{i=1,...,N} \le 5, \quad N = 40. \quad (13)$$

Let us solve problem (8) with a gradient method. Here, $\alpha$ is a scalar and for a given $\alpha$ the optimality condition for $J(\mathbf{x}, \alpha)$ gives obviously $\mathbf{x}^* = (\alpha, ..., \alpha) \in \mathbb{R}^N$. This would be the solution of a single-point optimization. Here $\alpha$ has a uniform PDF and we use a uniform sampling $I_M$ with $M = 40$. Minimizing the mean $\mu$ without any constraint on the other moments still produces a non robust optimum as the performance has large variability over the range of the operating condition $\alpha$. On the other hand, asking for the first moments to decrease under constraints on the other three produces a much flatter functional over the range of $\alpha$. Figure 3 shows the histograms of $J(\mathbf{x}^*, \alpha_{k=1,...,M})$ for these minimization. In particular, one sees how controlling more moments affects the left peak and the distribution of the frequencies.

## 3.6 Application to our 3D elastic problem

We now illustrate the multi-point discussion with moments on our 3D elasticity problem presented in section 2.1 where target displacement $u^*$ is assumed uncertain with its Probability Density Function (PDF) known. The target is a large vector with a size given by the number of the 3D mesh points. We consider centered Gaussian perturbations of $u^*$ with a standard deviation of 10%.

In our elasticity optimization problem, $O_{ad} = (\mathbb{R}^+)^N$. This means that each of the optimization variable needs to be positive. This brings additional constraints to those expressed earlier on the moments. This also means that one cannot solve anymore the single-point optimization by a sole inversion of a linear system and that, for instance, a projected gradient algorithm is necessary. This can be easily included during the iterations when minimizing the first moment under constraints on the other three.

Figure 4 shows the distribution of the ratio of the functional $J(\mathbf{x}^*, \alpha_k)$ to the Young modulus for $M = 1000$. One sees that the uncertainty on $u^*$ brings large variations in $J$. This variability shows the importance of a robust optimization approach for this problem. The impact of considering the variability during minimization through moment based functional can be easily demonstrated looking at the moment-based solutions of the minimization versus the solution of the single-point optimization for the target $u^*$ which corresponds to the deterministic situation. Table 1 features the solution of the optimization variables for this comparison. We see that the results are clearly different. Also the histograms in Figure 4 show that controlling all four
<figure><img src="image_3.png" /><figcaption>Figure 3: Upper: functional $J(\mathbf{x}^*, \alpha)$ given by Equation (13) vs. $\alpha$ for single-point optimization, mean-based and with constraints on moments two, three and four. Lower: histograms of $J(\mathbf{x}^*, \alpha_{k=1,...,M})$ for a single-point minimization and when controlling one and then all first four moments.</figcaption></figure>

moments greatly reduces the variability of the functional and therefore improves the robustness of the design. This can be confirmed considering the 4 moments showing the variability of the functional induced by the variability of the target displacement, considered as uncertain, for the single-point optimal solution and for the different optimal solutions obtained when controlling respectively the first, second and eventually all four moments. The first column in Table 2 show the first moment for $J(\mathbf{x}, \alpha)/E$ and the other columns show the second, third and fourth moments for $(J(\mathbf{x}, \alpha) - \mu)/E$. By extension, we denote these $\mu(J(\mathbf{x}, \alpha)/E)$, $\sigma((J(\mathbf{x}, \alpha) - \mu)/E)$, $\gamma((J(\mathbf{x}, \alpha) - \mu)/E)$ and $\kappa((J(\mathbf{x}, \alpha) - \mu)/E)$. The first column indicates that the third first designs give about the same average behavior while controlling all four moments clearly provide a different, and actually better, solution. Robustness improvements can be seen from the evolution of the second and fourth moments where again controlling all four moments is recommended. Also, the final density functions are slightly less symmetric when controlling more moments. However, one regain symmetry when controlling 4 moments instead of only two. As
we expressed in section 3.1, symmetric distributions should be preferred in robust designs.

<table><caption>Table 1: Optimization variables for minimization using the deterministic target $u^*$ (the single-point case) and when controlling 1, 2 and all 4 moments of $J$. The solutions are clearly different showing the impact of moment-based optimization. The robustness achieved with more moments accounted for in the design is well illustrated through the histograms in Figure 4.</caption><thead><tr><th></th><th>Single-point</th><th>1 moment</th><th>2 moments</th><th>4 moments</th></tr></thead><tbody><tr><th>x<sub>1</sub> (KN)</th><td>5500</td><td>5496</td><td>5266</td><td>7035</td></tr><tr><th>x<sub>2</sub> (kN)</th><td>4000</td><td>4000</td><td>3917</td><td>9262</td></tr><tr><th>x<sub>3</sub> (KN)</th><td>6000</td><td>5994</td><td>5518</td><td>6834</td></tr><tr><th>x<sub>4</sub> (KN)</th><td>4500</td><td>4503</td><td>4330</td><td>0</td></tr><tr><th>x<sub>5</sub> (KN)</th><td>4999</td><td>5005</td><td>5594</td><td>1772</td></tr></tbody></table>

<table><caption>Table 2: Distributions of the different moments for the optimal solutions obtained after a single-point optimization and when controlling first, second and eventually all four moments. This final design is superior both in performance and robustness.</caption><thead><tr><th></th><th>μ(J/E)</th><th>σ((J - μ)/E)</th><th>γ((J - μ)/E)</th><th>κ((J - μ)/E)</th></tr></thead><tbody><tr><th>single point</th><td>74.980</td><td>0.938</td><td>-0.005</td><td>2.422</td></tr><tr><th>1 moment</th><td>74.940</td><td>0.865</td><td>-0.018</td><td>2.180</td></tr><tr><th>2 moments</th><td>74.144</td><td>0.840</td><td>0.094</td><td>2.146</td></tr><tr><th>4 moments</th><td>64.138</td><td>0.765</td><td>0.032</td><td>1.685</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>Figure 4: Relative frequency histogram of $(J(\mathbf{x}^*, \alpha_{k=1,...,M}) - \mu)/E$ with $M = 1000$. Controlling all four moments reduces the variability of the functional and therefore improves the robustness of the design.</figcaption></figure>
# IV LINEAR MODELS AND BACKWARD UNCERTAINTY PROPAGATION: COVARIANCE MATRIX OF THE OPTIMIZATION VARIABLES

We showed how to quantify the impact of target state uncertainties in inverse design through multi-point analysis and showed how it can be monitored through the moments of the functional. But the multi-point analysis requires a sampling of a large dimensional space. Something which is a burden when this dimension is large. To avoid a sampling we would like to see if an alternative low-complexity approach based on the existence of a local linear relationship between the target state $u^*$ and the optimization variables $x$ can be applied to estimate the covariance matrix $Cov_x^2$ of the optimization variables. We compare this construction to a reference $Cov_x^1$ construction where a sampling is necessary.

As previously, we assume the PDF of the target solution $u^*$ known and therefore its covariance matrix $Cov_{u^*}$. Also, we assume $u^*$ is admissible and that there exists a set of optimization variables realizing the target.

As mentioned in (1), in our problem of interest the state equation is linear and to simplify the notation, and without any loss of generality, we consider the source term $f(x) = x$ considering that $x$ and $u$ have a same dimension. In our implementation, however, we use spatial localization terms given in (2) in order to reduce the size of $x$.

The construction can be applied in several useful situations. Let us describe two of them.

## 4.1 $Cov_x^1$ when a linear model $L_1 u = x$ is available

This is typically our case, for instance, with the elasticity equation and the optimization variables acting on the right-hand-side of the equation. In this case, we can express $Cov_x^1$ with $u \sim u^*$, through:

$$ 
\begin{aligned}
Cov_x^1 &= \mathbb{E}((x - \bar{x})(x - \bar{x})^T) \\
&= L_1 \mathbb{E}((u^* - \bar{u}^*)(u^* - \bar{u}^*)^T) L_1^T \\
&= L_1 Cov_{u^*} L_1^T.
\end{aligned}
$$

But, in our case, the solver is a black-box and we do not have $L_1$ in hand. In these cases, a linear model can be built through a machine learning approach using a database of simulation scenarios couples of $(x - u)$ we denote by $(X - U)$ generated by this code. This procedure can be applied in nonlinear situations as well. In these cases, as the uncertainty analysis is aimed to take place around a given/optimal solution, it is enough for the surrogate model to have a local domain of validity.

$L_1$ is solution of the first-order optimality condition for a functional $J(L_1)$:

$$ \nabla_{L_1} J(L_1) = 0 \Leftrightarrow UU^T L_1^T = U X^T, $$

where $U$ and $X$ gather the scenarios of simulations available in the database for $u$ and $x$ and we look for the best linear relation $X^T = U^T L_1^T$ linking the scenarios in $U$ to those in $X$ and, actually, vice versa. $L_1^T$, and therefore $L_1$, can be found as solution of a least-squares problem:

$$ J(L_1) = (X^T - U^T L_1^T, X^T - U^T L_1^T) = (X, X) - 2(XU^T, L_1) + (L_1, L_1)_{UU^T}. $$
The only requirement is for $U$ to have full rank in order for $UU^\top$ to be positive definite such that the scalar product $(L_1, L_1)_{UU^\top} \ge 0$. Our formulation presented in section II, and based on the introduction of a database of elementary solutions, can be seen as a particular case of machine learning taking advantage of the linearity to express $x = A^{-1}b(u^*)$. We recall that the components of $A$ and $b$ for $i, j = 1, ..., N$ are given by:

$$A_{i,j} = \langle u_i, u_j \rangle \quad \text{and} \quad b_i(u^*) = \langle (u_0 - u^*), u_i \rangle.$$

Again, the covariance matrix of the optimization variables reads:

$$\begin{aligned}
\mathrm{Cov}_{\mathbf{x}}^{1} &= \mathrm{E} \left( (\mathbf{x} - \bar{\mathbf{x}})(\mathbf{x} - \bar{\mathbf{x}})^{\top} \right) \\
&= A^{-1} \mathrm{E} \left( (b(u^*) - b(\bar{u}^*)) (b(u^*) - b(\bar{u}^*))^{\top} \right) A^{-\top} \\
&= A^{-1} \mathrm{Cov}_{b(u^*)} A^{-\top},
\end{aligned} \quad (14)$$

where it is easy to estimate $\mathrm{Cov}_{b(u^*)}$ knowing $\mathrm{Cov}_{u^*}$. The elementary $u_i$ and the homogeneous solutions are deterministic. The covariance matrix of $b(u^*)$ measures the covariance of the projection of the deviation between the uncertain target function $u^*$ from the homogeneous solution $u_0$ on the elementary solutions $u_i$:

$$\begin{aligned}
(\mathrm{Cov}_{b(u^*)})_{ij} &= \mathrm{E} \left( (b_i(u^*) - \overline{b_i(u^*)}) (b_j(u^*) - \overline{b_j(u^*)})^{\top} \right) \\
&= \mathrm{E} \left( \langle (\overline{u^*} - u^*), u_i \rangle \langle (\overline{u^*} - u^*), u_j \rangle^{\top} \right).
\end{aligned}$$

One interest of this formulation is that the size of this matrix is now five, the number of the optimization variables, instead of 14400, the number of the mesh grid points for $\mathrm{Cov}_{u^*}$. Also, if the perturbation on $u^*$ are centered, $\overline{u^*}$ is the deterministic target solution.

## 4.2 $\mathrm{Cov}_{\mathbf{x}}^2$ when a linear model $L_2\mathbf{x} = u$ is available

This is not directly the situation in this paper as the optimization variables intervene in the left-hand-side. However, this analysis is powerful when the state equation is linear, or if the Jacobian of the state with respect to the optimization can be estimated which is our situation where indeed the Jacobian is explicitly available knowing the elementary solutions as described in section (II).

After a deterministic optimization when the state is close to the target $(u \sim u^*)$, we assume the covariance matrix of $u$ close to the covariance matrix of $u^*$. Using the linear relationship $L_2\mathbf{x} = u$ we can write:

$$\begin{aligned}
\mathrm{Cov}_{u^*} &= \mathrm{E} \left( (u^* - \overline{u^*})(u^* - \overline{u^*})^{\top} \right) \\
&= \mathrm{E} \left( L_2 (\mathbf{x} - \bar{\mathbf{x}})(\mathbf{x} - \bar{\mathbf{x}})^{\top} L_2^{\top} \right) \\
&= L_2 \mathrm{E} \left( (\mathbf{x} - \bar{\mathbf{x}})(\mathbf{x} - \bar{\mathbf{x}})^{\top} \right) L_2^{\top} \\
&= L_2 \mathrm{Cov}_{\mathbf{x}}^2 L_2^{\top}.
\end{aligned}$$

Unlike the previous situation, to get $\mathrm{Cov}_{\mathbf{x}}^2$ we need to invert this expression. In cases where exact inversion is impossible a least-squares formulation can be used looking for $\mathrm{Cov}_{\mathbf{x}}^2$ minimizing:

$$\frac{1}{2} \langle L_2 \mathrm{Cov}_{\mathbf{x}}^2 L_2^{\top}, L_2 \mathrm{Cov}_{\mathbf{x}}^2 L_2^{\top} \rangle - \langle \mathrm{Cov}_{u^*}, L_2 \mathrm{Cov}_{\mathbf{x}}^2 L_2^{\top} \rangle.$$
First-order optimality condition with respect to $Cov_x$ gives:

$$L_2^\top L_2 \ Cov_x^2 L_2^\top L_2 - L_2^\top Cov_{u^*} L_2 = 0,$$

which leads to

$$Cov_x^2 = (L_2^\top L_2)^{-1} L_2^\top Cov_u L_2 (L_2^\top L_2)^{-1},$$

and eventually, to

$$Cov_x^2 = L_2^{-1} Cov_{u^*} L_2^{-\top} = (L_2^\top Cov_{u^*}^{-1} L_2)^{-1}. \quad (15)$$

Therefore, knowing $Cov_{u^*}$ one can access $Cov_x^2$ if a linear relationship, even locally valid, is available between $x$ and $u$. In practice, $Cov_{u^*}$ is diagonal as the target uncertainties are often assumed independent. The inversion of $Cov_{u^*}$ is therefore straightforward. The second expression only requires the inversion of a $N \times N$ matrix in our case. This expression in (15) is also interesting because it involves the inversion of a square matrix giving a least-squares sense to the inversion of rectangular matrix. An interesting application of this construction is using the state sensitivity with respect to optimization variables or the Jacobian $\mathcal{J} = (\partial u / \partial x)$ in a first-order linear relationship $\mathcal{J} \delta x = \delta u$. This permits to estimate the covariance matrix of the optimization variables perturbations even in the absence of a direct linear model linking $x$ and $u$. Actually, in situations where the state equation is linear (our case) we have $L_2 = \mathcal{J}$. In situations where there are only a few optimization variables, the Jacobian can simply be estimated by a finite differences approach approximating the $i^{th}$ line of the Jacobian by:

$$\mathcal{J}_{i=1,...,N} \sim \frac{(u(\mathbf{x}e_i + \varepsilon) - u(\mathbf{x}e_i - \varepsilon))}{2\varepsilon},$$

with $e_{i=1,...,N}$ the canonical basis of $\mathbb{R}^N$ in our case and $0 < \varepsilon << 1$. But, we have seen in section II that $u$ has a direct expression with respect to $x$ components $x_i$ through:

$$u = u_0 + \sum_{i=1}^{N} x_i u_i.$$

Therefore, the finite difference formula reduces to:

$$\mathcal{J}_{i=1,...,N} \sim \frac{(u(\mathbf{x}e_i + \varepsilon) - u(\mathbf{x}e_i - \varepsilon))}{2\varepsilon} = u_i.$$

This is another advantage of our direct formulation for the inversion as it makes evaluating $\mathcal{J}$ straightforward, not requiring any extra calculation. In addition, the Jacobian is exact in this case and is not affected by the finite differences approximation and the choice of $\varepsilon$. Actually, the Jacobian would be exact even using first order finite differences. This is because we recover the exact gradient expression in this case. Therefore, this second covariance estimation is well adapted to our situation and would be interesting if it gives nearly the same answer than the first construction where a sampling of the state variable variation domain is necessary.

Figure 5 shows the comparison of this two estimations of $Cov_x$. We see that the matrices are very close and especially the signs of all components agree.

Further comparisons are given in Figure 3 through different distances between the two covariance matrices versus the distance between the covariance matrices and a diagonal matrix with the diagonal given by mean values of the covariance matrices components (i.e. $\frac{1}{2N^2} \sum_{i,j=1}^N (Cov_x^1 + Cov_x^2)_{i,j}$). This is to provide an estimation of how close the two covariance matrices are versus a third fixed one.
<figure><img src="image_5.png" /><figcaption>Figure 5: Covariance Matrix $Cov_x^1$ and $Cov_x^2$ (in Newton$^2$) Visualization, left: $Cov_x^1$ from (14), right: $Cov_x^2$ from (15).</figcaption></figure>

<table><caption>Table 3: Different estimations of distances between the two covariance matrices. where $C = \frac{1}{2N^2} \sum_{i,j=1}^{N} (Cov_x^1 + Cov_x^2)_{i,j}$ is the mean values of the covariance matrices components and $I_N$ is the $N^{th}$ identity matrix. One notices that, regardless of the considered distance, $Cov_x^1$ and $Cov_x^2$ are about at the same distance of $CI_N$. Also, $Cov_x^1$ are $Cov_x^2$ are at a distance about ten times less except for the Kullback distance where two matrices are seen much closer. This is very interesting as this distance is particularly suitable for the comparison between probability density functions and covariance matrices.</caption><thead><tr><th>Distance d</th><th>d(Cov<sub>x</sub><sup>1</sup>, Cov<sub>x</sub><sup>2</sup>)</th><th>d(Cov<sub>x</sub><sup>1</sup>, CI<sub>N</sub>)</th><th>d(Cov<sub>x</sub><sup>2</sup>, CI<sub>N</sub>)</th></tr></thead><tbody><tr><td>Kullback distance</td><td>0.334</td><td>165.589</td><td>166.682</td></tr><tr><td>Euclidean distance</td><td>1.674</td><td>17.020</td><td>16.123</td></tr><tr><td>Log Euclidean distance</td><td>0.737</td><td>7.193</td><td>7.469</td></tr><tr><td>Riemannian distance</td><td>0.811</td><td>7.193</td><td>7.469</td></tr><tr><td>Wasserstein distance</td><td>0.389</td><td>4.098</td><td>4.097</td></tr><tr><td>Log-det distance</td><td>0.286</td><td>2.00</td><td>2.086</td></tr></tbody></table>
# V CONCLUDING REMARKS

Robust parametric optimization has been addressed through the control of the moments of the functional. Moments have been estimated using a multi-point formulation. This permits to address situations where operating parameters are not anymore single valued but defined through their PDF. Hence, the first four moments of the functional have been considered to go beyond the classical mean-variance based optimization. Subtleties of what should be the target for the skewness and the kurtosis of the probability density function of the performance of the system for a robust design have been discussed.

It has been shown that controlling second, third and fourth moments drastically improves the design while the time to the solution remains comparable to the single-point situation as all the extra calculations can be carried out in a fully parallel and independent manner. The cost is also reduced thanks to our construction of the optimal solution taking advantage of the presence of the linearity between the optimization variables and the solution of the state equation. As a consequence, the approach requires quite small additional coding and computational effort.

Low-complexity solutions for backward propagation of aleatory uncertainty in target data has also been presented. Derivation of the covariance matrix of the optimization parameters has been discussed through two formulations. These provide uncertainty quantification analysis for the inversion solution with confidence margins on the design parameters in very large design spaces. It has been shown that taking again advantage of existing linear relationships between the target state $u^*$ and the optimization variables x, samplings of the design and the target solution spaces can be avoided drastically reducing the computational cost of the approach.

We recall that the approach is especially suitable in situations where linear dependencies exist between the state and the optimization variables. However, if this is not the case, for instance if the thickness distribution is the optimization variable, the approach can still be applied to the linearized problem around a given optimization variable distribution. Indeed, consider a particular set of optimization variable $x^*$, and $u^*(x^*)$ solution of a nonlinear equation $F(u^*(x^*)) = 0$. The linearized problem reads $F(u^*(x^*) + v) \sim F(u^*(x^*)) + D_x F(u^*(x^*))$ $v = 0$ giving $D_x F(u^*(x^*))v = -F(u^*(x^*))$ and what presented in the paper can be applied as $v$ linearly depends on the right-hand-side. As a consequence, the outcome of the analysis presented here can be seen as stability analysis for this particular design.

# REFERENCES

## Publications

[1] R. Kalman. "A New Approach to Linear Filtering and Prediction Problems". In: *Transactions of the ASME - Journal of Basic Engineering* 82 (1960), pages 35-45.

[2] P. Hoel. "Introduction to Mathematical Statistics". John Wiley, London, 1971.

[3] H. Lindman. "Analysis of Variance in Complex Experimental Designs". Freeman, New York, 1974.

[4] A. Tarantola. "Inverse problem theory and methods for model parameter estimation". SIAM, N.Y., 1987.

[5] R. Ghanem and P. Spanos. "Stochastic Finite Elements: A Spectral Approach". Springer Verlag, New York, 1991.
[6] M. G. K. Ide P. Courtier and A. Lorenc. "Unified notation for data assimilation: operational, sequential and variational". In: *Journal of the Meteorological Society of Japan* 75/1B (1997), pages 181–189.

[7] K. Kapur and L. Lamberson. "Reliability in Engineering Design". John Wiley & Sons, New York, 1997.

[8] AIAA. "Guide for the Verification and Validation of Computational Fluid Dynamics Simulations". In: *AIAA G-077* (1998).

[9] G. Casella and R. Berger. "Statistical Inference, 2nd Ed". Duxbury Press, London, 2001.

[10] K. S. A. Conn and L. Vicente. "Introduction to Derivative-Free Optimization". SIAM, NY, 2002.

[11] R. O. O. J. T. Spooner M. Maggiore and K. Passino. "Stable Adaptive Control and Estimation for Nonlinear Systems: Neural and Fuzzy Approximator Techniques". John Wiley, New York, 2002.

[12] M. Schnellenbach-Held and K. Pfeffer. "Punching behavior of biaxial hollow slabs". In: *Elsevier Science* (2002).

[13] H.-J. Bungartz and M. Griebel. "Sparse Grids". In: *Acta Numerica* 13 (2004), pages 147–269.

[14] R. Ghanem and A. Doostan. "On the construction and analysis of stochastic models: characterization and propagation of the errors associated with limited data," in: *J. of Comput. Phys.* 217 (2006), pages 63–81.

[15] P. Jorion. "Value at Risk: The New Benchmark for Managing Financial Risk". McGraw-Hill, New York, 2006.

[16] X. Wan and G. Karniadakis. "Multi-element generalized polynomial chaos for arbitrary probability measures". In: *SIAM J. Sci. Comput.* 28/3 (2006), pages 901–928.

[17] H. A. V. d. V. W. H. Schilders and J. Rommes. "Model order reduction: Theory, research aspects and applications". Volume 13. Springer Math in Industry series, Berlin, 2008.

[18] C. Geuzaine and J.-F. Remacle. "Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities". In: *International Journal for Numerical Methods in Engineering* (2009).

[19] D. Xiu. "Numerical Methods for Stochastic Computations: A Spectral Method Approach". Princeton University Press, 2010.

[20] S. S. C. Schillings and V. Schulz. "Efficient shape optimization for certain and uncertain aerodynamic design". In: *Computers & Fluids* 46(1) (2011), pages 78–87.

[21] J. Caers. "Modeling Uncertainty in the Earth Sciences". Wiley-Blackwell, 2011.

[22] L. F. E. "Présentation et utilisation de cast3m". In: *ENSTA-LME* (2011).

[23] M. P. J. Chery B. Mohammadi and C. Joulain. "Plate rigidity inversion in southern California using interseismic GPS velocity field". In: *Geophys. J. Int.* 187/2 (2011), pages 783–796.

[24] B. Mohammadi. "Reduced sampling and incomplete sensitivity for low-complexity robust parametric optimization". In: *Int. J. Num. Meth. Fluids* 73/4 (2013), pages 307–323.

[25] B. Mohammadi. "Principal angles between subspaces and reduced order modeling accuracy in optimization". In: *Structural and Multidisciplinary Optimization* 50/2 (2014), pages 237–252.

[26] B. Mohammadi. "Uncertainty Quantification by geometric characterization of sensitivity spaces". In: *Compt. Meth. Appl. Mech. Eng.* 280 (2014), pages 197–221.

[27] B. Mohammadi. "Value at Risk for confidence level quantifications in robust engineering optimization". In: *Optimal Control: Applications and Methods* 35/2 (2014), pages 179–190.
[28] B. Mohammadi. "Ensemble Kalman Filters and geometric characterization of sensitivity spaces for Uncertainty Quantification in optimization". In: *Computer Methods in Applied Mech. & Eng.* 290 (2015), pages 228-249.

[29] B. Mohammadi. "Parallel reverse time integration and reduced order models". In: *J. of Computational Mathematics* 2 (2015), pages 17–33.

[30] B. Mohammadi. "Backward uncertainty propagation in shape optimization". In: *Int. J. for Numerical Methods in Fluids* 103(4)-DOI: 10.1002/fld.4077 (2016), pages 307–323.

[31] B. Mohammadi. "Controlling first four moments for robust optimization". In: *Optimization and Engineering* 18-3 (2017), pages 561–585.