<figure><img src="image_1.png" /></figure>

HAL
open science

# Multifidelity kriging meta-model based on discretization error bounds

Ludovic Pierre Jérôme Mell, Valentine Rey, Franck Schoefs

► To cite this version:

Ludovic Pierre Jérôme Mell, Valentine Rey, Franck Schoefs. Multifidelity kriging meta-model based on discretization error bounds. International Journal for Numerical Methods in Engineering, 2020, 10.1002/nme.6451. hal-02944220

HAL Id: hal-02944220

https://hal.science/hal-02944220

Submitted on 21 Sep 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Multifidelity kriging meta-model based on discretization error bounds

Ludovic Mell; Valentine Rey, Franck Schoefs

September 21, 2020

keywords: Reliability analysis; Probability of failure; Error estimation; kriging methods; Finite Element Method

## 1 Abstract

This paper presents an approach to build a multi-fidelity kriging metamodel from finite element computations on different meshes for structural reliability assessment. The proposed method takes advantage of the computation of bounds on the discretization error, which enables to guarantee the state (safe or failure) of each computation of the performance function. An algorithm to build the meta-model from the different levels of fidelity and estimate the failure probability is provided. Illustrations are presented on a two dimensional mechanical crack opening problem. Bounds on the failure probability are also post-processed.

## 2 Introduction

Uncertainties may arise from a lack of knowledge (for example, characterization of material properties can be improved) or from inherent variability (an example is the load due to wind flow on the blade of a wind turbine).

As a consequence, designing structures subjected to uncertainty in a deterministic framework becomes impractical. As aftermath in case of failure are sometimes dramatic, the reliability analysis becomes crucial and sensitive. This safety analysis is usually done in a probabilistic approach that aims at estimating quantities such as the probability of failure, reliability indexes, sensitivity factors, and so on.

In the context of virtual testing for numerical certification of structures, the probability of failure is usually computed using numerical simulations of a model of the structure. Different strategies can be considered. First, sampling techniques
that mainly rely on Monte Carlo or quasi Monte Carlo simulations [4] present the
advantage of being non intrusive. Besides, their efficiency does not depend on the
number of random variables, making them attractive. Their main disadvantage is
that the convergence rate against the number of samples is very low, so that they
are computationally expensive. To tackle this issue, reduction of variance or multi
level Monte Carlo approach have been proposed [19, 20]. Secondly, approximation
methods aim at searching for a random field in a given approximation space, the
Galerkin method [33] is an example of them. When the geometry of the structure
is random, methods such as XSFEM have been developed [36, 38]. However, those
techniques can be intrusive.

The reduction of computational costs can also be done using surrogate mod-
els. First, polynomial response surfaces were developed [47]. Then, kriging-based
metamodels [24, 45] became very popular in the field of safety analysis. The most
popular method was developed in [11]. It was then adapted to very low probability
using important sampling or subset sampling [12, 51] and finally to system relia-
bility [14]. The main asset of kriging is that it offers an indicator of the quality of
the meta-model. Indeed this method considers that the meta-model is a realiza-
tion of a random Gaussian process, the standard deviation of this process is then
computable at each point where the meta-model is evaluated.

Two components of the mechano-probabilistic problem may affect the accu-
racy of the reliability assessement : epistemic uncertainties on some parameters
(including model parameters) and the inaccuracy of the numerical solution to the
mechanical problem. To deal with epistemic uncertainties on parameter [54, 53]
propose a projection method to propagate through the mechanical model bounds
on these variable to obtain bounds on the probability of failure. Few people con-
sider the error due to the discretization method (mainly the finite element method).
Yet, this error may lead to a poor estimation of the probability of failure, as il-
lustrated in [41]. Performing numerical simulations on too coarse meshes might
in fact lead to considering that a realization belongs to the safety domain while it
actually belongs to the failure domain. Among works tackling this issue, [2] can be
cited. It uses the discretization error introduced by meshing using the Richardson
extrapolation in a non-deterministic analysis. In [30], the Richardson extrapo-
lation is also used to perform reliability analysis with the First Order Reliability
Method (FORM) [22]. In [16], the author estimates the probability of failure using
FORM while utilizing the discretization error to provide bounds on the probability
of failure.

In the family of kriging techniques, multi-fidelity kriging consists in using sev-
eral models with different precision to build a multi-fidelity meta-model. The main
objective of this approach is to obtain an improved trade-off between computa-
tional cost and precision by choosing to evaluate some samples on a low fidelity
model and some others on a high fidelity one. In [7], optimization in the context
of a fluid flow problem is performed with progressively fine grids. In [8, 39], multi-
fidelity kriging is performed on a flow problems. In the context of optimization,
[5, 6] adapts the stopping criterion of the iterative contact solver to use multi-
fidelity kriging. In [35], a multi-fidelity approach is used on a elasto-viscoplasticity
problem with different proper generalized decompositions. It is to be noted that
usually, the choice of levels of discretization or grids is done a priori and the
construction of the design of experiments is based on heuristics.

In this paper, we propose an adaptive multi-fidelity kriging approach to esti-
mate the probability of failure. This method uses the estimation of the discretiza-
tion error to choose a suitable mesh for each sampling point by noting that only
the sign of the limit state function is of importance : the reliability problem is
considered as a classification problem. In addition to the estimation of the proba-
bility of failure, the method gives lower and upper bounds on the exact probability
of failure.

The paper is organized as follow. In section 3, we define the mechanical prob-
lem and recall the kriging-based meta model approach used to estimate the failure
probability. In section 4, we present the approach we propose to build a multi-
fidelity kriging meta-model and its application to compute the probability of fail-
ure. In section 5 we apply our method to a crack-opening problem. Finally, section
6 concludes the paper.

3 Settings

In this section, the continuous mechanical problem and its discretized version
are first defined. Then, the discretization error and its estimates are presented.
Finally, we define the limit-state functions and formulate the reliability assessment
in terms of estimation of the probability of failure.

3.1 Continuous problem

Let $\mathbb{R}^d$ represent the physical space. Let us consider the static equilibrium of a (polyhedral) structure occupying the open domain $\Omega \subset \mathbb{R}^d$ and subjected to a given body force $f$ within $\Omega$, to a given traction force $\underline{F}$ on $\partial_F \Omega$ and to a given displacement field $\underline{u}_d$ on the complementary part $\partial_u \Omega \neq \emptyset$. Let us assume the structure undergoes small perturbations and that the material is linear elastic, characterized by Hooke's elasticity tensor $\mathbb{H}$. Let $\underline{u}$ be the unknown displacement field, $\underline{\varepsilon}(\underline{u})$ the symmetric part of the gradient of $\underline{u}$, $\underline{\sigma}$ the Cauchy stress tensor.

Two affine subspaces and a positive form are introduced:
* The affine subspace of kinematically admissible fields (KA-fields)

$$KA = \{ \underline{u} \in (\mathbb{H}^1(\Omega))^d, \underline{u} = \underline{u}_d \text{ on } \partial_u \Omega \} \quad (1)$$

and we note $KA^0$ the associated vectorial space.

* Affine subspace of statically admissible fields (SA-fields)

$$SA = \left\{ \underline{\tau} \in (L^2(\Omega))_{\text{sym}}^{d \times d}; \forall \underline{v} \in KA^0, \right. \\ \left. \int_{\Omega} \underline{\tau} : \underline{\varepsilon}(\underline{v}) d\Omega = \int_{\omega} \underline{f} \cdot \underline{v} d\Omega + \int_{\partial_F \Omega} \underline{F} \cdot \underline{v} dS \right\} \quad (2)$$

* Error in constitutive equation

$$e_{CR_{\Omega}}(\underline{u}, \underline{\sigma}) = \|\underline{\sigma} - \mathbb{H} : \underline{\varepsilon}(\underline{u})\|_{\mathbb{H}^{-1}, \Omega} \quad (3)$$

$$\text{where } \|\underline{x}\|_{\mathbb{H}^{-1}, \Omega} = \sqrt{\int_{\Omega} (\underline{x} : \mathbb{H}^{-1} : \underline{x}) d\Omega}$$

The continuous problems reads:

$$\left\{ \begin{array}{l} \text{Find a displacement field } \underline{u} \text{ and a stress field } \underline{\sigma} \text{ such that} \\ \underline{u} = \underline{u}_d \text{ on } \partial\Omega \cap \partial\omega \text{ and } \underline{\varepsilon}(\underline{u}) = \frac{1}{2}(\underline{\underline{\text{grad}}}(\underline{u}) + \underline{\underline{\text{grad}}}^T(\underline{u})) \text{ on } \Omega \\ \underline{\text{div}}(\underline{\sigma}) + \underline{f} = \underline{0} \text{ on } \Omega \text{ and } \underline{\underline{\sigma}}\underline{n} = \underline{F} \text{ on } \partial_F \Omega \\ \underline{\underline{\sigma}} = \mathbb{H} : \underline{\varepsilon}(\underline{u}) \text{ on } \Omega \end{array} \right. \quad (4)$$

The solution to this problem exists and is unique. We will denote the couple displacement field and stress field as the exact solution.

The mechanical problem can also be formulated as:

$$\text{Find } (\underline{u}_{ex}, \underline{\sigma}_{ex}) \in KA \times SA \text{ such that } e_{CR_{\Omega}}(\underline{u}_{ex}, \underline{\sigma}_{ex}) = 0$$

## 3.2 Discrete problem

Let $\Omega_H$ be a tessellation of $\bar{\Omega}$ by triangles. The finite element method consists in searching for a displacement field in the finite subspace $KA_H$ of KA where $KA_H$ reads:

$$KA_H = \{ \underline{u} \in (\mathbb{H}^1(\Omega))^d, \underline{u} = \underline{u}_d \text{ on } \partial_u \Omega_H \} \quad (5)$$
The discrete problem can be formulated as:

$$ 
\begin{gathered}
\underline{u}_H \in KA_H \\
\underline{\sigma}_H = \mathbb{H} : \underline{\varepsilon}(\underline{u}_H) \\
\int_{\Omega_H} \underline{\sigma}_H : \underline{\varepsilon}(\underline{v}_H) d\Omega = \int_{\Omega_H} \underline{f} \cdot \underline{v}_H d\Omega + \int_{\hat{\partial}_F \Omega_H} \underline{F} \cdot \underline{v}_H dS
\end{gathered}
\qquad (6) $$

The solution of this discrete problem exists and is unique. However, the discrete solution $\underline{u}_H$ usually does not coincide with the continuous exact solution $\underline{u}$.

## 3.3 Estimation of the discretization error

### 3.3.1 Generalities

Since the discrete displacement field $\underline{u}_H$ is not the exact solution, the finite element method introduces an error usually called discretization error. Results on the convergence of the FEM offer an *a priori* estimate of this error. It involves problem dependent constants that are not computable thus making this error estimate impractical. One can also use *a posteriori* error estimator that rely on a post-process of the finite element solution to derive an estimation of the discretization error [1]. In this paper, estimators based on the error in constitutive relation (3) were chosen as they provide guaranteed error bounds [28]. These estimators are based on the fundamental relation:

$$ 
\begin{gathered}
\forall (\hat{\underline{u}}, \hat{\underline{\sigma}}) \in KA(\Omega) \times SA(\Omega), \\
\|\underline{\varepsilon}(\underline{u}_{ex}) - \underline{\varepsilon}(\hat{\underline{u}})\|_{\mathbb{H},\Omega}^2 + \|\underline{\sigma}_{ex} - \hat{\underline{\sigma}}\|_{\mathbb{H}^{-1},\Omega}^2 = e_{CR\Omega}^2(\hat{\underline{u}}, \underline{\sigma})
\end{gathered}
\qquad (7) $$

We note $\|\underline{v}\|_{\Omega} = \|\underline{\varepsilon}(\underline{v})\|_{\mathbb{H},\Omega}$ the energy norm of the displacement. By choosing $\hat{\underline{u}} = \underline{u}_H \in KA(\Omega)$, we obtain the following upper bound for the error $e_{discr} = \underline{u}_{ex} - \underline{u}_H$:

$$ e_{discr} := \|\underline{e}_{discr}\|_{\Omega} \le e_{CR\Omega}(\underline{u}_H, \underline{\sigma}) \qquad (8) $$

Computing a statically admissible stress field is a complex task but several methods have been developed to compute $\underline{\sigma} \in SA(\Omega)$ (see [27], [37], [40] and [42]).

However, this global energetic information on the discretization error may be useless for practical application. Therefore, goal-oriented error aims at providing bounds on the error on a quantity of interest. In this paper, the quantity of interest is considered to be a linear form of the displacement field, defined by extractors [3].
### 3.3.2 Definition of the linear quantity of interest and the adjoint problem

In reliability, the quantity of interest is defined as a function of performance defining the domain of failure and of safety. This performance function is also called limit state function and is written as a margin between a resistance R and a solicitation S. The resistance is either deterministic or random. This limit state function may be written:

$$g = R - S \qquad (9)$$

The sign of this function defines the failure or not of the structure. To estimate bounds on the quantity g, the solicitation S needs to be defined by a linear functional $\tilde{L}$:

$$g = R - \tilde{L}(\underline{u}) = R - \int_{\Omega} (\underline{\sigma}_{\Sigma} : \underline{\varepsilon}(\underline{u}) + \underline{f}_{\Sigma}\underline{u}) d\Omega \qquad (10)$$

where $\underline{\sigma}_{\Sigma}$ and $\underline{f}_{\Sigma}$ are extractors. It is possible to treat non linear quantities of interest. For some of them, there are specific methods to calculate guaranteed bounds [49, 44]. Otherwise, this quantity would need to be linearized. It would result in the loss of the guarantee that $g_{ex}$ lies between estimated bounds. Finally, the probability of failure that will later be estimated is likely to be biased. We define the subspace of statically admissible fields for the adjoint problem:

$$\tilde{\mathcal{S}}\tilde{\mathcal{A}}(\omega) = \left\{ \underline{\tau} \in (L^2(\omega))_{\text{sym}}^{d \times d}; \forall \underline{v} \in \text{KA}^{00}(\omega), \int_{\omega} \underline{\tau} : \underline{\varepsilon}(\underline{v}) d\omega = \tilde{L}(\underline{v}) \right\} \qquad (11)$$

The adjoint problem set on $\Omega$ reads:

$$\text{Find } (\tilde{\underline{u}}_{ex}, \tilde{\underline{\sigma}}_{ex}) \in \text{KA}^0(\Omega) \times \tilde{\mathcal{S}}\tilde{\mathcal{A}}(\Omega) \text{ such that } e_{CR_{\Omega}}(\tilde{\underline{u}}_{ex}, \tilde{\underline{\sigma}}_{ex}) = 0 \qquad (12)$$

The solution to this problem exists and is unique.

The adjoint problem is usually solved using the finite element method. The mesh does not need to be the same as for the forward problem. The adjoint displacement field obtained by solving the adjoint problem with the FEM is $\tilde{\underline{u}}_{\tilde{H}}$.

We define the discretization error of the adjoint problem as:

$$e_{discr} = \|\tilde{\underline{\varepsilon}}_{discr}\|_{\Omega} = \|\tilde{\underline{u}}_{ex} - \tilde{\underline{u}}_{\tilde{H}}\|_{\Omega} \qquad (13)$$

### 3.3.3 Error estimation on a quantity of interest

We denote $g_{ex} = R - \tilde{L}(\underline{u}_{ex})$ the unknown exact value of the limit state function. $g_H = R - \tilde{L}(\underline{u}_H)$ is an approximation of this quantity of interest.
Bounds from [25, 26] can be applied to the quantity of interest $g_{ex}$:

$$|g_{ex} - g_H + I_{HH}| \le \frac{1}{2} e_{CR_{\Omega}}(\underline{u}_H, \underline{\hat{\sigma}}_H) e_{CR_{\Omega}}(\tilde{u}_{\tilde{H}}, \tilde{\underline{\hat{\sigma}}}_{\tilde{H}}) \quad (14)$$

where

$$I_{HH} = \int_{\Omega} \frac{1}{2} (\underline{\hat{\sigma}}_{\tilde{H}} + \mathbb{H} : \underline{\varepsilon}(\tilde{u}_{\tilde{H}})) : \mathbb{H}^{-1} : (\underline{\hat{\sigma}}_H - \mathbb{H} : \underline{\varepsilon}(\underline{u}_H)) d\Omega \quad (15)$$

and where $\underline{\hat{\sigma}}_{\tilde{H}} \in \tilde{S}A_H(\Omega)$.

Therefore, the error on the quantity of interest can be obtained from the product of the two errors in energy norm of both reference and adjoint problem. In this paper, we choose to solve both reference and adjoint problems on the same mesh so that it does not require an additional factorization of the stiffness matrix. Solving the finite element problem is in fact simplified to a multiple (double) right-hand side linear equation.

It is to be noted that the bounds provided here do not guarantee that $g_H$ is inside them. In fact, for the structure studied in 5, the finite element solution is outside the bounds. It would be absurd to use a value that is guaranteed to be false. Therefore, it was decided to use the middle of the bounds as the output of the finite element code :

$$g^m = g_H - I_{HH} \quad (16)$$

In fact, it was observed in [31] that $g^m$ is in most cases a better approximation of the true value $g_{ex}$. A bounding of $g_{ex}$ is therefore :

$$g^{-} \le g_{ex} \le g^{+} \quad (17)$$

where

$$g^{-} := g^{m} - \frac{1}{2} e_{CR_{\Omega}}(\underline{u}_{H}, \underline{\hat{\sigma}}_{H}) e_{CR_{\Omega}}(\tilde{u}_{\tilde{H}}, \tilde{\underline{\hat{\sigma}}}_{\tilde{H}}) \quad (18)$$

and

$$g^{+} := g^{m} + \frac{1}{2} e_{CR_{\Omega}}(\underline{u}_{H}, \underline{\hat{\sigma}}_{H}) e_{CR_{\Omega}}(\tilde{u}_{\tilde{H}}, \tilde{\underline{\hat{\sigma}}}_{\tilde{H}}) \quad (19)$$

## 3.4 Estimation of the probability of failure

For the sake of simplicity, in the previous subsection, the mechanical problem is presented as deterministic. In reality, the geometry, the loading or the behavior of the structure are random and as a consequence the displacement fields $\underline{u}$ and $\underline{u}_H$ are random fields.

We gather the $n$ random variables of the mechanical problem into a random vector $\{X\}$ with joint probability density $f_{\{X\}}(\{x\})$.
The exact probability of failure $P_{f,ex}$ is defined from the exact displacement solution $u_{ex}$ via the exact limit state function $g_{ex}$:

$$P_{f,ex} = \int_{g_{ex}(\{X\}) \le 0} f_{\{X\}}(\{x\})dx_1 \dots dx_n \quad (20)$$

The probability of failure computed using the finite element solution $u_H$ is found by replacing the limit state function $g_{ex}$ by $g_H$. As a consequence, the discretization error leads to an error on the estimation of the probability of failure since $P_{f,H} \neq P_{f,ex}$.

In addition, the form of the limit-state function is usually unknown and surrogates techniques aim at building an approximation $\hat{g}$ of $g$ using samples. In this paper, we will consider that $\hat{g}$ is built using kriging techniques. One can define the probability of failure computed using the surrogate model $\hat{g}$ built from the exact mechanical solutions :

$$\hat{P}_{f,ex} = \int_{\hat{g}_{ex}(\{X\}) \le 0} f_{\{X\}}(\{x\})dx_1 \dots dx_n \quad (21)$$

Therefore, the meta-modeling error leads to an error in the estimation of the probability of failure as $\hat{P}_{f,ex} \neq P_{f,ex}$.

Since the exact solution is almost never known, the samples used to construct the surrogate model are computed with finite element solutions so that the probability of failure is

$$\hat{P}_{f,H} = \int_{\hat{g}_H(\{X\}) \le 0} f_{\{X\}}(\{x\})dx_1 \dots dx_n \quad (22)$$

The only difference between integrals (20), (21) and (22) is the domain of integration that is defined by the limit state $g = 0$. Therefore, the assessment of the probability of failure can be seen as a classification problem for which only the limit state needs to be precisely characterized.

# 4 Multi-fidelity kriging for reliability analysis

In this section, we first explain the principle of kriging and its application to the estimation of failure probabilities $\hat{P}_{f,ex}$ and $\hat{P}_{f,H}$. Then, we propose a methodology to construct a kriging-based multifidelity surrogate model exploiting the discretization error bounds. The objective will be to classify evaluations of $g$ correctly ($g \le 0$ means the sample is failed and $g > 0$ means safe).
## 4.1 Kriging

Metamodeling consists in computing an estimation of a quantity of interest cheaper to evaluate than the initial function g, where g can either be g<sub>ex</sub> or g<sub>H</sub>.

Let consider that m observations (or computations) at x̄<sub>1</sub><sup>obs</sup> ... x̄<sub>m</sub><sup>obs</sup> of g are gathered in a vector of observations g<sup>obs</sup>.

Kriging [24] relies on the hypothesis that the objective function g is the realization of a stationary Gaussian process G

$$G(\underline{x}) = y(\underline{x}) + Z(\underline{x}) \qquad (23)$$

where x̄ ∈ ℝ<sup>n</sup>, y(x̄) is the mean of G(x̄) and Z(x̄) is a random variable distributed as a centered Gaussian. Three types of kriging derive from the assumption made on the form of y. In simple kriging, y(x̄) = β<sub>0</sub> where β<sub>0</sub> is the mean computed from the observations g<sup>obs</sup>. In ordinary kriging, y(x) = β where β is an unknown constant. In universal kriging, y(x̄) = f(x̄)<sup>T</sup>β where f(x̄) = [f<sub>1</sub>(x̄), ..., f<sub>k</sub>(x̄)] is a vector containing k basis functions generally chosen as polynomials. To remain general, universal kriging will be considered in this paragraph.

The estimation ˆg of g is chosen as a linear combination of observations:

$$\hat{g}(\underline{x}) = \underline{w}(\underline{x})^T \underline{g}^{\text{obs}} \qquad (24)$$

where w(x̄) are the kriging coefficients to be determined.

ˆg is searched as the best linear unbiased predictor, which reads :

$$\begin{cases} \mathbb{E}[\hat{g} - G] = 0 \\ \beta = \arg \min \mathbb{E}[(\hat{g} - G)^2] \end{cases} \qquad (25)$$

where E is the expectation.

Solving system (25) using (24) and (23) enables to determine the kriging coefficients. However, it requires the estimation of the covariance matrix of the data E(z̄<sub>obs</sub>z̄<sub>obs</sub><sup>T</sup>), where z̄<sub>obs</sub> = ḡ<sub>obs</sub> - [y(x̄<sub>1</sub><sup>obs</sup>), ..., y(x̄<sub>m</sub><sup>obs</sup>)]. Generally, the form of the correlation function is chosen a priori. In this paper, we use a Gaussian correlation function as it was observed in [29] to give a better approximation than other standard correlation functions for several test cases. Therefore:

$$\mathbb{E}(z_{obs} \tilde{z}_{obs}^T) = \sigma^2 \mathbf{R} = \sigma^2 \left[ \prod_{k=1}^{n} e^{-\left( \frac{x_k^{obs,i} - x_k^{obs,j}}{\theta_k} \right)^2} \right]_{i,j \in \{1, m\}} \qquad (26)$$

where σ<sup>2</sup> and θ = [θ<sub>1</sub>, ..., θ<sub>n</sub>] are respectively the variance of the process and the vector containing the fluctuation parameters. The fluctuation parameter, or correlation length, can be interpreted as a smoothness parameter: the higher
its value, the farther two samples are correlated. These parameters are called hyperparameters.

Several methods exist to estimate $\sigma^2$ and $\theta$: the variograms [32], the leave one-out method [43] or the maximum likelihood estimation [46]. In this paper, the maximum likelihood estimation is considered as it is asymptotically optimal [52, p. 124]. The likelihood that the data is an outcome of a Gaussian process can be written:

$$ \mathcal{L}(\sigma^2, \underline{\theta}) = \frac{1}{\sqrt{(2\pi\sigma^2)^m \det \mathbf{R}}} e^{-\frac{\underline{z}_{obs} \mathbf{R}^{-1} \underline{z}_{obs}}{2\sigma^2}} \quad (27) $$

where $\det \mathbf{R}$ is the determinant of the correlation matrix. Then $\hat{\theta}$ is determined as

$$ \hat{\underline{\theta}} = \operatorname*{argmin}_{\underline{\theta}} (-\log(\mathcal{L}(\sigma^2, \underline{\theta}))) \quad (28) $$

Once $\hat{\underline{\theta}}$ is known, the variance $\sigma^2$ is computed by calculating stationary points of $\mathcal{L}$, which leads to:

$$ \sigma^2 = \frac{\underline{z}_{obs} \mathbf{R}^{-1} \underline{z}_{obs}}{m} \quad (29) $$

Once the hyperparameters are determined, the optimization system (25) can be solved using the Lagrangian method.

In addition to the estimator $\hat{g}$, kriging meta-modeling provides an estimation of the uncertainty on $\hat{g}$ with the variance $s_{\hat{g}}^2$.

## 4.2 Adaptive kriging

Using a kriging metamodel allows evaluating the uncertainty on the output through $s_{\hat{g}}^2$. Therefore, the meta-model can be improved by adding observations on points corresponding to large uncertainty on $\hat{g}$. This strategy is developed for reliability analysis in [13] where the learning function chosen for enrichment is:

$$ \underline{x} \rightarrow U(\underline{x}) = \frac{|\hat{g}(\underline{x})|}{s_{\hat{g}}^2(\underline{x})} \quad (30) $$

Low $U$ value corresponds to samples close to the limit state ($\hat{g} = 0$) with high uncertainty (high $s_{\hat{g}}^2$ value). Therefore, the limit state function $g$ can be computed at the minimizer of the learning function and added to the observations, which results in the construction of a better meta-model. Then, the meta-model is used to estimate the probability of failure with Monte Carlo estimation [34]:

$$ \hat{P}_{f,H} \approx \frac{1}{n_{MC}} \sum_{i=1}^{n_{MC}} ind_{\hat{g}_H \le 0}(X_i) \quad (31) $$

Where $ind_{\hat{g}_H \le 0}$ is the indicator function of the failure domain $\hat{g}_H \le 0$.

Further details on the algorithm can be found in [13].
4.3 Multi-fidelity kriging

We propose a methodology to construct a multi-fidelity kriging-based meta-model using computations from different meshes corresponding to different levels of fidelity. The refinement strategy is based on guaranteeing the classification of samples used to build the surrogate. For a realization $\underline{x}_i$ of the random variables, the finite element code provides the values $g_H(\underline{x}_i) = R - \tilde{L}(\underline{u}_H)$ and the error estimation enables to obtain $g_m(\underline{x}_i) = R - \tilde{L}(\underline{u}_H) - I_{hh}$ and bounds on $g_{ex}(\underline{x}_i)$ using (17).

If the bounds have the same sign, the state of the point $\underline{x}_i$ (safe or failure) is guaranteed. If not, it means that the discretization error may lead to a wrong classification of the point. Thus, a finer mesh can be defined and used for a new finite element solving for the same random variable realization.
The algorithm we propose is given in Algorithm 1.

<caption>**Algorithm 1: MF-AK**</caption>

Choose a coarsest mesh size $h_{max}$ and a finest mesh size $h_{min}$;
Generate Monte Carlo (MC) population;
Generate $n_1$ samples for initialization;
**for** $i = 1..n_1$ **do**
  Evaluate $g_H$, $g_-$ and $g_+$ on $\underline{x}_i$;
  **while** $g_-g_+ < 0$ and $h_{min} \le h$ **do**
    Perform the remeshing strategy;
    Evaluate $g_H$, $g_-$ and $g_+$ on $\underline{x}_i$;
  **end**
**end**
Construct the meta-model $\hat{g}$ and evaluate it on the MC population;
Compute the learning function $U$;
Compute the failure probability with the meta-model $\hat{g}$;
Compute the coefficient of variation COV of the failure probability;
**while** $\min(U) \le 2$ or $\text{COV} > \zeta$ **do**
  **while** $\min(U) \le 2$ **do**
    Select the learning point $\underline{x}_{new} = \text{argmin}(U)$;
    Evaluate $g_H$, $g_-$ and $g_+$ on $\underline{x}_{new}$;
    **while** $g_-g_+ < 0$ and $h_{min} \le h$ **do**
      Perform the remeshing strategy;
      Evaluate $g_H$, $g_-$ and $g_+$ on $\underline{x}_{new}$;
    **end**
    Construct the meta-model $\hat{g}$ and evaluate it on the MC population;
    Compute the learning function $U$;
  **end**
  Compute the failure probability with the meta-model $\hat{g}$;
  Compute the coefficient of variation COV of the failure probability;
  **if** $\text{COV} > \zeta$ **then**
    Enrichment of the Monte Carlo population;
    Evaluate the meta-model on the population;
    Compute the learning function $U$;
  **end**
**end**

This algorithm requires two new parameters compared to classical AK-MCS : the size of the coarsest mesh and the size of the finest mesh. Moreover, a remeshing strategy has to be defined. Several methods can be proposed. The first strategy is to create a family of meshes ranked from coarsest to finest prior to the reliability
analysis. A second strategy would be to exploit the error maps provided by the error estimation procedure to generate adaptive mesh. Yet, for quantity of interest, this procedure is not straightforward, as illustrated in [10]. If the convergence rate with respect to the mesh size $\alpha$ is known, a third strategy can be defined. From an error estimation on a mesh of size $l_c$ and by doing the hypothesis that reducing the mesh size reduces the interval at a convergence rate $l_c^{2\alpha}$ but does not translate the error interval, then, the optimal mesh size $l_{opt}$ can be computed as

$$l_{c, opt} = l_c \left( \frac{2|g_m|}{e_{CR\Omega}(\underline{u}_H, \underline{\hat{\sigma}}_H)e_{CR\Omega}(\tilde{u}_{\tilde{H}}, \tilde{\hat{\sigma}}_{\tilde{H}})} \right)^{\frac{1}{2\alpha}} \quad (32)$$

The meta-model is updated using the learning function defined in (30), as it is done in AK-MCS. Once the meta-model is built, the probability of failure is computed through Monte Carlo sampling.

## 4.4 Bounds on the probability of failure

For several mechanical problems, the limit state function is known a priori as monotonic against random variables, which enables to compute bounds on the probability of failure as presented in [9]. Let assume that there is only one random variable $x$ and that the limit state function decreases against $x$ (for example see Figure 1). If the bounds calculated for a given realization $x^-$ guarantee that it is in the failure domain (respectively in the safe domain for $x^+$), then, any sample with a greater value (respectively lower value) of the realization of the random variable is known to be in the failure domain (respectively safe). For Figure 1, these two samples would be $x^- = x_1$ and $x^+ = x_4$. Indeed, even on the finest mesh, the state of $x_2$ and $x_3$ cannot be guaranteed. Bounds on $P_{f,ex}$ can then be computed using the cumulative probability function $F$ associated to the joint probability density function $f$ of the random variable $x$:

$$F(x^-) > P_{f,ex} > F(x^+) \quad (33)$$

# 5 Numerical assessment

## 5.1 Description of the structure

We consider a cracked plate of width $w = 7m$ and length $L = 16m$. The length of the crack tip is $a$, as illustrated on Figure 2. We suppose a linear elastic homogeneous isotropic behavior characterized by the Hooke tensor $\mathbf{H}$. We assume small perturbations, plane strains. Therefore, the mechanical problem is 2D ($d =$
<figure><img src="image_2.png" /><figcaption>Figure 1: Exact limit state and multifidelity metamodel built with two level of fidelity $g_H$ and $g_h$ (respectively mesh size 0.5m and 0.05m - Number indicates sample index</figcaption></figure>

2). $E = 210$ GPa is the Young's modulus and $\nu = 0.3$ is the Poisson ratio. The structure undergoes traction with $P = 1Pa$, which corresponds to mode I opening.

The stress intensity factor $K_I$ is computed using auxiliary solutions and integral over a crown (represented in orange in Figure 2) [48, 17]. $K_I$ is therefore a linear functional of the displacement $\underline{u}$ and this linear functional defines the loading of the adjoint problem. The Griffith criterion [21] is used as the limit state function $G = K_{lim} - K_I$ where $K_{lim}$ is the critical value. For this structure, empirical expressions of the stress intensity factor $K_I$ are provided in [50]. The following is guaranteed with at most 0.5% error:

$$K_{I,ex} = \left( \sqrt{\frac{2w}{\pi a} \tan \frac{\pi a}{2w}} \frac{0.752 + 2.02 \frac{a}{w} + 0.37 \left(1 - \sin \frac{\pi a}{2w}\right)^3}{\cos \frac{\pi a}{2w}} \right) P \sqrt{\pi a} \quad (34)$$

This expression will be used to compute the reference exact solution of the problem $g_{ex}$.
<figure><img src="image_3.png" /><figcaption>Figure 2: Cracked plate with loading of the reference problem (blue) and of the adjoint problem (orange)</figcaption></figure>

The length $a$ is modeled as a random variable with beta distribution bounded between 1.6m and 4.25m and of shape parameters both equal to 2.

## 5.2 Illustration of the influence of the discretization error on estimation of the probability of failure

Considering $K_{lim}$ deterministic, the exact probability of failure is computed from (34) knowing the cumulative distribution function of $a$. The lack of precision of formula (34) (0.5%) is propagated on the probability of failure. These two probabilities obtained are referred as mesh size $0^-$ and $0^+$ in Table 2. The middle of these bounds is $0^m$. We note $\epsilon = \frac{P_f - P_f(0^m)}{P_f(0^m)}$ the relative error.

The algorithm AK-MCS is used to estimate the failure probability for different mesh sizes $l$ and for two values of $K_{lim}: K_{lim} = 14Pa\sqrt{m}$ and $K_{lim} = 9Pa\sqrt{m}$ for considering two order of magnitude of $P_f$ namely $10^{-3}$ and $10^{-1}$. The parameters of AK-MCS are depicted in Table 1. All the metamodels are built using the DACE toolbox [29]. The same Monte Carlo population is used for every simulation.

The results in terms of probability of failure are gathered in Table 2.

Firstly, we observe that the probability of failure depends on the discretization size $l$. Secondly, as a too coarse mesh leads to an overestimation of the stiffness of the structure, the probability of failure is underestimated, which may be dramatic in the context of uncertainty qualification through virtual testing. Finally, to obtain less than 40% error on the probability of failure, a mesh size of $l = 0.5m$
<table><caption>Table 1: AK-MCS parameters</caption><tr><td>Stopping criterion on learning</td><td>2</td></tr><tr><td>Coefficient of Variation criterion ζ</td><td>10<sup>-2</sup></td></tr><tr><td>Size of initial Monte-Carlo population</td><td>1000</td></tr><tr><td>Type of initialization</td><td>Factorial experiment</td></tr><tr><td>Initial correlation length</td><td>0.55m</td></tr><tr><td>Number of samples to build initial metamodel</td><td>5</td></tr><tr><td>Maximal factor to multiply P<sub>MC</sub> between two iterations</td><td>5</td></tr></table>

<table><caption>Table 2: Probabilities of failure and relative error for different mesh sizes and resistance values</caption><thead><tr><th></th><th colspan="2">K<sub>lim</sub> = 9Pa√m</th><th colspan="2">K<sub>lim</sub> = 14Pa√m</th></tr><tr><th><i>l</i></th><th><i>P</i><sub>f</sub></th><th>ε</th><th><i>P</i><sub>f</sub></th><th>ε</th></tr></thead><tbody><tr><td>0<sup>-</sup></td><td>2.23 10<sup>-1</sup></td><td>0.02</td><td>5.27 10<sup>-3</sup></td><td>0.12</td></tr><tr><td>0<sup>m</sup></td><td>2.27 10<sup>-1</sup></td><td>0</td><td>5.98 10<sup>-3</sup></td><td>0</td></tr><tr><td>0<sup>+</sup></td><td>2.31 10<sup>-1</sup></td><td>0.02</td><td>6.72 10<sup>-3</sup></td><td>0.12</td></tr><tr><td>0.5m</td><td>1.40 10<sup>-1</sup></td><td>0.38</td><td>0</td><td>1</td></tr><tr><td>0.4m</td><td>1.71 10<sup>-1</sup></td><td>0.25</td><td>0</td><td>1</td></tr><tr><td>0.3m</td><td>1.82 10<sup>-1</sup></td><td>0.20</td><td>3.79 10<sup>-5</sup></td><td>0.99</td></tr><tr><td>0.2m</td><td>1.94 10<sup>-1</sup></td><td>0.14</td><td>1.10 10<sup>-4</sup></td><td>0.90</td></tr><tr><td>0.1m</td><td>2.11 10<sup>-1</sup></td><td>0.07</td><td>2.46 10<sup>-3</sup></td><td>0.59</td></tr><tr><td>0.05m</td><td>2.19 10<sup>-1</sup></td><td>0.03</td><td>3.91 10<sup>-3</sup></td><td>0.35</td></tr><tr><td>0.02m</td><td>2.24 10<sup>-1</sup></td><td>0.01</td><td>4.93 10<sup>-3</sup></td><td>0.17</td></tr></tbody></table>

would be sufficient for the case $K_{lim} = 9Pa\sqrt{m}$ while a mesh size $l = 0.05m$ would be needed for the case $K_{lim} = 14Pa\sqrt{m}$. This can be explained by two reasons: the discretization error does not depend linearly on the crack length $a$ and, for comparable precisions, the distribution of $a$ has an influence on the obtained value of failure probability. The latter is depicted by the fact that an error of 0.5% with the empirical formula on $K_I$ leads to a relative error of 12% on the failure probability for the case $K_{lim} = 14Pa\sqrt{m}$ and only 2% for the case $K_{lim} = 9Pa\sqrt{m}$. Therefore, as the conception point is unknown, doing an *a priori* choice of the discretization size is impossible. Moreover, using the same mesh for different safety studies may still lead to errors.

## 5.3 Multi-fidelity kriging

In this subsection, we apply the new Algorithm 1 with the same parameters as for AK-MCS (see Table 1). The size of the coarsest mesh is 0.5m and the size of
in Table 2. The reason is that the metamodel is built with the values $g_m$. It has been observed that $g_m$ is usually a better approximation than $g_H$.

Moreover, whatever the intermediate mesh sizes and population, the method is faster than AK-MCS on the finest mesh using $g_m$ (Simulation 9).

From Table 3, we can observe that on average, increasing the number of possible meshes improves the repartition of simulations between the different levels of fidelity. However, it also increases the average CPU time as some evaluations may not be useful to determine the sign of the limit state function. This is highlighted by simulation 5 on population 1 (S5P1) where only necessary levels of fidelity according to S1P1 are considered. The difference in CPU time is also attributable to a greater number of samples being evaluated during enrichment. It is to be noticed that S2P1 is shorter than simulation S3P1 and S4P1. The reason is a smaller number of calls to the finest mesh.

Figure 3 shows the limit state function against the crack length. It shows that samples close to the limit state are calculated on the finest mesh, which is paramount when building a multifidelity metamodel for reliability purposes. It is to be noted that every sample added during enrichment phase is calculated on the finest mesh.

As highlighted, it is not trivial to choose more than two levels of fidelity as the best repartition depends on the quality of bounds, the limit state function and samples for which the finite element solver will be called. It is therefore suggested to only use the coarsest and finest mesh.

This strategy is compared to the one exploiting the computation of an optimal mesh size $l_{opt}$. The results are presented in Table 4. This table is presented in a different way than Table 3. Each column corresponds to one sample added to the metamodel and each row within each column corresponds to the intermediate mesh sizes used for this sample. For example, for the strategy with optimal mesh size on population 1, the four first samples used to initialize the metamodel are guaranteed by the bounds on the coarsest mesh. However, the fifth sample is not guaranteed and is therefore calculated on a finer mesh (i.e. 0.21m) and its state is then guaranteed. During the enrichment phase, two samples were added to the metamodel and were both calculated on mesh size 0.059m and 0.05m (finest mesh). The second part of Table 4 presents simulation 8 from Table 3 and serves as a reference for comparaison with the optimal mesh size estimation strategy.

First, it is to be noted that whatever the Monte Carlo population, the optimal mesh size estimation strategy converges to the same probability of failure as the reference simulation. For most samples, the first optimal mesh size estimation allows to guarantee its state and the strategy finds a coarser mesh than the finest one. The finest mesh is always reached allowing to estimate the same probability of failure as the reference simulation. The finding of a good mesh compromise
<figure><img src="image_4.png" /><figcaption>Figure 3: Exact limit state and multifidelity metamodel built with two level of fidelity $g_H$ and $g_h$ (respectively mesh size 0.5m and 0.05m)</figcaption></figure>

allows to reduce CPU time on all populations.

In most engineering problems the mesh size convergence rate is not known a priori. It is therefore not possible to use the optimal mesh size estimation strategy to compute the probability of failure. In this case, it is recommended to only choose a coarse mesh and the fine one.

### 5.3.2 Comparison with a meta-model built by evofusion

In [15], authors propose an approach named evofusion to construct a meta-model from data for multiple levels of fidelity (at least two). Let $g_2$ and $g_1$ denote the functions estimated respectively from a high and low fidelity model. $g_2$ may be written as:

$$g_2 = g_1 + \underbrace{g_2 - g_1}_{g_{err}} \qquad (35)$$

The principle consists in creating a first kriging metamodel for $g_1$ and a second
<table><caption>Table 4: Comparison of strategy between two levels of fidelity (same as Set 8 in Table 3) and optimal mesh size estimation</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Pop.</th><th colspan="5">Mesh size (m)</th><th rowspan="2">P<sub>f</sub><br>(×10<sup>-3</sup>)</th><th rowspan="2">CPU time<br>(×10<sup>3</sup>s)</th></tr><tr><th>0.5 (7)</th><th>0.21 (1)</th><th>0.095 (1)</th><th>0.059 (1)</th><th>0.05 (1)</th></tr></thead><tbody><tr><td rowspan="5">With<br>optimal<br>mesh size</td><td>1</td><td>0.5 (7)</td><td>0.21 (1)</td><td>0.095 (1)</td><td>0.059 (1)</td><td>0.05 (1)</td><td>5.1</td><td>4.1</td></tr><tr><td>2</td><td>0.5 (7)</td><td>0.21 (1)</td><td>0.094 (1)</td><td>0.062 (1)</td><td>0.05 (1)</td><td>5.1</td><td>4.3</td></tr><tr><td>3</td><td>0.5 (7)</td><td>0.21 (1)</td><td>0.094 (1)</td><td>0.062 (1)</td><td>0.05 (1)</td><td>5.1</td><td>4.2</td></tr><tr><td>4</td><td>0.5 (7)</td><td>0.21 (1)</td><td>0.15 (1)</td><td>0.069 (1)</td><td>0.05 (2)</td><td>5.1</td><td>5.9</td></tr><tr><td>5</td><td>0.5 (7)</td><td>0.21 (1)</td><td>0.098 (1)</td><td>0.055 (1)</td><td>0.05 (1)</td><td>5.1</td><td>4.0</td></tr><tr><td rowspan="5">Without<br>optimal<br>mesh size</td><td>1</td><td>0.5 (7)</td><td>0.05 (3)</td><td></td><td></td><td></td><td>5.1</td><td>5.8</td></tr><tr><td>2</td><td>0.5 (7)</td><td>0.05 (3)</td><td></td><td></td><td></td><td>5.1</td><td>5.6</td></tr><tr><td>3</td><td>0.5 (7)</td><td>0.05 (3)</td><td></td><td></td><td></td><td>5.1</td><td>5.6</td></tr><tr><td>4</td><td>0.5 (7)</td><td>0.05 (2)</td><td></td><td></td><td></td><td>5.1</td><td>3.8</td></tr><tr><td>5</td><td>0.5 (7)</td><td>0.05 (3)</td><td></td><td></td><td></td><td>5.1</td><td>5.7</td></tr></tbody></table>

for $g_{err}$. The sum of the two metamodels is an evofused metamodel that is built from 2 levels of fidelity. The advantage of using such a metamodel is that the global trend may be captured thanks to low fidelity data. The error metamodel between level 1 and 2 may then be built thanks to a few calls to both $g_1$ and $g_2$. A measure of the uncertainty on the evofused meta-model is computed from the standard deviations of the low-fidelity metamodel and the error metamodel by considering that $g_{err}$ and $g_1$ are uncorrelated, which gives:

$$s_{\hat{g}_2}^2 = s_{\hat{g}_1}^2 + s_{g_{err}}^2 \qquad (36)$$

For readers familiar with co-kriging [23], equation (35) is very similar to the co-kriging formulation:

$$g_2 = \rho_1 \times g_1 + \delta_1 \qquad (37)$$

Where $\rho_1$ and $\delta_1$ are unknowns.

Actually, evofusion might be seen as the particular case of co-kriging where $\rho_1 = 1$ and $\delta_1 = g_{err}$.

In this paper, we consider an adaptation of evofusion to AK-MCS. For this purpose, a learning function needs to be introduced. The U function is adapted:

$$x \rightarrow U(x) = \frac{|\hat{g}_1(x) + \hat{g}_{err}(x)|}{s_{\hat{g}_1}^2(x) + s_{\hat{g}_{err}}^2(x)} \qquad (38)$$

Given $x_i$ with $U(x_i) < 2$, a refinement strategy is to use the decomposition of $s_{\hat{g}}^2$ as a sum for every level of fidelity. $g_1(x_i)$ could be simulated and added to $\hat{g}_1$ thus
reducing $s_{\tilde{g}}^2(x_i)$ (as $s_{\tilde{g}_1}^2(x_i) = 0$). If $s_{\tilde{g}}^2(x_i)$ is still smaller than 2 then a call to the high-fidelity model is necessary.

The low fidelity model is built from the finite element solution for a mesh size of 0.5m (no error estimation is performed) and the high fidelity model is built from $g^m$ (obtained after error estimation) on a mesh size of 0.05m : the same high fidelity is chosen as the finest mesh from 5.3.1 to be able to compare results.

To build an initial metamodel, 5 samples are selected with a factorial experiment over the whole space on which $x^i$ is defined. $g_1$ is called for each of them. Three samples are selected among them to calculate $g_2$ and be able to compute $g_{err}$. These samples are selected by picking the three values closest to the limit state $g_1 = 0$.

Parameters of the evofusion metamodel are shown in Table 5.

<table><caption>Table 5: Evofusion parameters</caption><tr><td>Stopping criterion on learning</td><td>2</td></tr><tr><td>Coefficient of Variation criterion ζ</td><td>10<sup>-2</sup></td></tr><tr><td>Size of initial Monte-Carlo population</td><td>1000</td></tr><tr><td>Type of initialization</td><td>Factorial experiment</td></tr><tr><td>Initial correlation length for g<sub>1</sub> and for g<sub>er</sub></td><td>0.55m</td></tr><tr><td>Number of sample to build initial metamodel</td><td>5</td></tr><tr><td>Maximal factor to multiply P<sub>MC</sub> between two iterations</td><td>5</td></tr></table>

**Results** The obtained probabilities of failure and CPU time for both evofusion and multi-fidelity kriging are given in Table 6.

As shown in Table 6, evofusion does not reduce the number of evaluations needed to estimate the probability of failure. Actually, in simulations for which the method converges, the number of FEM evaluations is still higher than the reference simulation of multifidelity kriging.

For some Monte-Carlo populations, the method using evofusion does not converge. The metamodel and the samples used to build it were inspected. When too many samples are added to the metamodel, it becomes likely that two of them have very close input values. In order to fit those two samples, the metamodel needs to bend sharply thus making the correlation length drop. Finally, the metamodel uncertainty increases and the number of samples needed to improve its quality plumets: the method never converges. It is possible that the trigger of this problem is the learning function $U$. Actually, it was seen to call for too many enrichment in [18].
<table><caption>Table 6: Comparison between evofusion and multi-fidelity kriging</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">Mesh size (m)</th><th rowspan="2">P<sub>f</sub></th><th rowspan="2">CPU time</th></tr><tr><th>Pop.</th><th>0.5</th><th>0.05</th></tr></thead><tbody><tr><td rowspan="5">Kriging</td><td>1</td><td>7</td><td>3</td><td rowspan="5">5.1 for all</td><td>6</td></tr><tr><td>2</td><td>7</td><td>3</td><td>6</td></tr><tr><td>3</td><td>7</td><td>3</td><td>6</td></tr><tr><td>4</td><td>6</td><td>2</td><td>4</td></tr><tr><td>5</td><td>7</td><td>3</td><td>6</td></tr><tr><td rowspan="5">Evofusion</td><td>1</td><td>28</td><td>5</td><td>5.1</td><td>6</td></tr><tr><td>2</td><td>128</td><td>5</td><td>5.1</td><td>11</td></tr><tr><td>3</td><td colspan="3">Error<br>metamodel<br>not converged</td></tr><tr><td>4</td><td colspan="3">Low fidelity metamodel<br>not converged</td></tr><tr><td>5</td><td>97</td><td>6</td><td>5.1</td><td>13</td></tr></tbody></table>

### 5.3.3 Bounds on the probability of failure

Bounds on $P_f$ can be computed using samples that are guaranteed to be in the failure domain as highlighted in 4.4. Table 7 shows bounds added to simulations from Table 3. Bounds can be thin, in the same range as the one provided in Table 2 with the empirical formula. However, for some cases there is no guaranteed sample close to the limit state. Nothing guarantees a priori that bounds of good quality will be computed with this method. The best practice would be to add with expertise a few samples close enough to the limit state and with guaranteed sign. In any case, this post-process is extremely cheap and offer an estimation of safety margins.

## 6 Conclusion

In this article, we presented the construction of a multi-fidelity kriging-based meta-model for the estimation of the probability of failure. By exploiting discretization error estimators, it is possible to ensure the state (safe or failure) of the points used to build the meta-model. Therefore, the correct classification of those points is guaranteed. It allows to define a strategy to build the meta-model from computations on different mesh sizes thus adapting the discretization to the objective. Such a strategy allows to reduce total computational cost compared to a strategy using a unique mesh size and to focus expensive computations on critical points. Results using a kriging metamodel were compared to the ones using evofusion, a more
<table><caption>Table 7: Bounds on the probability of failure depending on the refinement strategy and Monte Carlo population (P1 to P5)</caption><thead><tr><th rowspan="2">Set</th><th colspan="5">P<sub>f</sub><sup>-</sup> (×10<sup>-3</sup>)</th><th colspan="5">P<sub>f</sub><sup>+</sup> (×10<sup>-2</sup>)</th></tr><tr><th>P1</th><th>P2</th><th>P3</th><th>P4</th><th>P5</th><th>P1</th><th>P2</th><th>P3</th><th>P4</th><th>P5</th></tr></thead><tbody><tr><td>1</td><td>3.8</td><td>0</td><td>0</td><td>0</td><td>0</td><td>5.7</td><td>16</td><td>16</td><td>16</td><td>16</td></tr><tr><td>2</td><td>3.8</td><td>0</td><td>0</td><td>0</td><td>4.0</td><td>5.7</td><td>16</td><td>16</td><td>16</td><td>16</td></tr><tr><td>3</td><td>3.8</td><td>0</td><td>0</td><td>0</td><td>0</td><td>16</td><td>16</td><td>16</td><td>16</td><td>16</td></tr><tr><td>4</td><td>4.1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>16</td><td>16</td><td>16</td><td>16</td><td>16</td></tr><tr><td>5</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>16</td><td>16</td><td>16</td><td>16</td><td>16</td></tr><tr><td>6</td><td>4.1</td><td>0</td><td>0</td><td>0</td><td>4.0</td><td>5.7</td><td>16</td><td>16</td><td>16</td><td>16</td></tr><tr><td>7</td><td>3.8</td><td>0</td><td>0</td><td>0</td><td>0</td><td>5.7</td><td>16</td><td>16</td><td>16</td><td>16</td></tr><tr><td>8</td><td>3.8</td><td>0</td><td>0</td><td>0</td><td>0</td><td>5.7</td><td>16</td><td>16</td><td>16</td><td>16</td></tr></tbody></table>

evolved multifidelity kriging-based metamodel. The initial metamodel appears to be more effective than evofusion as it allows to calculate the same probability of failure with smaller CPU time. When the convergence rate is known, the estimation of an optimal mesh size for classification is proposed. Results obtained were quite similar, thus highlighting the very good performance of the initial strategy.

Future work will consist in extending this approach to more complex reliability problems. Indeed, an increased number of random variables, three dimensional problems with non linear mechanics or a non linear quantity of interest will be considered.

**Acknowledgments** The author would like to thank Florent Pled for the helpful and fruitful discussions. This work was carried out within the project MUSCAS (MUlti-SCAle Stochastic computation for MRE) granted by WEAMEC, West atlantic Marine Energy Community with the support of Région Pays de la Loire and in partnership with Chantiers de l'Atlantique.

## References

[1] Mark Ainsworth and J Tinsley Oden. A posteriori error estimation in finite element analysis. *Computer methods in applied mechanics and engineering*, 142(1-2):1-88, 1997.

[2] Kenneth F Alvin. Method for treating discretization error in nondeterministic analysis. *AIAA journal*, 38(5):910-916, 2000.
[3] R. Becker and R. Rannacher. A feed-back approach to error control in finite element methods: Basic analysis and examples. *Journal of Numerical Mathematics*, 4:237-264, 1996.

[4] Russel E. Caflisch. Monte carlo and quasi-monte carlo methods. *Acta Numerica*, 7:1-49, 1998.

[5] Nicolas Courrier, P-A Boucard, and Bruno Soulier. The use of partially converged simulations in building surrogate models. *Advances in Engineering Software*, 67:186-197, 2014.

[6] Nicolas Courrier, Pierre-Alain Boucard, and Bruno Soulier. Variable-fidelity modeling of structural analysis of assemblies. *Journal of Global Optimization*, 64(3):577-613, 2016.

[7] A. Dadone and B. Grossman. Progressive optimization of inverse fluid dynamic design problems. *Computers & Fluids*, 29(1):1-32, 2000.

[8] Jouke de Baar, Stephen Roberts, Richard Dwight, and Benoit Mallol. Uncertainty quantification for a sailing yacht hull, using multi-fidelity kriging. *Computers & Fluids*, 123:185-201, 2015.

[9] E. de Rocquigny. Structural reliability under monotony: Properties of form, simulation or response surface methods and a new class of monotonous reliability methods (mrm). *Structural Safety*, 31(5):363-374, 2009.

[10] Pedro Díez and Giovanni Calderón. Remeshing criteria and proper error representations for goal oriented h-adaptivity. *Computer methods in applied mechanics and engineering*, 196(4-6):719-733, 2007.

[11] B Echard, N Gayton, and M Lemaire. Ak-mcs: an active learning reliability method combining kriging and monte carlo simulation. *Structural Safety*, 33(2):145-154, 2011.

[12] B. Echard, N. Gayton, M. Lemaire, and N. Relun. A combined importance sampling and kriging reliability method for small failure probabilities with time-demanding numerical models. *Reliability Engineering & System Safety*, 111:232-240, 2013.

[13] Benjamin Echard. *Assessment by kriging of the reliability of structures subjected to fatigue stress*. PhD thesis, 2012.

[14] William Fauriat and Nicolas Gayton. Ak-sys: An adaptation of the ak-mcs method for system reliability. *Reliability Engineering & System Safety*, 123:137-144, 2014.
[15] Alexander IJ Forrester, Neil W Bressloff, and Andy J Keane. Optimization using surrogate models and partially converged computational fluid dynamics simulations. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 462(2071):2177-2204, 2006.

[16] L Gallimard. Error bounds for the reliability index in finite element reliability analysis. *International Journal for Numerical Methods in Engineering*, 87(8):781-794, 2011.

[17] L. Gallimard and J. Panetier. Error estimation of stress intensity factors for mixed-mode cracks. *International Journal for Numerical Methods in Engineering*, 68(3):299-316, 2006.

[18] B Gaspar, AP Teixeira, and C Guedes Soares. A study on a stopping criterion for active refinement algorithms in kriging surrogate models. *Safety and reliability of complex engineered systems*, pages 1219-1227, 2015.

[19] M. B. Giles. Multilevel Monte Carlo path simulation. *Oper. Res.*, 56(3):607-617, 2008.

[20] M. B. Giles. Multilevel monte carlo methods. *Acta Numerica*, 24:259, 2015.

[21] AA Griffith. The phenomena of rupture and flow in solids. *Philosophical Transactions of the Royal society of London*, 221:163-198, 1921.

[22] Abraham M Hasofer and Niels C Lind. Exact and invariant second-moment code format. *Journal of the Engineering Mechanics division*, 100(1):111-121, 1974.

[23] Marc C Kennedy and Anthony O'Hagan. Predicting the output from a complex computer code when fast approximations are available. *Biometrika*, 87(1):1-13, 2000.

[24] DG Krige. A statistical approach to some mine valuation and allied problems on the witwatersrand. *Mémoire de D.E.A.*, 1951.

[25] P. Ladevèze. Upper error bounds on calculated outputs of interest for linear and nonlinear structural problems. *Comptes Rendus Académie des Sciences - Mécanique, Paris*, 334(7):399-407, 2006.

[26] P. Ladevèze. Strict upper error bounds on computed outputs of interest in computational structural mechanics. *Computational Mechanics*, 42(2):271-286, 2008.
[27] P. Ladevèze and D. Leguillon. Error estimate procedure in the finite element method and application. *SIAM Journal of Numerical Analysis*, 20(3):485-509, 1983.

[28] P. Ladevèze and J. P. Pelle. *Mastering Calculations in Linear and Nonlinear Mechanics*. Springer, New York, 2004.

[29] Søren Nymand Lophaven, Hans Bruun Nielsen, and Jacob Søndergaard. *Aspects of the matlab toolbox DACE*. Citeseer, 2002.

[30] Sankaran Mahadevan and Ramesh Rebba. Inclusion of model errors in reliability-based optimization. *Journal of Mechanical Design*, 128(4):936-944, 2006.

[31] Gouranga Mallik, Martin Vohralík, and Soleiman Yousef. Goal-oriented a posteriori error estimation for conforming and nonconforming approximations with inexact solvers. *Journal of Computational and Applied Mathematics*, 366:112367, 2020.

[32] Georges Matheron. *Les variables régionalisées et leur estimation: une application de la théorie des fonctions aléatoires aux sciences de la nature*. Masson et CIE, 1965.

[33] Hermann G. Matthies and Andreas Keese. Galerkin methods for linear and nonlinear elliptic stochastic partial differential equations. *Computer Methods in Applied Mechanics and Engineering*, 194(12-16):1295-1331, April 2005.

[34] Nicholas Metropolis and Stanislaw Ulam. The monte carlo method. *Journal of the American statistical association*, 44(247):335-341, 1949.

[35] Stéphane Nachar, Pierre-Alain Boucard, David Néron, and Felipe Bordeu. Coupling multi-fidelity kriging and model-order reduction for the construction of virtual charts. *Computational Mechanics*, 2019.

[36] Anthony Nouy, Alexandre Clement, Franck Schoefs, and N Moës. An extended stochastic finite element method for solving stochastic partial differential equations on random domains. *Computer Methods in Applied Mechanics and Engineering*, 197(51-52):4663-4682, 2008.

[37] N. Parés, P. Díez, and A. Huerta. Subdomain-based flux-free a posteriori error estimators. *Computer Methods in Applied Mechanics and Engineering*, 195(4-6):297-323, 2006.
[38] Olivier Pasqualini, F Schoefs, Mathilde Chevreuil, and M Cazuguel. Measurements and statistical analysis of fillet weld geometrical parameters for probabilistic modelling of the fatigue capacity. *Marine Structures*, 34:226–248, 2013.

[39] P Perdikaris, D Venturi, JO Royset, and GE Karniadakis. Multi-fidelity modelling via recursive co-kriging and gaussian-markov random fields. *Proc. R. Soc. A*, 471(2179):20150018, 2015.

[40] F. Pled, L. Chamoin, and P. Ladevèze. On the techniques for constructing admissible stress fields in model verification: Performances on engineering examples. *International Journal for Numerical Methods in Engineering*, 88(5):409–441, 2011.

[41] Mohsen Rashki. Low-cost finite element method-based reliability analysis using adjusted control variate technique. page 10, 2017.

[42] Valentine Rey, Pierre Gosselet, and Christian Rey. Study of the strong prolongation equation for the construction of statically admissible stress fields: Implementation and optimization. *Computer Methods in Applied Mechanics and Engineering*, 268(0):82–104, 2014.

[43] Shmuel Rippa. An algorithm for selecting a good value for the parameter c in radial basis function interpolation. *Advances in Computational Mathematics*, 11(2-3):193–210, 1999.

[44] Marcus Rüter and Erwin Stein. Goal-oriented a posteriori error estimates in linear elastic fracture mechanics. *Computer methods in applied mechanics and engineering*, 195(4-6):251–278, 2006.

[45] Jerome Sacks, William J. Welch, Toby J. Mitchell, and Henry P. Wynn. Design and analysis of computer experiments. *Statistical Science*, 4(4):409–423, 11 1989.

[46] Jerome Sacks, William J Welch, Toby J Mitchell, and Henry P Wynn. Design and analysis of computer experiments. *Statistical science*, pages 409–423, 1989.

[47] Franck Schoefs. Sensitivity approach for modelling the environmental loading of marine structures through a matrix response surface. *Reliability Engineering & System Safety*, 93(7):1004–1017, 2008.

[48] M. Stern, E. B. Becker, and R. S. Dunham. A contour integral computation of mixed-mode stress intensity factors. *International Journal of Fracture*, 12(3):359–368, 1976.
[49] T Strouboulis, I Babuška, DK Datta, K Copps, and SK Gangaraj. A posteriori estimation and adaptive control of the error in the quantity of interest. part i: A posteriori estimation of the error in the von mises stress and the stress intensity factor. *Computer Methods in Applied Mechanics and Engineering*, 181(1-3):261-294, 2000.

[50] Hiroshi Tada, Paul C Paris, and George R Irwin. The stress analysis of cracks. *Handbook, Del Research Corporation*, 1973.

[51] Cao Tong, Zhili Sun, Qianli Zhao, Qibin Wang, and Shuang Wang. A hybrid algorithm for reliability analysis combining kriging and subset simulation importance sampling. *Journal of Mechanical Science and Technology*, 29(8):3183-3193, 2015.

[52] Larry Wasserman. *All of statistics: a concise course in statistical inference*. Springer Science & Business Media, 2013.

[53] Jinhao Zhang, Mi Xiao, Liang Gao, and Sheng Chu. A combined projection-outline-based active learning kriging and adaptive importance sampling method for hybrid reliability analysis with small failure probabilities. *Computer Methods in Applied Mechanics and Engineering*, 344:13 – 33, 2019.

[54] Jinhao Zhang, Mi Xiao, Liang Gao, and Junjian Fu. A novel projection outline based active learning method and its combination with kriging metamodel for hybrid reliability analysis with random and interval variables. *Computer Methods in Applied Mechanics and Engineering*, 341:32 – 52, 2018.