Uncertainty and sensitivity analysis of functional risk curves based on Gaussian processes

TL;DR

This paper introduces a Gaussian process-based method for efficiently constructing functional risk curves and analyzing the impact of input uncertainties, with confidence bounds and sensitivity analysis, for complex physical models.

Contribution

It presents a novel approach using Gaussian processes to build risk curves with confidence bounds and explores sensitivity analysis methods for input parameters.

Findings

Gaussian process regression effectively approximates risk curves.

Confidence bounds quantify approximation uncertainty.

PLI sensitivity indices reveal input parameter effects.

Abstract

A functional risk curve gives the probability of an undesirable event as a function of the value of a critical parameter of a considered physical system. In several applicative situations, this curve is built using phenomenological numerical models which simulate complex physical phenomena. To avoid cpu-time expensive numerical models, we propose to use Gaussian process regression to build functional risk curves. An algorithm is given to provide confidence bounds due to this approximation. Two methods of global sensitivity analysis of the models' random input parameters on the functional risk curve are also studied. In particular, the PLI sensitivity indices allow to understand the effect of misjudgment on the input parameters' probability density functions.

Equations

$$Y=G(a,X)=G(a,X_1,\dots ,X_d),$$

$$\forall a>0\text{mbox}POD(a)=P\left(Y>s\mid a\right),$$

$$\forall a>0\text{mbox}Frag(a)=P\left(Y>s\mid a\right),$$

$$P\left(Y>s|a\right)=\Phi \left(\frac{a-\alpha }{\beta }\right),$$

$$\tilde{Y}(a)={\beta }_0+{\beta }_{1}a+ϵ,$$

$$Y(a,X)={\beta }_0+{\beta }_{1}a+Z_{{\sigma }^2,\theta }(a,X),$$

$$\forall x\left(Y(a,x)\mid y_N\right)=\mathcal{N}\left(\hat{Y}(a,x),{\sigma }_Y^2(a,x)\right)$$

$$\Psi (a)=P\left(Y_N(a,X)>s|a\right),$$

$$\hat{\Psi }(a)=E_X\left[1-\Phi \left(\frac{s-\hat{Y}(a,X)}{{\sigma }_Y(a,X)}\right)\right].$$

$$\Psi (a)=P_X(G(a,X)>s|a)=\int 1_{G(a,x)>s}f(x)dx,$$

$$\Psi (a)\approx {\Psi }_{MC}(a)=\frac{1}{n}\sum _{i=1}^n{1}_{G(a,x^{(i)})>s},$$

$${\Psi }_{MC}(a)\approx {\Psi }_{MC,GP}(a)=\frac{1}{n}\sum _{i=1}^n{1}_{Y_N(a,X)>s}.$$

$$\sqrt{n}\left({\Psi }_{MC,GP}(a)-{\Psi }_{GP}(a)\right)⟶\mathcal{N}\left(0,{\Psi }_{GP}(a)(1-{\Psi }_{GP}(a))\right),$$

$${\Psi }_{GP}(a)=\int 1_{Y_N(a,X)>s}f(x)dx.$$

$${\psi }_{MC,GP}^{(j)}(a)=\frac{1}{n}\sum _{i=1}^n{1}_{y^{(j)}(a,x^{(i)})>s}.$$

$$\frac{1}{m\times n_{CLT}}\sum _{\stackrel{j=1,\dots ,m}{k=1,\dots ,n_{CLT}}}{\psi }_{MC,GP}^{(j,k)}(a),$$

$$\text{mbox}Var(Y)=V=\sum _{i=1}^d{V}_i+\sum _{1\le i<j\le d}{V}_{ij}+\dots +V_{1\dots d}$$

$$S_i=\frac{V_i}{V},S_{ij}=\frac{V_{ij}}{V},\dots ,\text{mbox}and{T}_i=S_i+\sum _{j\ne i}{S}_{ij}+\sum _{j<k,j\ne i,k\ne i}{S}_{ijk}+\dots .$$

$$\begin{array}[]{rcl}\Psi_{X}(a)&=&P_{X}(Y>s\mid a,X)=1-\Phi\left(\frac{s-\widehat{Y}(a,X)}{\sigma_{Y}(a,X)}\right)\;,\\ \Psi_{X_{i}}(a)&=&P_{X}(Y>s\mid a,X_{i})\;,\\ \Psi_{X_{-i}}(a)&=&P_{X}(Y>s\mid a,X_{-i})\;,\\ D&=&E_{X}\|\Psi(a)-\Psi_{X}(a)\|^{2}\\ &=&E_{X}\left[\int\left(\Psi(a)-\Psi_{X}(a)\right)^{2}\,da\right],\end{array}$$

$$\begin{array}[]{rcl}S_{i}^{\mbox{\tiny FRC}}&=&\displaystyle\frac{E_{X}\|\Psi(a)-\Psi_{X_{i}}(a)\|^{2}}{D}\;,\\ T_{i}^{\mbox{\tiny FRC}}&=&\displaystyle\frac{E_{X}\|\Psi_{X}(a)-\Psi_{X_{-i}}(a)\|^{2}}{D}\;.\end{array}$$

$$P_X\left(G(a,x)>s\right)=\int 1_{G(a,x)>s}f(x)dx,$$

$$P_{X_i,\delta }\left(G(a,x)>s\right)=\int 1_{G(a,x)>s}\frac{f_{i,\delta }(x_i)}{f_i(x_i)}f(x)dx,$$

$$S_{i,\delta}=\left\{\begin{array}[]{l}\displaystyle\frac{P_{X_{i},\delta}\left(G(a,X)>s\right)-P_{X}\left(G(a,X)>s\right)}{P_{X_{i},\delta}\left(G(a,X)>s\right)}\\ \mbox{if}\quad P_{X_{i},\delta}\left(G(a,X)>s\right)\geq P_{X}\left(G(a,X)>s\right),\\ \\ \displaystyle\frac{P_{X_{i},\delta}\left(G(a,X)>s\right)-P_{X}\left(G(a,X)>s\right)}{P_{X}\left(G(a,X)>s\right)}\\ \mbox{if}\quad P_{X_{i},\delta}\left(G(a,X)>s\right)<P_{X}\left(G(a,X)>s\right).\\ \end{array}\right.$$

$$\Psi (a)=E_X\left[1-\Phi \left(\frac{s-\hat{Y}(a,X)}{{\sigma }_Y(a,X)}\right)\right]=\int 1-\Phi \left(\frac{s-\hat{Y}(a,x)}{{\sigma }_Y(a,x)}\right)f(x)dx.$$

$${\Psi }_{i,\delta }(a)=\int 1-\Phi \left(\frac{s-\hat{Y}(a,x)}{{\sigma }_Y(a,x)}\right)\frac{f_{i,\delta }(x_i)}{f_i(x_i)}f(x)dx.$$

$$S_{i,\delta}(a)=\left\{\begin{array}[]{l}\displaystyle\frac{\Psi_{i,\delta}(a)-\Psi(a)}{\Psi_{i,\delta}(a)}\quad\mbox{if}\quad\Psi_{i,\delta}(a)\geq\Psi(a),\\ \\ \displaystyle\frac{\Psi_{i,\delta}(a)-\Psi(a)}{\Psi(a)}\quad\mbox{if}\quad\Psi_{i,\delta}(a)<\Psi(a).\\ \end{array}\right.$$

$$KL(f_{i,\delta }(x_i),f_i(x_i))={\int }_{-\infty }^{+\infty }{f}_{i,\delta }(x_i)\log \frac{f_{i,\delta }(x_i)}{f_i(x_i)}d{x}_i.$$

Full text

††journal: Reliability Engineering & System Safety

1 Introduction

In industrial practice, the estimation of a functional risk curve (FRC) is often required as a quantitative measure of the safety of a system [1]. A FRC (also called “functional risk criterion” or “fragility curve” in some application domains [2]) gives the probability of an undesirable event as a function of the value of a critical parameter of a considered physical system. The estimation of this curve sometimes relies on deterministic phenomenological computer models which simulate complex physical phenomena. Uncertain input parameters of this computer code are modeled as random variables while one of its scalar output variables becomes the studied random variable of interest.

As it is based on a probabilistic modeling of uncertain physical variables and their propagation through a numerical model, this problem can be directly related to the uncertainty management methodology in numerical simulation [1, 3]. This methodology proposes a generic framework of modeling, calibrating, propagating and prioritizing uncertainty sources through a numerical model (or computer code). Indeed, investigation of complex computer code experiments has remained an important challenge in all domains of science and technology, in order to carry on simulations as well as predictions, uncertainty analysis or sensitivity studies [4, 5]. In this paper, the physical numerical model $G$ is expressed as

[TABLE]

with $a\in\mathbb{R}$ the critical input parameter of the FRC and $X\in{\mathbb{R}}^{d}$ the random input vector of dimension $d$ and $Y\in\mathbb{R}$ a scalar model output. In this paper, parameters $a$ and $X$ are considered as independent. However the presented results hold when $a$ and $X$ are dependent and the distribution of $X$ conditionally to $a$ is known.

However, standard uncertainty treatment techniques require many model evaluations and a major algorithmic difficulty arises when the computer code under study is too time expensive to be directly used. For cpu-time expensive models, one solution consists in replacing the numerical model by a mathematical approximation, called a response surface or a metamodel. Several statistical tools based on numerical design of experiments, efficient algorithms for uncertainty propagation and metamodeling concepts will then be useful [4]. In this paper, Gaussian process regression [6] is used as a metamodeling technique (surrogate model of the computer code) and applied in the particular context of a FRC as a quantity of interest. Numerous studies have shown that this model provides a powerful statistical framework to compute an efficient predictor of a deterministic computer code response [7, 8].

Associated to the estimation of the quantity of interest after the uncertainty propagation step, the sensitivity analysis step is performed to determine those parameters that mostly influence the model response [9]. In particular, global sensitivity analysis methods take into account the overall uncertainty ranges of the input parameter (see [10] for a recent review). Several works have focused on estimating global sensitivity indices (especially the variance-based ones, known as the Sobol’ indices) using Gaussian process models [11, 12, 8]. In this paper, we present global sensitivity indices attached to the whole FRC (firstly introduced in [13] in a particular context), while showing how to develop them with a Gaussian process model. We also focus on a recently developed method, the Perturbed-Law based sensitivity Indices (PLI) of [14], which seems promising in terms of computational efficiency and interpretation.

Two examples of FRC largely used in some industrial safety practices are given in the next section. The third section describes the Gaussian process way to model and estimate a FRC. The fourth section develops the sensitivity indices adapted to FRC.

2 Motivating examples and classical methods

2.1 Probability of detection curves

In several industries, the Probability Of Detection (POD) curve is a standard tool to evaluate the performance of Non Destructive Testing (NDT) procedures. The goal is to assess the inspection capability for the detection of flaws in the inspected structure. For instance, new aeronautic regulations require to perform appropriate damage tolerance assessments to address the potential for failure due to critical parts of material and manufacturing variability. This imposes enhanced requirement on POD sizing. Another example is the eddy current non destructive examination process which is used to ensure integrity of steam generator tubes in nuclear power plants [15, 16, 13]. Experimental campaigns and simulated experiments are used in order to provide POD curves, which can be used to demonstrate the performance of inspection process to regulatory authorities.

In practice, the high costs of the implementation of experimental POD campaigns combined with the continuous increase in the complexity of configuration make them sometimes unaffordable. To overcome this problem, it is possible to resort to numerical simulation of NDT processes. This approach is commonly called MAPOD for “Model Assisted Probability of Detection” [17, 18, 13].

More precisely, in the POD context, the problem is formulated as follows. Given a threshold $s>0$, a flaw is considered to be detected when $Y>s$ with $Y$ the signal amplitude and $s$ a detection threshold. Therefore the one dimensional POD curve is denoted by:

[TABLE]

where $a$ is the POD parameter of interest (for example the size flaw) and $X$ are the other input parameters. When simulated experiments are used to build the POD (see for example [13]), $Y$ is a scalar output of a numerical model (see Eq. 1)), where $a$ is determined by its bounds and $X$ is a random vector defined by its joint probability density function (pdf).

In Sections 3 and 4 of this paper, we will illustrate the mathematical methods for FRC building and analysis on a case study on the eddy current examination of the wear induced on steam generator tubes by anti-vibration bars [15]. For details, the reader is referred to [16, 13]. The study involves the finite element computer code C3D [19]. The NDT process under study is the examination of anti-vibration bars (AVB) subject to wear of steam generator tubes inside a nuclear power plant [15]. The computer code used simulates by finite elements the eddy-current propagation phenomena. The critical input parameter $a$ is the maximum of the two flaw sizes while five other input parameters (vector $X$) of the computer code are also considered (with their associated pdf): E (pipe thickness), ebav1 and ebav2 (lengths of the gap between the AVB and the first and second flaws), h11 and h12 (first and second flaw heights). The output $Y$ of the computer code is the signal amplitude, and the threshold of detection $s$ is fixed at a specific value.

2.2 Seismic fragility curves

The “fragility curve” is a popular functional risk criterion, commonly used in many engineering fields [2]. It describes the probability that the actual damage to a structure exceeds a damage threshold, when the structure is assigned to a specified load intensity. For instance, the seismic fragility curve, which concerns systems subject to earthquakes, is of particular interest in nuclear safety studies (see [20]).

In the case of seismic risk assessment the load is usually expressed as a scalar characteristic of a seismic signal, typically the horizontal Peak Ground Acceleration (PGA), which is a common choice in civil engineering. Several parameters and phenomena, distinct from the PGA, also influence the load: for each PGA value, the occurrence of the damage event is random. The fragility curve may then be interpreted as the cumulative distribution function of the “structural capability”, i.e. the maximum load the structure under investigation can bear without damage [21, 22]. Fragility curves are useful tools in structural analysis as they provide more complete information than the usual failure probability (established for a reference value of the load only).

In standard practice, the assessment is made either following an approach entirely based on available expertise or by the statistical analysis of actually observed or simulated data. As actual damage data may be scarce due to the rarity of severe earthquakes liable to generate damages on highly safe structures, observations are generated by mock-up or (most often) numerical experiments.

In this context, a fragility function is expressed as

[TABLE]

where $Y$ is the variable which characterizes the structural damage, $s=Y_{\mbox{\tiny max}}$ is the maximal capacity of the structure and $a$ is an intensity measure characterizing the ground motion severity (e.g. peak ground acceleration). When the fragility curve is estimated by numerical simulation (see for example [20]), $Y$ is the numerical model output (see Eq. (1)), $a$ is determined by its bounds and the random vector $X$ is defined by its joint pdf.

2.3 Classical methods

The standard practice for establishing a fragility curve consists in assuming a lognormal distribution for $Y$ [23]. The classical approach for determining a POD curve is called the Berens model and relies on assuming a simple linear relation (with Gaussian residuals) between $\log(Y)$ and $a$ [24]. In fact, both cases are equivalent and lead to the following form for the FRC:

[TABLE]

where $\Phi$ is the standard Gaussian distribution function, and $\alpha$ and $\beta$ are two constants to be estimated (by the maximum likelihood estimation for example) as functions of the $Y$-data sample and $s$. The parameter $a$ can also be expressed in terms of a logarithm in this expression. We note that if the Berens method is applied in a FRC framework (e.g. with $Y$ the maximal displacement of a building and $a$ the pick ground acceleration of a seismic signal) the resulting PoD curve is the classical FRC.

As shown in [13], this model can easily be improved by considering a Box-Cox transformation [25] instead of a simple logarithmic one on $Y$: $\tilde{Y}=\frac{Y^{\lambda}-1}{\lambda}$. $\lambda$ is determined by maximum likelihood as the real number that offers the best linear regression of $\tilde{Y}$ regarding the parameter $a$. The model is now based on $\tilde{Y}$ and is defined as

[TABLE]

with $\epsilon$ the model error such as $\epsilon\sim\mathcal{N}\left(0,\sigma_{\epsilon}^{2}\right)$. The maximum likelihood method provides the estimators $\hat{\beta_{0}},\hat{\beta_{1}}\;\text{and}\;\hat{\sigma_{\epsilon}}$. Hence the model implies the following result: $\forall a>0,\quad\tilde{Y}(a)\sim\mathcal{N}\left(\hat{\beta_{0}}+\hat{\beta_{1}}a,\hat{\sigma_{\epsilon}}^{2}\right)$. Then the value of the FRC can be estimated from $\hat{\beta_{0}},\hat{\beta_{1}}\;\text{and}\;\hat{\sigma_{\epsilon}}$, as its confidence interval from the property of the maximum likelihood estimators. Note that it is possible to relax the Gaussian hypothesis on $\epsilon$ by using non-parametric approaches [26, 13]

3 Gaussian-process based functional risk curve estimation

Classical methods described in Section 2.3 use a simple model between $Y$ and $a$ without modeling the functional dependence of $Y$ on the other uncertain variables $X$. As in our case data $y_{N}$ (sample of $Y$ of size $N$) are obtained from numerical models, the input variables $X$ are controlled and can be introduced in the FRC determination process. However, with cpu-time expensive numerical models, model evaluations can be somewhat limited and only small samples of $Y$ can be obtained. As shown in [27, 21, 13], the use of a metamodel, which is a mathematical approximation of the computer code [4], is useful. We highlight that it is important in a safety framework to evaluate accurately the uncertainty on the surrogate model. We propose in this section a general formulation for FRC that uses Gaussian process regression as a metamodel [28, 6]. We propose in this section a general formulation for FRC that uses the Gaussian process regression metamodel [28, 6]. Through the Gaussian assumption, this metamodel provides a convenient framework to perform uncertainty quantification on the quantity of interest estimate. This formulation has been introduced in [13] without the full mathematical algorithm which is developed in the present paper.

Let $X$ be the random vector of influential and uncertain parameters of the computer model $G(\cdot)$ and $a$ be the parameter of interest (the abscissa of the FRC). The prior knowledge on $G(a,X)$ is modeled by $Y(a,X)$ and defined as follows .

[TABLE]

where $Z_{\sigma^{2},\theta}$ is a centered Gaussian process. We make the assumption that $Z_{\sigma^{2},\theta}$ is second order stationary with variance $\sigma^{2}$ and a parametric covariance function (which corresponds to the kernel in other statistical learning methods). As covariance choice, one can cite the popular Matérn 5/2 parametrized by its lengthscale $\theta$). Thanks to the maximum likelihood method, we can estimate the values of the so far-unknown parameters: $\beta_{0},\beta_{1},\sigma^{2}$ and $\theta$ (see for instance [7] for more details).

Gaussian process regression (also known as the kriging process) provides an estimator of $G(a,x)$ which is called the kriging predictor and written $\widehat{Y}(a,x)$. In addition to the kriging predictor, the kriging variance $\sigma_{Y}^{2}(a,x)$ quantifies the uncertainty induced by estimating $Y(a,x)$ with $\widehat{Y}(a,x)$. The predictive distribution is given by $Y(a,x)$ conditioned by $y_{N}$:

[TABLE]

where $\widehat{Y}(a,x)$ (the kriging mean) and $\sigma_{Y}^{2}(a,x)$ (the kriging variance) can both be explicitly estimated.

Obtaining the FRC $\Psi(a)$ consists in replacing $Y=G(a,X)$ by its Gaussian process metamodel (7) in (2). Hence we can estimate the value of $\Psi(a)$, for $a>0$ from:

[TABLE]

where $Y_{N}(a,X)$ is a Gaussian process which follows the distribution $(Y(a,X)\mid y_{N})$ and $s$ is a fixed threshold value. Two sources of uncertainty have to be taken into account in (8): the first coming from the parameter $X$ and the second coming from the Gaussian distribution in (7).

From (8), the following estimate for $\Psi(a)$ can be deduced:

[TABLE]

This equation corresponds to the mean FRC with respect to the Gaussian process metamodel. The expectation in (9) is estimated using a classical Monte Carlo integration procedure.

In order to estimate the uncertainty on the FRC estimation, we start from its integral expression:

[TABLE]

where $f(x)$ is the joint pdf of $X$ (independent on $a$). The first uncertainty source on $\Psi(a)$ comes from the numerical evaluation of the integral inside the $\Psi(a)$ definition. This evaluation is done by the Monte Carlo method:

[TABLE]

where $(x^{(i)})_{i=1,\dots,n}$ is a sample of the random variable $X$. The required size of this sample can be too large to use this formula in practice (case of a costly evaluation of the code). We then replace the code by its Gaussian process approximation $Y_{N}(a,X)$:

[TABLE]

Let us recall that $Y_{N}(a,X)$ is a Gaussian process with known mean and variance. To compute the integration error, the central limit theorem is used and gives for $n\rightarrow\infty$:

[TABLE]

where

[TABLE]

$\Psi_{GP}(a)$ corresponds to the FRC for the process $Y_{N}(a,X)$. The algorithm uses a significantly large value of $n$ in order to deal with a valid Gaussian approximation.

The second uncertainty source on $\Psi(a)$ estimation comes from Gaussian process approximation. We simulate $m$ realizations $(y^{(j)}(a,x))_{j=1,\dots,m}$ from $({Y(a,x)|y_{N}})$ in order to evaluate the variability of $\Psi_{MC,GP}(a)$ which comes from Gaussian process approximation. Thus, we compute:

[TABLE]

Then, for each $\psi_{MC,GP}^{(j)}(a)$, $j=1,\dots,m$, we compute the Monte Carlo error using the central limit theorem cited before in (13) where the variance of the Gaussian distribution in (13) is estimated with $\psi_{MC,GP}^{(j)}(a)(1-\psi_{MC,GP}^{(j)}(a))$, $j=1,\dots,m$.

From (13), we compute for each $j=1,\dots,m$ a sample of size $n_{CLT}$. From this double Monte Carlo method, we obtain a sample $\left[\psi_{MC,GP}^{(j,k)}(a)\right]_{\stackrel{{\scriptstyle j=1,\dots,m}}{{k=1,\dots,n_{CLT}}}}$ of size $m\times n_{CLT}$. Then the FRC estimate is given by :

[TABLE]

and the uncertainty is evaluate from $\left[\psi_{MC,GP}^{(j,k)}(a)\right]_{\stackrel{{\scriptstyle j=1,\dots,m}}{{k=1,\dots,n_{CLT}}}}$ by using classical quantile estimators.

From this double Monte Carlo method, we are able to compute the FRC uncertainty due to the Gaussian process metamodel and due to the numerical integration process. Figure 1 illustrates the algorithm on a POD curve estimation in an NDT application (see section 2.1). A Gaussian process metamodel has been built on $N=100$ numerical simulations of eddy current examinations of steam generator tubes [13]. We visualize the confidence interval induced by the Monte Carlo (MC) estimation ($n=10,000$, $n_{CLT}=100,000$), the one induced by the Gaussian process approximation ($m=3,000$) and the total confidence interval (including both approximations: GP+MC). It could be noticed that from the double Monte-Carlo sample $\left[\psi_{MC,GP}^{(j,k)}(a)\right]_{\stackrel{{\scriptstyle j=1,\dots,m}}{{k=1,\dots,n_{CLT}}}}$, one can evaluate with an ANOVA decomposition the contribution of the Gaussian process and the parameters $X$ on the FRC estimate mean squared error (see [12]).

Choosing $m$, the number of conditional Gaussian process realizations to simulate, is a difficult issue of this process. As explained in [12] and other papers, there are some computational constraints (in memory and cpu time) when simulating conditional Gaussian processes on large samples. The choice of $m=3000$ has been chosen as a sufficient value (that we avoid to be too large) in order to estimate a $95\%$ lower bound (due to the Gaussian process error). Then $n$ has been tuned in order to have a smaller error induced by the Monte Carlo sampling process than the one induces by the Gaussian process.

4 Global sensitivity analysis of functional risk curves

The objective of sensitivity analysis is to determine those input variables that mostly influence the model response [9, 29, 30]. Global sensitivity analysis methods take into account the overall uncertainty of the input parameters. Previous works on the POD [31] and seismic fragility context [32] have only considered sensitivity analysis on the variance and distribution of the model output variable $Y$. In this section, we propose two global sensitivity indices attached to the whole FRC: one already introduced in [13] (aggregated Sobol’ indices) and a new one based on a completely different idea (perturbed-law based indices). Therefore, the influence of the parameter of interest $a$, as it is the FRC abscissa, will not be considered as an input variable whose sensitivity index has to be computed.

4.1 Aggregated Sobol’ indices

We first focus on the variance-based sensitivity formulation [33], which is one the most popular tools that has been proven to be robust and easily interpretable [10]. The associated sensitivity indices are often called Sobol’ indices.

In the case of independent inputs, using the Sobol-Hoeffding decomposition [34], the variance of the numerical model $Y=G(X_{1},\ldots,X_{d})$ is decomposed into the following sum:

[TABLE]

with $V_{i}=\mbox{Var}[E(Y|X_{i})]$, $V_{ij}=\mbox{Var}[E(Y|X_{i},X_{j})]-V_{i}-V_{j}$, etc. Then, $\forall i,j=1\ldots d,\;i<j$, the Sobol’ indices are given by

[TABLE]

The first-order Sobol’ index $S_{i}$ measures the individual effect of the input $X_{i}$ on the variance of the output $Y$, while the total Sobol’ index $T_{i}$ measures the $X_{i}$ effect and all the interaction effects between $X_{i}$ and the other inputs (as the second-order effect $S_{ij}$). $T_{i}$ can be reinterpreted as $T_{i}=\displaystyle 1-\frac{V_{-i}}{V}$ with $V_{-i}=\mbox{Var}[E(Y|X_{-i})]$ and $X_{-i}$ is the vector of all inputs except $X_{i}$. Each Sobol’ index is then interpreted in terms of percentage of total variance explanation.

In order to define similar sensitivity indices for the whole FRC (Eq. (9)), [13] defines the following quantities:

[TABLE]

with $\|.\|$ the euclidean norm. As in [35] (see also [36]) which deals with functional model outputs, our objective is to obtain a single sensitivity index for each input by averaging the variability of the function, here the FRC depending on $a$. The specificity here is that we deal with a FRC function. Therefore, [13] defines the FRC aggregated Sobol’ indices by:

[TABLE]

$S_{i}^{\mbox{\tiny FRC}}$ (resp. $T_{i}^{\mbox{\tiny FRC}}$) gives the first-order (resp. total) effect of $X_{i}$ on the mean FRC.

These FRC aggregated Sobol’ indices are easily computed by the metamodels. In particular, the kriging metamodel allows to replace $P_{X}(Y>s\mid a)$ by the expectation $\displaystyle E_{X}\left[1-\Phi\left(\frac{s-\widehat{Y}(a,X)}{\sigma_{Y}(a,X)}\right)\right]$ in the FRC expressions of (18). In order to compute the Sobol’ indices, we use the classical pick-freeze formulas (see for example [34]). The confidence intervals due to the finite Monte Carlo sample are then obtained by bootstrap. Let us noticed that it is also possible to estimate the error on the aggregated Sobol’ indices due to the kriging error as in [12] for the classical Sobol’ indices. This computationally heavy process has not been applied here.

The sensitivity analysis results on our use case are given in Figure 2 with the FRC aggregated Sobol’ indices of the five physical input parameters (called E, ebav1, ebav2, h11 and h12). ebav1 is the main influent parameter on the FRC (the POD curve of Figure 1), followed by h12 and ebav2. $E$ and h11 have no influence.

The FRC aggregated Sobol’ indices quantify the sensitivity of each input on the overall curve. However, we could be interested in the sensitivities on the FRC at a specific value of the parameter of interest $a$. This can be directly done by replacing $Y$ by $\Psi_{X}(a)$ in the equations (16) and (17). If we are now interested by the sensitivities at a specific probability value of the FRC, we have to study the inverse function of the FRC: $\Psi_{X}^{-1}(p)$ with $p$ a given probability. Similarly to the previous case, the Sobol’ indices can be obtained by replacing $Y$ by $\Psi_{X}^{-1}(p)$ in the equations (16) and (17). Figure 3 gives these sensitivity indices on our data for $p=0.90$ which corresponds to the quantile of $a$ at order $90\%$ (noted $a_{90}$). We see that the influences of the inputs on $a_{90}$ are very close to the ones on the overall POD curve. It could be explained by the fact that the different inputs have rather linear effects on the output of the model (the amplitude).

4.2 Perturbed-law based sensitivity indices

We propose now to quantify the impact on the FRC of a perturbation of the input parameters pdfs by answering to the following question: “What would be the FRC if the pdf of the $i^{\mbox{\tiny th}}$ input $X_{i}$ had been modified?”. In this approach, all the input parameters are modeled by random variables, and their input probability densities are supposed to be unknown. Let us remark that a negligible sensitivity index of an input will not allow to fix this input, but just to say that its pdf has no influence on the FRC.

We use the so-called Perturbed-Law based sensitivity Indices (PLI) measures recently introduced in [14] (see also [37]). We start from the integral-form of the FRC:

[TABLE]

where $f(x)$ is the joint pdf of $X$. Modifying the pdf $f_{i}(x_{i})$ of $X_{i}$ gives us $f_{i,\delta}(x_{i})$, the perturbed pdf of $X_{i}$. After this perturbation, the FRC, denoted $P_{X_{i},\delta}(\cdot)$ instead of $P_{X}(\cdot)$, can be written as:

[TABLE]

The PLI measures only consist of the comparison of the FRCs before and after the perturbation, and are defined by:

[TABLE]

A negative $S_{i,\delta}$ means that the FRC is smaller after the perturbation, while a positive $S_{i,\delta}$ means that the FRC has increased. The estimation of $P_{X_{i},\delta}\left(G(a,X)>s\right)$ is based on reverse importance sampling [38]. Asymptotical properties of the estimators give also some confidence intervals on the PLI measures.

In [14], the numerical model is directly used to estimate the PLI measures by large Monte Carlo samples. We propose here to estimate the PLI measures using the Gaussian process metamodel and integrating its error in the estimates. The mean FRC which we consider is given by (9):

[TABLE]

After the perturbation of the pdf of the $i^{\mbox{\tiny th}}$ input, the mean FRC is given by:

[TABLE]

PLI measures are then given by:

[TABLE]

The last element of the PLI method is the definition of the perturbations $f_{i,\delta}(x_{i})$. [14] choose to perturb a statistical characteristic (for example the mean, or the variance, or a quantile, …) of $X_{i}$ in order to be “as close as possible” to the initial pdf $f_{i}(x_{i})$. The dissimilarity measure between $f_{i,\delta}(x_{i})$ and $f_{i}(x_{i})$, which contains the required properties, is the Kullback-Leibler divergence:

[TABLE]

It implies that $f_{i,\delta}(X_{i})$ and $f_{i}(X)$ have the same definition domain (as a consequence, the input domain bounds cannot be changed). The Kullback-Leibler divergence has the advantage to be easily minimized [39], in order to obtain explicit solutions for $f_{i,\delta}(X_{i})$ for classical pdf (e.g. Gaussian).

Figure 4 gives some examples of perturbations for the uniform pdf on $[0,1]$ (mean is $0.5$ and variance is $1/12$).

Figure 5 gives some examples of perturbations for the standard Gaussian pdf (mean is [math] and variance is $1$).

Figure 6 gives an example of a PLI-based sensitivity analysis on our NDT test case (see section 2.1) which aims to estimate POD curves and the associated sensitivity indices to its five physical input parameters (called E, ebav1, ebav2, h11 and h12), the defect size $a$ being fixed here at a given value. The input parameters pdfs are all uniform on $[0,1]$. Their non-perturbed mean is then $1/2$. Then, the mean of each input is modified in the range of $\delta\in[0.1,0.9]$. For each $\delta$ value and each influential parameter, Figure 6 shows the PLI estimates $\hat{S}_{\delta}$. Large absolute values of $\hat{S}_{\delta}$ imply a large impact of the perturbation on the FRC. Moreover, the sign of $\hat{S}_{\delta}$ indicates if the new probability has decreased (negative case) or increased (positive case). Confidence intervals shown in Figure 6 are obtained using a Gaussian asymptotic property of the estimates $\hat{S}_{i,\delta}$ [14]. In Figure 6, one can see for instance that increasing the h12 input mean largely increases the FRC, while decreasing the ebav1 input mean largely decreases the FRC. In contrary, perturbations of the $E$ mean have no effect on the FRC.

In order to have a global view of the FRC sensitivity for different defect sizes $a$, Figure 7 gives the PLI graphs for each input, with variations of $a$ and the $\delta$ perturbation. Clearly, the inputs ebav1 and h12 have strong impacts on the FRC, but especially in the $a<0.3$ range. In contrary, the inputs $E$ and h11 have no influence. These elements are key points for the engineers to understand the effects of the physical parameters on the POD curve.

5 Conclusion

In this paper, several methods of uncertainty and sensitivity analysis of model outputs have been developed on the particular quantity of interest of the functional risk curve (FRC). A FRC provides the probability of an undesirable event as a function of the value of a critical input parameter of a considered physical system. Focus has been on the use of the Gaussian process metamodels in order to build FRCs from numerical experiments. This approach is useful when the computer model is expensive to evaluate in such a way that only a small sample of the model output can be obtained. In addition to the mean risk curve, the metamodel allows to obtain the confidence bands via conditional Gaussian process simulations.

One important advantage of using a metamodel is that sensitivity analysis is facilitated. We have defined two kinds of sensitivity indices related to the FRC as the quantity of interest. First, we have formulated the FRC aggregated Sobol’ indices, which are variance-based measures. Second, based on perturbation of the pdf of each model input variable, sensitivity indices (called PLI) of the model inputs on the FRC are also proposed. These allow to understand the effect of misjudgment on the pdf of each input parameters.

An example taken from simulated NDE inspections highlights the added value of the FRC, which in this context correspond to POD curves. The two proposed sensitivity analysis methods provide strongly complementary information to experts interpreting the NDE results. However, FRCs are used in many other engineering contexts, e.g. in the seismic fragility assessment (as shown in section 2.2) and in the evaluation of hydraulic works reliability subject to extremely high water levels [2]. These tools would also be useful in separating the effects of stochastic input variables and epistemic input parameters. For instance, this problem is highlighted in [40] and [41].

Finally, two main mathematical perspectives are identified from this work. First, adaptive designs can be developed: in the spirit of [42], it would consist of using the Gaussian process model in order to define some SUR (“Sequential Uncertainty Reduction”) criteria on the FRC as the quantity of interest. Optimizing such a criterion would provide a new set of input parameter values, which would be run with the computer code in order to decrease the POD confidence interval. Second, [43] have started a preliminary study that considers FRC as a random distribution function. This framework allows to deal with stochastic computer codes instead of the deterministic ones of this paper (see for example Kleijnen [44]).

Acknowledgments

Part of this work has been backed by French National Research Agency (ANR) through project ByPASS ANR-13-MONU-0011. We acknowledge three anonymous reviewers for their deep analysis of the paper and the useful comments they provide. We are also grateful to Alberto Pasanisi for giving helpful discussions and Andreas Schumm for his help with the English language.