Global sensitivity analysis in the context of imprecise probabilities (p-boxes) using sparse polynomial chaos expansions

TL;DR

This paper extends global sensitivity analysis to models with imprecise probabilities using p-boxes, introducing imprecise Sobol' indices computed efficiently via sparse polynomial chaos expansions.

Contribution

It proposes a novel algorithm for calculating imprecise Sobol' indices with low computational cost using augmented space and sparse PCE.

Findings

Validated the method on analytical examples

Achieved significant computational efficiency

Compared results with brute-force Monte Carlo

Abstract

Global sensitivity analysis aims at determining which uncertain input parameters of a computational model primarily drives the variance of the output quantities of interest. Sobol' indices are now routinely applied in this context when the input parameters are modelled by classical probability theory using random variables. In many practical applications however, input parameters are affected by both aleatory and epistemic (so-called polymorphic) uncertainty, for which imprecise probability representations have become popular in the last decade. In this paper, we consider that the uncertain input parameters are modelled by parametric probability boxes (p-boxes). We propose interval-valued (so-called imprecise) Sobol' indices as an extension of their classical definition. An original algorithm based on the concepts of augmented space, isoprobabilistic transforms and sparse polynomial chaos expansions is devised to allow for the computation of these imprecise Sobol' indices at extremely low cost. In particular, phantoms points are introduced to build an experimental design in the augmented space (necessary for the calibration of the sparse PCE) which leads to a smart reuse of runs of the original computational model. The approach is illustrated on three analytical and engineering examples which allows one to validate the proposed algorithms against brute-force double-loop Monte Carlo simulation.

Tables

Table 1. Table 1: SDOF oscillator – definition of the parametric p-boxes in the input vector.

$X_i$	Distribution		Mean value	Standard deviation
$r$	Gaussian		$[0.49,0.51]$	$0.05$
$F_1$	Gaussian		$[0.8,1.2]$	$0.2$
$t_1$	Gaussian		$[0.95,1.05]$	$0.2$
$c_1$	Gaussian		$1$	$0.1$
$c_2$	Gaussian		$0.1$	$0.01$
$m$	Gaussian		$1$	$0.05$

Table 2. Table 2: SDOF oscillator – imprecise first-order Sobol’ indices estimated by augmented PCE ( n p ​ h = 10 subscript 𝑛 𝑝 ℎ 10 n_{ph}=10 ).

	Reference	Augmented PCE
		$N=30$	$N=50$	$N=100$	$N=200$
$S_r$	$[0.220,\mathrm{\hspace{0.17em}0.307}]$	$[0.312,\mathrm{\hspace{0.17em}0.313}]$	$[0.223,\mathrm{\hspace{0.17em}0.305}]$	$[0.220,\mathrm{\hspace{0.17em}0.307}]$	$[0.220,\mathrm{\hspace{0.17em}0.307}]$
$S_{F_1}$	$[0.308,\mathrm{\hspace{0.17em}0.459}]$	$[0.373,\mathrm{\hspace{0.17em}0.375}]$	$[0.306,\mathrm{\hspace{0.17em}0.450}]$	$[0.308,\mathrm{\hspace{0.17em}0.459}]$	$[0.308,\mathrm{\hspace{0.17em}0.460}]$
$S_{t_1}$	$[0.215,\mathrm{\hspace{0.17em}0.413}]$	$[0.214,\mathrm{\hspace{0.17em}0.215}]$	$[0.224,\mathrm{\hspace{0.17em}0.415}]$	$[0.216,\mathrm{\hspace{0.17em}0.414}]$	$[0.215,\mathrm{\hspace{0.17em}0.413}]$
$S_{c_1}$	$[0.017,\mathrm{\hspace{0.17em}0.034}]$	$[0.030,0.030]$	$[0.017,\mathrm{\hspace{0.17em}0.032}]$	$[0.017,\mathrm{\hspace{0.17em}0.034}]$	$[0.017,\mathrm{\hspace{0.17em}0.034}]$
$S_{c_2}$	$[0.000,\mathrm{\hspace{0.17em}0.000}]$	$[0.001,\mathrm{\hspace{0.17em}0.001}]$	$[0.000,\mathrm{\hspace{0.17em}0.001}]$	$[0.000,\mathrm{\hspace{0.17em}0.000}]$	$[0.000,\mathrm{\hspace{0.17em}0.000}]$
$S_m$	$[0.003,\mathrm{\hspace{0.17em}0.006}]$	$[0.008,0.009]$	$[0.003,\mathrm{\hspace{0.17em}0.006}]$	$[0.003,\mathrm{\hspace{0.17em}0.006}]$	$[0.003,\mathrm{\hspace{0.17em}0.006}]$
${\widehat{err}}_{gen}$		$1.05\cdot {10}^{-1}$	$2.39\cdot {10}^{-4}$	$1.63\cdot {10}^{-5}$	$1.65\cdot {10}^{-6}$

Equations

$$\mathcal{M}:\mathbf{x}\in {\mathcal{D}}_{\mathbf{X}}\subset {\mathbb{R}}^M\mapsto y\in {\mathcal{D}}_Y\subset \mathbb{R}.$$

$$Y\approx {\mathcal{M}}^{(\text{P})}\left(\mathbf{X}\right)=\sum _{\mathbf{\alpha }\in \mathcal{A}}{a}_{\mathbf{\alpha }}{\psi }_{\mathbf{\alpha }}\left(\mathbf{X}\right),$$

$$⟨{\varphi }_1,{\varphi }_2{⟩}_i={\int }_{{\mathcal{D}}_{X_i}}{\varphi }_1(x){\varphi }_2(x)f_{X_i}(x)\text{d}x,$$

$$⟨P_j^{(i)},P_k^{(i)}⟩={\int }_{{\mathcal{D}}_{X_i}}{P}_j^{(i)}(x)P_k^{(i)}(x)f_{X_i}(x)\text{d}x={\delta }_{jk},$$

$${\psi }_{\mathbf{\alpha }}\left(\mathbf{x}\right)=\prod _{i=1}^M{P}_{{\alpha }_i}^{(i)}\left(x_i\right).$$

$$\mathbb{E}\left[{\psi }_{\mathbf{\alpha }}\left(\mathbf{X}\right){\psi }_{\mathbf{\beta }}\left(\mathbf{X}\right)\right]\stackrel{\text{def}}{=}{\int }_{{\mathcal{D}}_{\mathbf{x}}}{\psi }_{\mathbf{\alpha }}\left(\mathbf{x}\right){\psi }_{\mathbf{\beta }}\left(\mathbf{X}\right)f_{\mathbf{X}}\left(\mathbf{x}\right)\text{d}\mathbf{x}={\delta }_{\mathbf{\alpha }\mathbf{\beta }},$$

$${\mathcal{A}}_q^{M,p}=\left\{\mathbf{\alpha }\in {\mathbb{N}}^M:||\mathbf{\alpha }|{|}_q\le p\right\},||\mathbf{\alpha }|{|}_q={\left(\sum _{i=1}^M{\alpha }_i^q\right)}^{1/q},$$

$$\hat{\mathbf{a}}=\arg \underset{\mathbf{a}\in {\mathbb{R}}^{|\mathcal{A}|}}{\min }\frac{1}{N}\sum _{i=1}^N{\left({\mathcal{Y}}^{(i)}-\sum _{\mathbf{\alpha }\in \mathcal{A}}{a}_{\mathbf{\alpha }}{\psi }_{\mathbf{\alpha }}\left({\mathbf{\chi }}^{(i)}\right)\right)}^2,$$

$$\hat{\mathbf{a}}={\left({\mathbf{F}}^{\text{T}}\mathbf{F}\right)}^{-1}{\mathbf{F}}^{\text{T}}\mathcal{Y}.$$

$$F_{ij}={\psi }_j\left({\mathbf{\chi }}^{(i)}\right),i=1,\dots ,N,j=1,\dots ,n_{\mathcal{A}}.$$

$$\begin{gathered} \mathcal{M}(\mathbf{x})={\mathcal{M}}_0+\sum _{i=1}^M{\mathcal{M}}_i(x_i)+\sum _{1\le i<j\le M}{\mathcal{M}}_{ij}(x_i,x_j)+\dots \\ +\sum _{1\le i_1<\dots <i_s\le M}{\mathcal{M}}_{i_1\dots i_s}\left(x_{i_1},\dots ,x_{i_s}\right)+\dots +{\mathcal{M}}_{1,2,\dots ,M}(\mathbf{x}), \end{gathered}$$

$${\int }_{{\mathcal{D}}_{X_k}}{\mathcal{M}}_{i_1\dots i_s}(x_{i_1},\dots ,x_{i_s})f_{X_k}\left(x_k\right)\text{d}{x}_k=0,1\le i_1<\dots <i_s\le M,k\in \left\{i_1,\dots ,i_s\right\}.$$

$${\mathcal{M}}_0={\int }_{{\mathcal{D}}_{\mathbf{X}}}\mathcal{M}(\mathbf{x})f_{\mathbf{X}}\left(\mathbf{x}\right)\text{d}\mathbf{x}.$$

$${\int }_{{\mathcal{D}}_{\mathbf{X}}}{\mathcal{M}}_{i_1\dots i_s}(x_{i_1},\dots ,x_{i_s}){\mathcal{M}}_{j_1\dots j_t}(x_{j_1,\dots ,j_t})f_{\mathbf{X}}\left(\mathbf{x}\right)\text{d}\mathbf{x}=0,\left\{i_1,\dots ,i_s\right\}\ne \left\{j_1,\dots ,j_t\right\}.$$

$$D=\text{Var}\left[\mathcal{M}(\mathbf{X})\right]={\int }_{{\mathcal{D}}_{\mathbf{X}}}{\mathcal{M}}^2(\mathbf{x})f_{\mathbf{X}}\left(\mathbf{x}\right)\text{d}\mathbf{x}-{\mathcal{M}}_0^2.$$

$$D=\sum _{i=1}^M{D}_i+\sum _{1\le i<j\le M}{D}_{ij}+\dots +D_{1,2,\dots ,M},$$

$$D_{i_1,\dots ,i_s}={\int }_{{\mathcal{D}}_{\mathbf{X}}}{\mathcal{M}}_{i_1\dots i_s}^2(x_{i_1},\dots ,x_{i_s})f_{\mathbf{X}}\left(\mathbf{x}\right)\text{d}\mathbf{x},1\le i_1<\dots <i_s\le M,s=1,\dots ,M.$$

$$S_{i_1\dots i_s}=\frac{D_{i_1\dots i_s}}{D}=\frac{D_{i_1\dots i_s}}{{\sum }_{i=1}^M{D}_i+{\sum }_{1\le i<j\le M}{D}_{ij}+\dots +D_{1,2,\dots ,M}},$$

$$\sum _{i=1}^M{S}_i+\sum _{1\le i<j\le M}{S}_{ij}+\dots +S_{1,2,\dots ,M}=1.$$

$$S_i^{(T)}=\sum _{\left\{i_1,\dots ,i_s\right\}\supset \left\{i\right\}}\frac{D_{i_1\dots i_s}}{D}.$$

$$\mathcal{I}_{i_{i},\ldots,i_{s}}=\left\{\boldsymbol{\alpha}\in\mathcal{A}:\begin{array}[]{l}\alpha_{k}>0\quad\text{for all}\quad k\in\left\{i_{1},\ldots,i_{s}\right\}\\ \alpha_{k}=0\quad\text{for all}\quad k\not\in\left\{i_{1},\ldots,i_{s}\right\}\end{array}\right\}.$$

$$\begin{gathered} {\mathcal{M}}^{(\text{P})}(\mathbf{x})=a_0+\sum _{i=1}^N\left(\sum _{\mathbf{\alpha }\in {\mathcal{I}}_i}{a}_{\mathbf{\alpha }}{\psi }_{\mathbf{\alpha }}(x_i)\right)+\sum _{1\le i_1<i_2\le M}\left(\sum _{\mathbf{\alpha }\in {\mathcal{I}}_{i_1,i_2}}{a}_{\mathbf{\alpha }}{\psi }_{\mathbf{\alpha }}(x_{i_1},x_{i_2})\right)+\dots \\ +\sum _{1\le i_1<\dots <i_s\le M}\left(\sum _{\mathbf{\alpha }\in {\mathcal{I}}_{i_1,\dots ,i_s}}{a}_{\mathbf{\alpha }}{\psi }_{\mathbf{\alpha }}(x_{i_1},\dots ,x_{i_s})\right)+\dots +\sum _{\mathbf{\alpha }\in {\mathcal{I}}_{1,2,\dots ,M}}{a}_{\mathbf{\alpha }}{\psi }_{\mathbf{\alpha }}(x_1,\dots ,x_M). \end{gathered}$$

$${\mathcal{M}}_{i_1\dots i_s}(x_{i_1},\dots ,x_{i_s})=\sum _{\mathbf{\alpha }\in {\mathcal{I}}_{i_1,\dots ,i_s}}{a}_{\mathbf{\alpha }}{\psi }_{\mathbf{\alpha }}(x_{i_1},\dots ,x_{i_s}).$$

$$D_{i_1\dots i_s}=\sum _{\mathbf{\alpha }\in {\mathcal{I}}_{i_1,\dots ,i_s}}{a}_{\mathbf{\alpha }}^2.$$

$$S_{i_1\dots i_s}^{(\text{P})}=\sum _{\mathbf{\alpha }\in {\mathcal{I}}_{i_1,\dots ,i_s}}{a}_{\mathbf{\alpha }}^2/\sum _{\mathbf{\alpha }\in \mathcal{A},\mathbf{\alpha }\ne 0}{a}_{\mathbf{\alpha }}^2.$$

$$S_i^{(T)(\text{P})}=\sum _{\left\{i_1,\dots ,i_s\right\}\supset \left\{i\right\}}{S}_{i_1,\dots ,i_s}^{(\text{P})}.$$

$${\underline{F}}_X^{(p)}(x)=\underset{\mathbf{\theta }\in {\mathcal{D}}_{\mathbf{\Theta }}}{\min }{F}_X(x|\mathbf{\theta }),{\overline{F}}_X^{(p)}(x)=\underset{\mathbf{\theta }\in {\mathcal{D}}_{\mathbf{\Theta }}}{\max }{F}_X(x|\mathbf{\theta }).$$

$${\underline{S}}_{i_1\dots i_s}=\underset{\mathbf{\theta }\in {\mathcal{D}}_{\mathbf{\Theta }}}{\min }{S}_{i_1\dots i_s}\left(\mathbf{\theta }\right),{\overline{S}}_{i_1\dots i_s}=\underset{\mathbf{\theta }\in {\mathcal{D}}_{\mathbf{\Theta }}}{\max }{S}_{i_1\dots i_s}\left(\mathbf{\theta }\right),$$

$${\underline{S}}_j^{(T)}=\underset{\mathbf{\theta }\in {\mathcal{D}}_{\mathbf{\Theta }}}{\min }{S}_j^{(T)}\left(\mathbf{\theta }\right),{\overline{S}}_j^{(T)}=\underset{\mathbf{\theta }\in {\mathcal{D}}_{\mathbf{\Theta }}}{\max }{S}_j^{(T)}\left(\mathbf{\theta }\right).$$

$$W={\mathcal{M}}^{(\text{aug})}\left(\mathbf{X},\mathbf{\Theta }\right),$$

Full text

Keywords: uncertainty quantification – - global sensitivity analysis – - probability-boxes – - imprecise Sobol’ indices – sparse polynomial chaos expansions

1 Introduction

Computational simulation tools, such as finite element models (FEM), are the most popular approach to model complex systems and processes in modern engineering. Such simulation tools map a set of input parameters describing the system and its environmental and operational conditions through a computational model to a so-called quantity of interest (QoI), *e.g. *performance indicators. The input parameters are often not perfectly known, due to noisy measurements, expert judgement, or intrinsic variability. Hence, each variable is modelled by an uncertainty model, as an example by a probability distribution. The uncertainty propagates through the computational model and results in an uncertain QoI.

Sensitivity analysis (SA) examines the impact of the uncertainty in the input vector of a computational model onto the uncertainty in the QoI. This is of importance in practice where the relation between the input variables and the QoI is implicitly defined by a complex computer code. A large number of SA methods can be found in the literature, as reported in the reviews in Saltelli et al. (2008); Helton et al. (2006); Xu and Gertner (2008). The methods can be separated broadly into two classes, namely local and global SA methods. Local SA examines how small variations in the input variables affect the QoI, whereas global SA focuses on the entire variation of the input variables. In this paper, we look at global SA methods.

Traditionally, global SA methods are developed in the context of probability theory, *i.e. *the uncertainty in the input variables is modelled by probability distributions. However, a common situation in practice is to have only scarce or incomplete knowledge to characterize the uncertainty. This leads to the concept of epistemic uncertainty (representing the lack of data and lack of knowledge) alongside with aleatory uncertainty (natural variability) as a source of uncertainty. For such cases, imprecise probability frameworks have been developed as a generalization of probability theory. Among these frameworks are Bayesian hierarchical models Gelman et al. (2009), probability-boxes Ferson and Ginzburg (1996), and fuzzy distributions Möller and Beer (2004).

The number of publications for imprecise sensitivity analysis (ISA) remains small. Li & Mahadevan Li and Mahadevan (2016b, a), Sankararaman & Mahadevan Sankararaman and Mahadevan (2013), and Krzykacz-Hausmann Krzykacz-Hausmann (2006) describe global SA methods in the presence of Bayesian hierarchical models. Pinching and imprecise Sobol’ indices are discussed in Oberguggenberger and Fellin (2005); Oberguggenberger et al. (2009). Pinching of probability-boxes is also elaborated in Ferson and Troy Tucker (2006). Helton et al. Helton et al. (2006) discuss variance-based algorithm in the context of evidence theory.

SA methods, and in particular ISA methods, usually require a large number of evaluations of the computational model for different realizations of the input vector. Hence, such analyses become intractable when the computational model is an expensive-to-evaluate function. Meta-models, a.k.a. surrogate models and emulators, are a popular solution to lower the computational costs by approximating the computational model by a cheap-to-evaluate analytical function. Examples include polynomial chaos expansions (PCE) (Xiu and Karniadakis, 2002; Ghanem and Spanos, 2003), Kriging (a.k.a. Gaussian process models) (Sacks et al., 1989; Santner et al., 2003), support vector machines (Gunn, 1998; Smola and Schölkopf, 2006), and artificial neural networks in the context of reliability analysis (Hurtado and Alvarez, 2001; Schueremans and Van Gemert, 2005).

A naive use of meta-models consists in the implementation of a sequential algorithm: first, calibrate a meta-model and subsequently conduct SA (or ISA) on the approximation. There are, however, more elaborate ways to use meta-models in the context of SA. As an example, based on polynomial chaos expansions (PCE) Ghanem and Spanos (2003), Sudret Sudret (2006, 2008) computes Sobol’ indices and Sudret & Mai Sudret and Mai (2015) compute derivative-based global sensitivity measures (DGSM) analytically, see also Alexanderian et al Alexanderian et al. (2012) and Le Gratiet et al Le Gratiet et al. (2017). Konakli & Sudret Konakli and Sudret (2016) also obtain analytical expressions for Sobol’ indices when the surrogate model is a low-rank tensor representation. However, the meta-model-based approaches are to-date limited to the probabilistic context.

In this paper, the standard Sobol’ indices Sobol’ (1993) are first extended to the context of probability-boxes in order to account for mixed epistemic and aleatory uncertainty in global sensitivity analysis. Moreover, making use of the PCE-based Sobol’ indices Sudret (2008) as well as a recently published variation of PCE Schöbi and Sudret (2015), a novel approach is proposed to compute imprecise Sobol’ indices efficiently. Applying the proposed approach, ISA are expected to become tractable for realistic engineering problems.

This paper is organized as follows. Section 2 recapitulates the basics of PCE and Sobol’ indices in the standard setting of pure probabilistic modelling of input variables. Section 3 defines probability-boxes and derives the definition of imprecise Sobol’ indices. In Section 4, a smart PCE model is introduced, which is later used to compute the imprecise Sobol’ indices in Section 5. Three applications illustrate the workings of the proposed approach in Section 6. The paper is terminated with conclusions in Section 7.

2 Polynomial chaos expansions for global sensitivity analysis

2.1 Computational model

Let us define a computational model $\mathcal{M}$ as a mapping from the $M$-dimensional input space $\mathcal{D}_{\boldsymbol{X}}$ to the one-dimensional output space $\mathcal{D}_{Y}$:

[TABLE]

Herein, we assume that the model predicts a scalar quantity of interest $y$ once it is fed with a vector of input parameters $\boldsymbol{x}$. Extensions to vector-valued models are straightforward by applying the proposed methodology to each component separately.

Due to uncertainties in the input vector, the latter is represented by a random vector $\boldsymbol{X}$ with joint cumulative distribution function (CDF) $F_{\boldsymbol{X}}$. The components of $\boldsymbol{X}=\left(X_{1},\ldots,X_{M}\right)^{\textsf{T}}$ are assumed independent for the sake of simplicity throughout this paper. The model response becomes a random variable $Y$ obtained by propagating the uncertainties in the input random vector $\boldsymbol{X}$ through the computational model $\mathcal{M}$.

Uncertainty propagation is typically conducted by sampling-based algorithms such as Monte Carlo simulation (Robert and Casella, 2004). When the computational model is an expensive-to-evaluate function, the estimation of $Y$ is costly due to the repeated evaluations of the computational model. Hence, a number of techniques have been proposed in the literature for surrogating $\mathcal{M}$. In the following section, Polynomial Chaos Expansions (Xiu and Karniadakis, 2002; Ghanem and Spanos, 2003; Soize and Ghanem, 2004) are introduced.

2.2 Sparse Polynomial Chaos Expansions

2.2.1 Polynomial approximation

Polynomial Chaos Expansions (PCE) approximate the computational model $\mathcal{M}$ by a series of multivariate polynomials that are orthogonal to the distributions of the input variables (Ghanem and Spanos, 2003; Soize and Ghanem, 2004):

[TABLE]

where $\psi_{\boldsymbol{\alpha}}\left(\boldsymbol{X}\right)$ are multivariate orthonormal polynomials indexed by multi-indices $\boldsymbol{\alpha}=\left(\alpha_{1},\ldots,\alpha_{M}\right)$, ${\mathcal{A}\subset\mathbb{N}^{M}}$ is a finite set of such indices (so-called truncation scheme), and $a_{\boldsymbol{\alpha}}\in\mathbb{R}$ are expansion coefficients. To derive the orthonormal polynomial basis, one proceeds as follows. Since the components of $\boldsymbol{X}$ are assumed independent, its joint probability density function (PDF) is the product of its marginals. Then, a functional inner product for each marginal PDF $f_{X_{i}}$ is defined by:

[TABLE]

for any two functions $\left\{\phi_{1},\phi_{2}\right\}$ such that the integral exists. For each variable $X_{i}$, an orthonormal polynomial basis can be constructed that satisfies (Xiu and Karniadakis, 2002):

[TABLE]

where $\left\{P_{j}^{(i)},P_{k}^{(i)}\right\}$ are two univariate polynomials in the $i^{\text{th}}$ variable and $\delta_{jk}$ is the Kronecker symbol which is equal to $1$ for $j=k$ and equal to [math] otherwise. Xiu & Karniadakis (Xiu and Karniadakis, 2002) summarize orthogonal bases for classical probability distributions.

The multivariate polynomials are then obtained by tensor product of the univariate ones:

[TABLE]

Through the tensor product construction, the orthogonality property is kept, meaning that:

[TABLE]

where $\delta_{\boldsymbol{\alpha}\boldsymbol{\beta}}$ is the Kronecker delta with $\delta_{\boldsymbol{\alpha}\boldsymbol{\beta}}=1$ for $\boldsymbol{\alpha}=\boldsymbol{\beta}$ and $\delta_{\boldsymbol{\alpha}\boldsymbol{\beta}}=0$ otherwise.

2.2.2 Truncation schemes

The efficiency of PCE models greatly depends on the choice of the set of multi-indices $\mathcal{A}$ in Eq. (2). A popular strategy to choose $\mathcal{A}$ consists in upper-bounding the total degree of polynomials $\sum_{i=1}^{M}\alpha_{i}$ to a maximum value $p$ or in using hyperbolic index sets (Blatman and Sudret, 2010):

[TABLE]

where $0<q\leq 1$ is a parameter and $p$ is the maximal total degree of the univariate polynomials retained.

2.2.3 Computation of coefficients $a_{\boldsymbol{\alpha}}$

After defining the set of candidate polynomials, the coefficients $a_{\boldsymbol{\alpha}}$ in Eq. (2) are to be determined. Based on a finite set of realizations of the input vector $\boldsymbol{X}$, denoted by $\boldsymbol{\mathcal{X}}=\left\{\boldsymbol{\chi}^{(1)},\ldots,\boldsymbol{\chi}^{(N)}\right\}$ and called experimental design, and the corresponding response values $\mathcal{Y}=\left\{\mathcal{Y}^{(1)},\ldots,\mathcal{Y}^{(N)}\right\}\stackrel{{\scriptstyle\text{def}}}{{=}}\left\{\mathcal{M}\left(\boldsymbol{\chi}^{(1)}\right),\ldots,\mathcal{M}\left(\boldsymbol{\chi}^{(N)}\right)\right\}$, the coefficients are calculated by discretized least-squares minimization Berveiller et al. (2006); Sudret (2007):

[TABLE]

where $n_{\mathcal{A}}=|\mathcal{A}|$ is the number of polynomials retained in the truncated expansion. The optimal coefficients $\widehat{\boldsymbol{a}}$ can be obtained by solving the linear system:

[TABLE]

In this expression, the information matrix ${\bf F}$ contains the values of all selected polynomials at all points of the experimental design, namely:

[TABLE]

2.2.4 Polynomial selection

Typically for smooth functions, a small number of polynomials is sufficient to accurately represent the output of the computational model. Thus, a further reduction of the set of predictors is possible. A number of methods are available in the literature, including the least absolute shrinkage and selection operator (LASSO) (Tibshirani, 1996), least angle regression (LAR) (Efron et al., 2004), orthogonal matching pursuit (OMP) (Pati et al., 1993; Mallat and Zhang, 1993), and Bayesian compressive sensing (Sargsyan et al., 2014). Based on LAR, Blatman & Sudret Blatman and Sudret (2011) proposed the LARS algorithm which allows for an efficient selection of a small number of polynomials out of the candidate set $\mathcal{A}_{q}^{M,p}$ (see also Eq. (7)). For a thorough discussion of the LARS algorithm, the reader is referred to Blatman and Sudret (2011); Blatman (2009).

2.3 Sobol’ indices

The so-called Sobol’ decomposition represents the computational model by a series of summands of increasing dimension (Sobol’, 1993):

[TABLE]

where $\mathcal{M}_{0}$ is a constant (mean value of the function) and where it is imposed that the integral of each summand $\mathcal{M}_{i_{1},\ldots,i_{s}}\left(x_{i_{1}},\ldots,x_{i_{s}}\right)$ over any of its arguments is zero:

[TABLE]

The mean value of the function can be computed as:

[TABLE]

Due to the constraints imposed in Eq. (12), the summands are orthogonal to each other so that (Homma and Saltelli, 1996):

[TABLE]

Considering now that the input parameters are modelled by independent random variables, the total variance of the computational model is defined as:

[TABLE]

By using the properties of Eq. (14) in Eq. (11), one obtains the following decomposition of the variance:

[TABLE]

where the partial variances are computed as:

[TABLE]

Making use of this decomposition, the Sobol’ indices are defined as the relative partial variances:

[TABLE]

and by virtue of Eq. (16), it holds that:

[TABLE]

The Sobol’ indices measure the amount of the total variance coming from the uncertainties in a set of input parameters. In practice, it is common to compute the first order indices $S_{i}$, which measure the influence of each parameter taken separately. Higher order indices account for the interactive contributions to the total variance. In this sense, the total sensitivity indices $S_{i}^{(T)}$ are defined as:

[TABLE]

In other words, $S_{i}^{(T)}=1-S_{-i}$, where $S_{-i}$ is the sum of all $S_{i_{1}\ldots i_{s}}$ that do not include index $i$.

2.4 PCE-based Sobol’ indices

Let us consider now that an appropriate approximation of the computational model $\mathcal{M}$ exists in the form of a truncated PCE, as seen in Eq. (2). Then, a set of multi-indices $\mathcal{I}_{i_{1},\ldots,i_{s}}$ can be defined so that only the indices $i_{1},\ldots,i_{s}$ are non-zero, i.e.:

[TABLE]

The elements of Eq. (2) can be reordered according to the decomposition in Eq. (11) (Sudret, 2008):

[TABLE]

Then, the summands of Eq. (22) can be identified as summands of the Sobol’ decomposition:

[TABLE]

Due to the property of uniqueness, it can be concluded that Eq. (22) is indeed the Sobol’ decomposition of the PCE model. Furthermore, due to the orthonormality of the PC basis, the partial variance $D_{i_{1}\ldots i_{s}}$ is readily obtained from Eq. (23) as:

[TABLE]

The corresponding polynomial-chaos-based Sobol’ indices, denoted by $S^{(\text{P})}_{i_{1}\ldots i_{s}}$, are then defined, as originally proposed in Sudret (2006), as:

[TABLE]

The total PC-based sensitivity indices are then computed analogously to Eq. (20):

[TABLE]

As seen in this section, the computation of the Sobol’ indices for a PCE model boils down to working with the coefficients of the PC-expansion, which is a simple operation at almost no additional computational cost. Hence, PCE-based Sobol’ indices are obtained by simple post-processing of an existing PCE model.

3 Imprecise Sobol’ indices

3.1 Free and parametric p-boxes

Sobol’ indices can be computed by the procedure described in the above section as soon as the input parameters are modelled by well-defined probability density functions. In practice however, probability theory is often not fully appropriate to characterize the uncertainty in the system. Probability-boxes (p-boxes) pose an intuitive generalization of probability theory. Instead of representing $X$ by a single CDF, lower and upper bounds denoted by $\underline{F}_{X}$ and $\overline{F}_{X}$ are introduced (Ferson and Ginzburg, 1996; Ferson and Hajagos, 2004). For any value of $x\in\mathcal{D}_{X}$, the true-but-unknown CDF lies within these bounds such that $\underline{F}_{X}(x)\leq F_{X}(x)\leq\overline{F}_{X}(x)$, which form a so-called p-box.

In the literature, two types of p-boxes can be found, namely free and parametric p-boxes. Free p-boxes are defined directly though the formulation of the bounds $\left\{\underline{F}_{X},\overline{F}_{X}\right\}$. Parametric p-boxes (a.k.a distributional p-boxes) are defined as cumulative distribution functions with interval-valued hyper-parameters: $F_{X}(x)=F_{X}(x|\boldsymbol{\theta})$, where $\boldsymbol{\theta}\in\mathcal{D}_{\boldsymbol{\Theta}}\subset\mathbb{R}^{n_{\boldsymbol{\theta}}}$, and $\mathcal{D}_{\boldsymbol{\Theta}}$ is the interval domain of the distribution parameters of dimension $n_{\boldsymbol{\theta}}=|\boldsymbol{\theta}|$. Precisely, $\mathcal{D}_{\boldsymbol{\Theta}}=\left[\underline{\theta}_{1},\overline{\theta}_{1}\right]\times\ldots\times\left[\underline{\theta}_{n_{\boldsymbol{\theta}}},\overline{\theta}_{n_{\boldsymbol{\theta}}}\right]$ denotes a hyper-rectangular domain, and $\underline{\theta}_{k}$ and $\overline{\theta}_{k}$ denote the lower and upper bounds of the interval $\left[\underline{\theta}_{k},\overline{\theta}_{k}\right]$ for the $k$-th parameter of the CDF of $X$. Boundary curves for parametric p-boxes can be defined as for free p-boxes as follows:

[TABLE]

Figure 1 shows the intrinsic difference between free and parametric p-boxes. In both cases, only non-decreasing realizations of the true-but-unknown CDF are feasible. However, free p-boxes are a more general description of uncertainty than parametric p-boxes. This can be seen by the stair-case realizations, which are possible in free but not in parametric p-boxes (see also Realization #3 in Figure 1(A)).

In this paper, the focus lies on the parametric p-boxes due to their clear separation of aleatory and epistemic uncertainty: the distribution family represents the aleatory uncertainty, whereas the interval-valued hyper-parameters represent the epistemic uncertainty.

3.2 Imprecise Sobol’ indices

In the context of parametric p-boxes, uncertainty quantification is typically conducted by nested Monte Carlo simulations (nMCS) algorithms where an outer loop samples the parameters $\boldsymbol{\theta}$ and an inner loop samples the conditional probability distribution $F_{\boldsymbol{X}}\left(\boldsymbol{x}|\boldsymbol{\theta}\right)$ (Eldred and Swiler, 2009; Chowdhary and Dupuis, 2013). Taking this idea into the context of sensitivity analysis, the Sobol’ indices can be computed for each realization $\boldsymbol{\theta}^{(j)}$ and the corresponding $F_{\boldsymbol{X}}\left(\boldsymbol{x}\middle|\boldsymbol{\theta}^{(j)}\right)$. Each realization can result in a different set of Sobol’ indices. Finally, the bounds of the imprecise Sobol’ indices can be determined by analysing the set of realizations.

In the sequel, we define imprecise Sobol’ indices $S_{i_{1}\ldots i_{s}}$ as a pair of bounds obtained by minimizing (resp. maximizing) the classical Sobol’ indices over the range of the hyper-parameters $\boldsymbol{\theta}\in\mathcal{D}_{\boldsymbol{\Theta}}$:

[TABLE]

where $S_{i_{1}\ldots i_{s}}\in\left[\underline{S}_{i_{1}\ldots i_{s}},\,\overline{S}_{i_{1}\ldots i_{s}}\right]$ are the interval-valued Sobol’ indices and $S_{i_{1}\ldots i_{s}}\left(\boldsymbol{\theta}\right)$ is the corresponding Sobol’ index for the conditional input distribution $F_{\boldsymbol{X}}\left(\boldsymbol{x}|\boldsymbol{\theta}\right)$. Similarly to Eq. (26), the bounds of the total Sobol’ indices can be obtained by:

[TABLE]

3.3 PCE-based indices

From the definitions above, it can be seen that computing the two bounds of imprecise Sobol’ indices requires solving two optimization problems over the hyper-parameter space $\mathcal{D}_{\boldsymbol{\Theta}}$. In order to speed up the process, we propose to use polynomial chaos expansions as presented in Section 2.2.

The use of the PCE meta-models would reduce the computational effort of computing Sobol’ indices conditioned on a particular $\boldsymbol{\theta}$ but does not remove the repeated computation of the variance decomposition for many $\boldsymbol{\theta}\in\mathcal{D}_{\boldsymbol{\Theta}}$ to solve the optimization problems. In other words, when adopting a brute force procedure, a large number of PCE models should be calibrated, each of them with its own experimental design, own evaluations of the costly computational model $\mathcal{M}$, etc. Hence in the following, we introduce a smart PCE based on a single experimental design that is sufficient to estimate any $S_{i}\left(\boldsymbol{\theta}\right)$ where $\boldsymbol{\theta}\in\mathcal{D}_{\boldsymbol{\Theta}}$.

4 Augmented polynomial chaos expansions for parametric distributions

4.1 Augmented space

Due to the clear separation of aleatory and epistemic uncertainty in the formulation of parametric p-boxes, these sources of uncertainty can be treated as separate entities, *i.e. *$\boldsymbol{X}$ and $\boldsymbol{\Theta}$. Hence, the QoI can be interpreted as a function of the augmented input vector $\left(\boldsymbol{X},\boldsymbol{\Theta}\right)$. The corresponding map then reads:

[TABLE]

where the augmented computational model is based on the original computational model as $\mathcal{M}^{(\text{aug})}\left(\boldsymbol{x},\boldsymbol{\theta}\right)\equiv\mathcal{M}\left(\boldsymbol{x}\right)$.

The components of the augmented input vector are dependent on each other because $\boldsymbol{X}$ depends on $\boldsymbol{\Theta}$ as a result of the hierarchical formulation of the parametric p-box. The dependencies are visualized in Figure 2(A). In this paper, however, the input variables to the computational model are assumed to be independent. In order to obtain independent variables in the augmented input vector, an isoprobabilistic transform is used that maps $\left(\boldsymbol{X},\boldsymbol{\Theta}\right)$ to $\boldsymbol{V}$, which shall be a vector with independent components. Assuming that the components $X_{i}$ are independent and that the parameter intervals in $\boldsymbol{\Theta}$ are independent too (*i.e. *$\mathcal{D}_{\boldsymbol{\Theta}}$ is a hyper-rectangular domain), the parametric p-box $\boldsymbol{X}$ can be recast as:

[TABLE]

where $\boldsymbol{C}=\left(C_{1},\ldots,C_{M}\right)$ is a vector of independent uniform distributions and where $C_{i}\sim\mathcal{U}(0,1)$ describes the CDF value of $X_{i}$. Then, the augmented vector $\boldsymbol{V}=\left(\boldsymbol{C},\boldsymbol{\Theta}\right)$ is a vector with independent components that allows to define explicitly any realization of $\boldsymbol{X}$ through Eq. (31).

The corresponding probabilistic graphical model transform is visualized in Figure 2(B). Note that $\boldsymbol{x}$ is used in lower case due to the deterministic transform of $\boldsymbol{V}=\left(\boldsymbol{C},\boldsymbol{\Theta}\right)\mapsto\boldsymbol{X}$. In Figure 2(A), the transform of $\boldsymbol{\Theta}$ to $\boldsymbol{X}$ is not deterministic, hence, the upper case $\boldsymbol{X}$.

4.2 Augmented PCE

The augmented computational model $\mathcal{M}^{(\text{aug})}$ can be formulated as a function of vector $\boldsymbol{V}$ as follows:

[TABLE]

where $T$ denotes the isoprobabilistic transform from $\boldsymbol{V}$ to $\left(\boldsymbol{X},\boldsymbol{\Theta}\right)$. When interpreting the epistemic intervals in $\boldsymbol{\Theta}$ as uniform distributions, a sparse PCE model can be calibrated to meta-model the augmented computational model (analogously to Eq. (2)):

[TABLE]

where the dimensionality is $M=|\boldsymbol{V}|=|\boldsymbol{C}|+|\boldsymbol{\Theta}|\equiv|\boldsymbol{X}|+|\boldsymbol{\Theta}|$. The meta-model should be trained with an experimental design $\boldsymbol{\mathfrak{V}}=\left\{\boldsymbol{\mathfrak{v}}^{(1)},\ldots,\boldsymbol{\mathfrak{v}}^{(N)}\right\}$ and the corresponding responses:

[TABLE]

by using Eq. (32). In order to efficiently calibrate a sparse PCE meta-model, the reader is referred to Section 2.2 for further details. Note that in this setup, the input vector $\boldsymbol{V}$ is sampled in the same way for the $\boldsymbol{C}$- and the $\boldsymbol{\Theta}$-components.

4.3 Phantom points

4.3.1 Isoprobabilistic transform

The number of dimensions in the augmented input space $\mathcal{D}_{\boldsymbol{V}}$ is larger than the one of the original input space $\mathcal{D}_{\boldsymbol{X}}$ due to the epistemic uncertainty. However, Eq. (31) leads to an interesting feature of the augmented PCE. For a given realization of $\boldsymbol{x}\in\mathcal{D}_{\boldsymbol{X}}$, there is an infinite number of realizations $\boldsymbol{v}$ that satisfy Eq. (31). In other words, the inverse operation of Eq. (31) is non-unique. This characteristics can be exploited when constructing the experimental design of the augmented PCE model.

Consider again the experimental design $\boldsymbol{\mathfrak{V}}$. Each sample $\boldsymbol{\mathfrak{v}}^{(j)}$ can be transformed into an equivalent $\boldsymbol{\chi}^{(j_{0})}$ using Eq. (31). A number of realizations $\boldsymbol{\mathfrak{v}}=\left(\boldsymbol{\mathfrak{c}},\boldsymbol{\tau}\right)$ can be generated by sampling $\boldsymbol{\tau}\in\mathcal{D}_{\boldsymbol{\Theta}}$ and inserting them into the inverse isoprobabilistic transform to obtain realizations $\boldsymbol{\mathfrak{c}}$ that all correspond to the single point $\boldsymbol{\chi}^{(j_{0})}$ in the original experimental design. In particular, for each dimension $i=1,\ldots,M$ and experimental design sample $j=1,\ldots,N$, $n_{ph}$ samples can be generated from $\boldsymbol{\theta}_{i}\in\mathcal{D}_{\boldsymbol{\Theta}_{i}}$ and collected in $\left\{\boldsymbol{\tau}_{i}^{(j)(k)},k=1,\ldots,n_{ph}\right\}$. Then, the inverse transform consists of evaluating simple CDFs and reads:

[TABLE]

where each experimental design sample point is $\boldsymbol{\chi}^{(j)}=\left(\chi_{1}^{(j)},\ldots,\chi_{M}^{(j)}\right)$. The new components in dimension $i$ of vector $\boldsymbol{\mathfrak{v}}^{(j)}$ reads then:

[TABLE]

When pursuing this procedure for all $i=1,\ldots,M$ dimensions of the input vector, new samples of the experimental design are generated. These samples are called phantom points (Schöbi and Sudret, 2015) because they contribute to the experimental design $\boldsymbol{\mathfrak{V}}$ without increasing the number of model evaluations $\mathcal{M}$, which is once and for all $N$. Indeed, the $n_{ph}$ points in Eq. (35) all correspond to a single $\boldsymbol{\chi}^{(j_{0})}$ in the original experimental design, meaning a single run $\mathcal{Y}^{(j_{0})}=\mathcal{M}\left(\boldsymbol{\chi}^{(j_{0})}\right)$. This has a large effect on the global efficiency when considering that model evaluations are dominating the total computational costs. The experimental design with phantom points and the corresponding responses then read:

[TABLE]

which consists of $n_{ph}\times N$ samples at a cost of exactly $N$ evaluations of the computational model $\mathcal{M}$.

For illustration purposes, consider a parametric p-box defined by a Gaussian distribution family with interval-valued mean and standard deviation. As an example, it is assumed that ${\mu_{X}\in\left[-1,\,1\right]}$ and ${\sigma_{X}\in\left[0.5,\,1.0\right]}$. Figure 3 shows the bounds of the parametric p-box (see Eq. (27)) as well as the experimental sample point $\chi^{(0)}=0$. Furthermore, the dotted and dashed-dotted blue lines show two realizations of the CDF in the parametric p-box, namely $\left\{\mu_{X}=-0.5,\,\sigma_{X}=1\right\}$ and $\left\{\mu_{X}=0.5,\,\sigma_{X}=0.5\right\}$, respectively. The phantom points corresponding to $\chi^{(0)}$ are then computed through the CDF of the Gaussian distribution:

[TABLE]

which are highlighted by the red diamonds in Figure 3.

As a summary, a parametric p-box $F_{x}\left(x\middle|\boldsymbol{\theta};\,\boldsymbol{\theta}\in\mathcal{D}_{\boldsymbol{\Theta}}\right)$ can be equivalently represented by $\left(n_{\boldsymbol{\theta}}+1\right)$ independent random variables in the augmented space. The extended model $\mathcal{M}^{(\text{aug})}$ in the augmented space is formulated as a function of a $n_{\text{aug}}$-dimensional input vector, where $n_{\text{aug}}=M+\sum_{i=1}^{M}n_{\boldsymbol{\theta}}^{(i)}$, $n_{\boldsymbol{\theta}}^{(i)}$ being the number of parameters $\boldsymbol{\theta}$ characterizing the $i^{\text{th}}$ p-box. However, using the “phantom points” trick, an experimental design of arbitrary size $n_{ph}\times N$ can be computed in the augmented input space at the cost of only $N$ evaluations of the original model $\mathcal{M}$. This leads to a remarkably efficient computational procedure. Further computational details are given below.

4.3.2 Case of unbounded distributions

In the case of unbounded distributions, the phantom points can be generated as described above. However, the CDF used in Eq. (34) is a highly non-linear function due to its mapping of an unbounded to a bounded domain, *i.e. *$X\in\left[-\infty,\,\infty\right]\mapsto C\in\left[0,\,1\right]$. In order to avoid the non-linearity in this mapping and hence expecting more accurate meta-models when using low-degree polynomials, an auxiliary variable $\xi$ is introduced replacing the uniform distribution $C$. Then, an auxiliary isoprobabilistic transform is formulated for the mapping $\left(\xi,\boldsymbol{\Theta}\right)\mapsto X$, which shall be less non-linear. In the following, such an auxiliary transform is derived for the case of Gaussian, lognormal, Gumbel, and Weibull parametric p-boxes.

Gaussian p-box. Consider a Gaussian parametric p-box where the mean value $\mu_{X}$ and standard deviation $\sigma_{X}$ are given in intervals. Instead of using the uniform distribution $C$, a standard Gaussian variable $\xi\sim\mathcal{N}(0,\,1)$ can be used in the independent vector $\boldsymbol{V}=\left(\xi,\mu_{X},\sigma_{X}\right)$. Then, the isoprobabilistic transform reads:

[TABLE]

This transform is simpler than the inverse CDF formulation in Eq. (31) because the unbounded variable $X$ is modelled by means of the unbounded variable $\xi$. In fact, the transform features a linear mapping from $X$ to $\xi$.

Lognormal p-box. Considering a lognormal parametric p-box with interval-valued parameters $\left(\lambda_{X},\zeta_{X}\right)$, such that $\log X\sim\mathcal{N}\left(\lambda_{X},\zeta_{X}\right)$, then:

[TABLE]

where $\xi$ follows a standard Gaussian distribution, and $\left(\lambda_{X},\,\zeta_{X}\right)$ are obtained by the usual equations from the mean value $\mu_{X}$ and standard deviation $\sigma_{X}$. Note that a lognormal p-box defined by lower and upper bounds on mean $\mu_{X}$ and standard deviation $\sigma_{X}$ leads to a non-rectangular domain in the $\left(\lambda_{X},\,\zeta_{X}\right)$-plane. This is, however, not an issue since the experimental designs will be sampled in the $\left(\mu_{X},\,\sigma_{X}\right)$-space.

Gumbel p-box. The CDF of a Gumbel distribution is defined as:

[TABLE]

where $\alpha=\mu_{X}-\beta\gamma_{e}$ and $\beta=\sigma_{X}\sqrt{6}/\pi$ are the distribution hyper-parameters ($\mu_{X}$ and $\sigma_{X}$ denote the mean and standard deviation) and $\gamma_{e}=0.5772\ldots$ is the Euler-Mascheroni constant. Then, a standard Gumbel variable is denoted by $\varpi\sim\mathcal{GU}(\alpha=0,\beta=1)$. It follows that a distribution with arbitrary hyper-parameter values $\left(\alpha,\beta\right)$ can be formulated in terms of a standard Gumbel distribution as follows:

[TABLE]

which reduces to:

[TABLE]

Hence, a standard Gumbel distribution can be transformed into an arbitrary Gumbel distribution by:

[TABLE]

which resembles the construction for Gaussian variables in Eq. (39) with different transformation parameters. Accordingly, a Gumbel p-box can be represented by three random variables, including two uniform and one standard Gumbel.

Weibull p-box. The CDF of a Weibull distribution is defined as:

[TABLE]

where $x\in\left[0,\,\infty\right]$ and $\left(\alpha,\beta\right)$ are its hyper-parameters. Then, a standard Weibull distribution is denoted by $\varpi\sim\mathcal{W}(\alpha=1,\beta=1)$, which also corresponds to a standardized exponential distribution. It follows that a distribution with arbitrary hyper-parameter values $\left(\alpha,\beta\right)$ can be formulated in terms of an exponential distribution as follows:

[TABLE]

Accordingly, a Weibull p-box can be represented by three random variables including two uniform and one exponential.

4.3.3 Bounded distributions

When dealing with distributions that are bounded on two sides, *i.e. *$X\in\left[\underline{x},\,\overline{x}\right]$, the phantom points cannot be generated by the approach detailed in the previous section. The problem is illustrated in Figure 4, which shows a parametric p-box consisting of a uniform distribution with interval-valued distribution bounds, *i.e. *$X\sim\mathcal{U}\left(a,b\right)$, where $a\in[1,\,2]$ and $b\in[3,\,4]$. Consider the experimental design point $\chi^{(0)}=3.5$ and the two realizations with $\left\{a=1.2,\,b=3.8\right\}$ and $\left\{a=1,\,b=3.25\right\}$. For the first realization, the phantom point for $\chi^{(0)}$ reads:

[TABLE]

For the second realization, however, the phantom point lies outside of the support interval $X\in[1,\,3.25]$. Hence, the construction of a phantom point is not meaningful in this case. For $\chi^{(0)}=3.5$ in particular, a phantom point is only feasible when $a\in\left[1,\,2\right]$ and $b\geq 3.5$ at the same time.

In other words, the generation of phantom points is submitted to constraints in case of bounded parametric distributions. In order to generate phantom points efficiently, these constraints must be taken into account. Therefore, Eq. (34) is transformed into:

[TABLE]

where $X_{\boldsymbol{\tau}_{i}^{(j)(k)}}$ denotes a conditional distribution defined within the bounds $X_{\boldsymbol{\tau}_{i}^{(j)(k)}}\in\left[\underline{x}_{\boldsymbol{\tau}_{i}^{(j)(k)}},\ \overline{x}_{\boldsymbol{\tau}_{i}^{(j)(k)}}\right]$. When $\chi_{i}^{(j)}$ is outside the interval bounds, no phantom point is generated for the realization $\boldsymbol{\tau}_{i}^{(j)(k)}$.

5 Imprecise Sobol’ indices by PCE

5.1 Reordering

When an augmented PCE model is available as in Eq. (33), the imprecise Sobol’ indices can be obtained by post-processing this augmented PCE model. The imprecise Sobol’ indices shall be depending on $\boldsymbol{C}$ and conditioned on $\boldsymbol{\Theta}$. In order to do so, the augmented PCE model is recast as follows:

[TABLE]

where $\boldsymbol{\alpha}=\left(\boldsymbol{\alpha}_{\boldsymbol{C}},\boldsymbol{\alpha}_{\boldsymbol{\Theta}}\right)$ forms the same index set as in Eq. (33) and $\psi_{\boldsymbol{\alpha}}\left(\boldsymbol{V}\right)=\psi_{\boldsymbol{\alpha}_{\boldsymbol{C}}}\left(\boldsymbol{C}\right)\cdot\psi_{\boldsymbol{\alpha}_{\boldsymbol{\Theta}}}\left(\boldsymbol{\Theta}\right)$ are the corresponding multivariate orthonormal polynomials. With this notation, we emphasize that any polynomial of $\left(\boldsymbol{C},\boldsymbol{\Theta}\right)$ can be cast as a product of two polynomials of $\boldsymbol{C}$ and $\boldsymbol{\Theta}$, respectively, due to the initial tensor product construction. Then, the computational model can be further rearranged, for a given value $\boldsymbol{\theta}$ as a function of $\boldsymbol{C}$ only:

[TABLE]

where $a_{\boldsymbol{\alpha},\boldsymbol{\theta}}(\boldsymbol{\theta})=a_{\boldsymbol{\alpha}}\ \psi_{\boldsymbol{\alpha}_{\boldsymbol{\Theta}}}\left(\boldsymbol{\theta}\right)$ is a new coefficient dependent on the value of $\boldsymbol{\theta}$. Note that several $\boldsymbol{\alpha}_{\boldsymbol{\Theta}}$ may correspond to the same $\boldsymbol{\alpha}_{\boldsymbol{C}}$ in $\boldsymbol{\alpha}$. Hence, Eq. (49) can be further rewritten to a model with proper variance decomposition as in Eq. (22):

[TABLE]

where $i$ marks the $1,\ldots,n_{\mathcal{A}}$ multi-indices in $\mathcal{A}$, $\mathbb{I}$ is the indicator function with $\mathbb{I}=1$ for $\boldsymbol{\alpha}_{\boldsymbol{C}}^{(i)}=\boldsymbol{\alpha}_{\boldsymbol{C}}^{*}$ and $\mathbb{I}=0$ otherwise, $\mathcal{A}_{\boldsymbol{C}}^{*}$ is the set of unique multi-indices $\boldsymbol{\alpha}_{\boldsymbol{C}}$ in $\mathcal{A}$, and $\boldsymbol{\alpha}^{(i)}=\left(\boldsymbol{\alpha}_{\boldsymbol{C}}^{(i)},\boldsymbol{\alpha}_{\boldsymbol{\Theta}}^{(i)}\right)$. Rearranging the terms to obtain the same structure as in Eq. (22) leads to the following coefficients depending on $\boldsymbol{\theta}$:

[TABLE]

The simplified decomposition then reads:

[TABLE]

5.2 Optimization

Finally, in order to obtain the extreme values of the Sobol’ indices, the optimizations in Eq. (28) use the reformulated model in Eq. (52). Then, the bounds of imprecise Sobol’ indices can be computed by solving two optimization problems:

[TABLE]

[TABLE]

In a nutshell, once the augmented PCE has been calibrated, the conditional Sobol’ indices $S_{i_{1}\ldots i_{s}}\left(\boldsymbol{\theta}\right)$ are simply ratios of multivariate orthogonal polynomials in $\boldsymbol{\theta}$. The two optimization problems in Eqs. (53) and (54) can then be solved at reasonable cost through global optimization algorithms, such as genetic or simulated annealing algorithms.

5.3 Connection to Bayesian hierarchical models

The imprecise Sobol’ indices computed here for parametric p-boxes can be extended to Bayesian hierarchical models with little modifications. When interpreting $\boldsymbol{\theta}$ as distributions rather than interval-valued quantities, Eqs. (53) and (54) converge to a single equation describing a distribution of Sobol’ indices:

[TABLE]

As a direct consequence, the distribution of any Sobol’ index $S_{i_{1}\ldots i_{s}}^{(\text{P})}$ due to the epistemic uncertainty in $\boldsymbol{\theta}$ could be obtained at low cost by Monte Carlo sampling of Eq. (55).

6 Applications

6.1 Step-by-step illustration of algorithm

6.1.1 Problem statement

A simple two-dimensional function is considered here to visualize the algorithm for computing the imprecise Sobol’ indices. The analytical function is defined as:

[TABLE]

where $x_{1}$ and $x_{2}$ are modelled by parametric p-boxes. The input variables are modelled by Gaussian distributions with interval-valued mean and standard deviation, namely ${\mu_{i}\in\left[-1,1\right]}$ and ${\sigma_{i}\in\left[0.5,1.0\right]}$. As seen in the sequel, this function allows us to derive analytically the conditional Sobol’ indices as well as their bounds.

6.1.2 Imprecise Sobol’ indices

Each input parameter may be cast as $X_{i}=\mu_{i}+\sigma_{i}\cdot\xi_{i}$ where $\xi_{i}$ is a standard normal variable. Then, the augmented computational model can be written as:

[TABLE]

The variance of the response variable conditioned on the values of $\boldsymbol{\mu}$ and $\boldsymbol{\sigma}$ then reads:

[TABLE]

Hence, the first order Sobol’ indices can be computed as a function of $\boldsymbol{\theta}=\left(\boldsymbol{\mu},\boldsymbol{\sigma}\right)$:

[TABLE]

The analytical derivation of the extreme Sobol’ indices results in the following bounds for the first order Sobol’ indices:

[TABLE]

As an example, the minimal first order Sobol’ index for $X_{1}$ is obtained by setting $\mu_{2}=0$. Similarly to the first order indices, the total order Sobol’ indices can be written as:

[TABLE]

The bounds of the total order Sobol’ indices read:

[TABLE]

As an example, the minimal total order Sobol’ index for $X_{1}$ is obtained by setting $\mu_{1}=\pm 1$, $\sigma_{1}=0.5$, $\mu_{2}=0$, and $\sigma_{2}=1$. Due to the symmetry of the problem definition, *i.e. *computational model $f_{1}^{(\text{aug})}$ and input p-boxes $X_{i}$, the Sobol’ indices are symmetrical in $\left(\mu_{i},\sigma_{i}\right)$. Note that in this example, the first and total order Sobol’ indices are not identical due to the interactive term $\left(\sigma_{1}\cdot\sigma_{2}\right)\cdot\xi_{1}\cdot\xi_{2}$ in Eq. (57).

6.1.3 Augmented PCE-based Sobol’ indices

Using the Matlab-based uncertainty quantification framework UQLab (Marelli and Sudret, 2014, 2015), the augmented PCE of $f_{1}^{(\text{aug})}$ and the parametric p-boxes results in the following set of multi-indices $\mathcal{A}$:

[TABLE]

which contains ten vectors $\boldsymbol{\alpha}$ and where the indices in each row correspond to the augmented input vector $\left(\mu_{1},\sigma_{1},\xi_{1},\mu_{2},\sigma_{2},\xi_{2}\right)$. Splitting the multi-index set $\mathcal{A}$ to the contributions of the epistemic ($\mu_{1},\,\sigma_{1},\,\mu_{2},\,\sigma_{2}$) and aleatory ($\xi_{1},\,\xi_{2}$) variables leads to:

[TABLE]

Note that according to Eq. (57), only four elements should be part of $\mathcal{A}$. However, due to the normalization to a standard space (uniform distributions to $[-1,1]$ and Gaussian to standard normal space), an additional six elements are present.

$\mathcal{A}_{\boldsymbol{C}}^{*}$ summarizes the set of unique row vectors $\boldsymbol{\alpha}_{\boldsymbol{C}}$ in $\mathcal{A}_{\boldsymbol{C}}$ and contains the following four elements:

[TABLE]

The final coefficients are composed as follows (see also Eq. (51)):

[TABLE]

where Eqs. (66) and (67) can be simplified to:

[TABLE]

The variance decomposition of the PCE model then reads (see also Eq. (52)):

[TABLE]

where $\boldsymbol{C}=\left(\xi_{1},\xi_{2}\right)$. By optimization on $\boldsymbol{\theta}$, the PC-based imprecise Sobol’ indices can be computed. In this particular case, we obtain the analytical results of Eqs. (60) and (62) exactly. Indeed, the initial model $f_{1}$ is polynomial, as well as the isoprobabilistic transforms. In the end, a 10-term PCE in the augmented space is exact.

6.2 Single-degree-of-freedom oscillator

6.2.1 Problem statement

Consider the non-linear undamped single-degree-of-freedom (SDOF) oscillator sketched in Figure 5 (Schueremans and Van Gemert, 2005; Echard et al., 2011, 2013). The corresponding computational model reads:

[TABLE]

where $m$ is the mass, $\left\{c_{1},c_{2}\right\}$ are the spring constants of the primary and secondary springs, $r$ is the displacement at which the secondary spring yields, $t_{1}$ is the duration of the loading, $F_{1}$ is the amplitude of the force and $\omega_{0}=\sqrt{\frac{c_{1}+c_{2}}{m}}$ is the natural frequency of the oscillator.

The input vector is modelled by a mix of probabilistic variables and parametric p-boxes accounting for the different levels of knowledge. The description of the input variables is provided in Table 1. It is assumed that the spring stiffnesses and the mass are well-known. Hence $\left\{c_{1},c_{2},m\right\}$ are modelled by regular (“precise”) CDFs. On the other side, knowledge on $\left\{r,F_{1},t_{1}\right\}$ is supposed to be scarce. Hence, these variables are modelled by parametric p-boxes. As seen in Table 1, the parametric p-box is characterized by a distribution function with interval-valued mean value but constant standard deviation.

6.2.2 Analysis

In order to estimate the imprecise Sobol’ indices, an augmented PCE model is trained for an augmented input space of dimension $M_{\text{aug}}=6+3=9$. The number of model evaluations $N$ is varied to investigate the influence of the experimental design size. Note that the number of phantom points is set to $n_{ph}=10$ for any $N$ in order to take full advantage of the phantom points. The sparse PCE is built with hyperbolic index sets ($q=0.75$) and degree-adaptive LARS with a maximal total polynomial degree $p=10$. The bounds of the imprecise Sobol’ indices are then obtained through Matlab’s genetic optimization algorithm.

A reference solution is obtained by using augmented PCE with a large sample size, *i.e. *$N=1000$ and $n_{ph}=10$. This results in an accurate meta-model and hence an accurate estimate of the imprecise Sobol’ indices.

6.2.3 Results

The estimate for the first-order Sobol’ indices $S_{i}$ are summarized in Table 2. The imprecise Sobol’ indices are shown for the cases of $N=\left\{30,\,50,\,100,\,200\right\}$. The larger the experimental design, the closer the estimates of bounds of the Sobol’ indices to the reference values. In fact, an experimental design of $N=50$ is sufficient to obtain results with sufficient accuracy in this example.

In order to estimate the accuracy of the augmented PCE models, the relative generalization error is derived as the mean squared error between prediction and exact response value relative to the variance of the response variable. Based on a large validation set of $n=10^{6}$ samples generated from the augmented input vector $\boldsymbol{V}$, the error estimate reads:

[TABLE]

where $\mathbb{W}=\left\{w_{1},\ldots,w_{n}\right\}$ are the exact response values corresponding to the validation set $\mathbb{V}=\left\{\boldsymbol{v}_{1},\ldots,\boldsymbol{v}_{n}\right\}$, $w_{i}^{(\text{P})}$ are the corresponding values predicted with the augmented PCE model, and $\mu_{w}=\frac{1}{n}\sum_{i=1}^{n}w_{i}$ is the mean value of the samples in $\mathbb{W}$.

The values of the relative generalization error are listed in Table 2. As indicated for the bounds of the imprecise Sobol’ indices, the augmented PCE model is accurate for large experimental designs. Values smaller than $\widehat{err}_{gen}\approx 10^{-4}$ indicate a high prediction accuracy.

Figure 6(A) illustrates the imprecise Sobol’ indices obtained by $N=200$ and $n_{ph}=10$ in a bar plot. The solid (black) part of the bars represents the lower bound of the first order Sobol’ indices whereas the hollow bar represents the upper bound. For comparison to a purely probabilistic approach (*i.e. *no epistemic uncertainty), the parametric p-boxes are “pinched” into classical probability distributions. The CDFs use the central value of the intervals as hyper-parameter: $r\sim\mathcal{N}\left(0.50,\,0.05\right)$, $F_{1}\sim\mathcal{N}\left(1.0,\,0.2\right)$, and $t_{1}\sim\mathcal{N}\left(1.00,\,0.2\right)$. The resulting Sobol’ indices are shown by blue markers in the bars. As expected, the imprecise Sobol’ indices envelope the probabilistic estimates due to the fact that the “probabilistic case” is one possible realization of the parametric p-boxes.

A different visualization of the imprecise Sobol’ indices is shown in Figure 6(B). On the horizontal axis, the central value of the Sobol’ indices is displayed, which is computed by $\left(\overline{S}_{i}+\underline{S}_{i}\right)/2$. Hence, a variable placed on the right hand side of the plot has a large influence on the uncertainty in the response variable $Y$. On the vertical axis, the epistemic uncertainty of the Sobol’ indices is shown defined as $\overline{S}_{i}-\underline{S}_{i}$. A variable placed at the upper end of the plot has a large epistemic uncertainty attached to the Sobol’ indices. In this example, $F_{1}$ has the largest overall influence, whereas $t_{1}$ has the largest epistemic uncertainty in the estimate. Variables whose representative point lies close to the lower left corner are unimportant ($c_{1}$, $c_{2}$, and $m$ in the current example).

6.3 Simply supported truss

6.3.1 Problem statement

Hurtado Hurtado (2013) introduced the two-dimensional truss structure shown in Figure 7 for the purpose of reliability analysis. The truss is subjected to seven loads $P_{i}$ which are modelled with parametric p-boxes. The loads are defined by lognormal distributions with mean value $\mu_{P_{i}}\in[95,105]~{}{\rm kN}$ and standard deviation $\sigma_{P_{i}}\in[13,17]~{}{\rm kN}$. The geometry of the structure and the material properties are assumed constant. The modulus of elasticity is $E=200\cdot 10^{9}~{}{\rm Pa}$, whereas the cross section area varies along the different bars: $A=0.00535~{}{\rm m}^{2}$ for bars marked by $\bullet$, $A=0.0068~{}{\rm m}^{2}$ for bars marked by $\circ$, and $A=0.004~{}{\rm m}^{2}$ for the remaining bars. The quantity of interest is the deflection at mid-span, denoted by $u_{4}$.

6.3.2 Analysis

The deflection of the truss is computed by a finite element model (FEM), which is implemented in a Matlab-based in-house code. The FEM interprets each bar as a bar element, whereas the loads are modelled as point loads at the intersections of the corresponding bars as indicated in Figure 7.

In order to estimate the imprecise Sobol’ indices, an augmented PCE model is computed based on an experimental design of $N=100$ samples generated by Latin-hypercube sampling in the unit-cube and subsequent isoprobabilistic transform to domain $\mathcal{D}_{\boldsymbol{V}}$. The original computational model requires a seven-dimensional input vector describing the loading forces. Due to the definition of the parametric p-boxes, the augmented PCE model requires a 21-dimensional input vector ($M_{\text{aug}}=21$). The number of phantom points $n_{ph}$ is varied to examine the influence of these points onto the accuracy of the computed Sobol’ indices. The bounds of the Sobol’ indices are obtained by using Matlab’s built-in genetic algorithm with default settings as an optimization tool.

The maximum deflection computed by the FE model corresponds to a monotone function of the load parameters. Hence, the extreme cases of the first order Sobol’ indices can be obtained by setting one variable to $\left\{\mu_{P_{i}}=105~{}\text{kN},\,\sigma_{P_{i}}=17~{}\text{kN}\right\}$ and all others to $\left\{\mu_{P_{i}}=95~{}\text{kN},\,\sigma_{P_{i}}=13~{}\text{kN}\right\}$ or vice versa. The reference solution of the Sobol’ indices are then obtained by PCE and a large experimental design.

6.3.3 Results

The imprecise first-order Sobol’ indices are visualized in Figure 8 as a function of $n_{ph}$. The filled and hollow bars represent $\underline{S}_{i}$ and $\overline{S}_{i}$, respectively. The estimated values (plotted with orange bars) are compared to the reference solution (black bars). As $n_{ph}$ increases, the extreme values of the Sobol’ indices converge to the true values without adding additional FEM runs. In this example, $n_{ph}=5$ is sufficient to obtain reliable results. This impression is confirmed by the relative generalization error which is provided in the caption of the figures. In the case of $n_{ph}=3$ and the corresponding relative generalization error of $\widehat{err}_{gen}=8.4\cdot 10^{-4}$, the imprecise PCE-based analysis provides usable estimates of the imprecise Sobol’ indices. In total, the full set of Sobol’ indices together with their ranges of epistemic uncertainty is obtained for a total cost of 100 finite element analyses of the truss structure, meaning several orders of magnitude less than any Monte Carlo-based double loop approach.

Similarly as for the SDOF oscillator example, Figure 9 shows the impact of aleatory vs. epistemic uncertainty for the truss structure. This plot illustrates nicely the symmetry of the imprecise Sobol’ indices with respect to the central load $P_{4}$, *i.e. *the points representing $P_{3}$ and $P_{5}$, $P_{2}$ and $P_{6}$, $P_{1}$ and $P_{7}$ are identical. Moreover, the amount of imprecision is proportional to the absolute impact for all seven loads in this example.

7 Conclusions

Sensitivity analysis is a popular tool to analyse the input-output relationship of computational models. In this paper, imprecise Sobol’ indices are discussed as a measure of global sensitivity in the context of parametric probability-boxes for modelling the model input parameters. The clear separation of aleatory and epistemic uncertainty in probability-boxes allows for an intuitive interpretation of the resulting interval-valued Sobol’ indices. Therein, the absolute value represents the importance of each parameter with respect to the aleatory uncertainty whereas the range of the interval represents the epistemic uncertainty.

In order to speed up the computations of imprecise Sobol’ indices, an augmented polynomial chaos expansion (PCE) model is introduced. The imprecise Sobol’ indices are obtained by simply post-processing the augmented PCE model. In particular, the extreme values of the Sobol’ indices are computed by optimizing on the coefficients of the augmented PCE model without rerunning the computational model. Hence, the imprecise Sobol’ indices can be evaluated efficiently. Furthermore, the introduction of phantom points counteracts the increased complexity of the augmented PCE model due to the augmented input vector.

The proposed algorithm allows for analysing realistic engineering problems as shown with the single-degree-of-freedom oscillator and the truss structure examples. A small experimental design is sufficient to estimate the imprecise PCE-based Sobol’ indices accurately. Moreover, the “phantom points trick” allows for a further increase in meta-model accuracy and at the same time a more efficient estimation of the Sobol’ indices.