Reliability sensitivity estimation with sequential importance sampling

Structural Safety 75 (2018) 24–34

Contents lists available at ScienceDirect

Structural Safety journal homepage: www.elsevier.com/locate/strusafe

Reliability sensitivity estimation with sequential importance sampling ⁎

T

Iason Papaioannou , Karl Breitung, Daniel Straub Engineering Risk Analysis Group, Technische Universität München, Arcisstr. 21, 80290 München, Germany

A R T I C LE I N FO

A B S T R A C T

Keywords: Reliability analysis Sensitivity analysis Simulation method Sequential importance sampling

In applications of reliability analysis, the sensitivity of the probability of failure to design parameters is often crucial for decision-making. A common sensitivity measure is the partial derivative of the probability of failure with respect to the design parameter. If the design parameter enters the definition of the reliability problem through the limit-state function, i.e. the function defining the failure event, then the partial derivative is given by a surface integral over the limit-state surface. Direct application of standard Monte Carlo methods for estimation of surface integrals is not possible. To circumvent this difficulty, an approximation of the surface integral in terms of a domain integral has been proposed by the authors. In this paper, we propose estimation of the domain integral through application of a method termed sequential importance sampling (SIS). The basic idea of SIS is to gradually translate samples from the distribution of the random variables to samples from an approximately optimal importance sampling density. The transition of the samples is defined through the construction of a sequence of intermediate distributions, which are sampled through application of a resample-move scheme. We demonstrate effectiveness of the proposed method in estimating reliability sensitivities to both distribution and limit-state parameters with numerical examples.

1. Introduction In reliability analysis, the interest is in the assessment of the performance of an engineering system through evaluating its probability of failure. Let X denote a continuous random vector of dimension n modeling the system variables that are expected to present an uncertain behavior. The failure event can be defined as the collection of the outcomes of X for which the so-called limit-state function (LSF) g (x ) takes non-positive values. Consider now a vector θ that collects deterministic parameters that enter the definition of the reliability problem. The vector θ may contain parameters of the joint probability density function (PDF) f (.) of X , denoted by θ d , as well as parameters of the LSF g (.) , denoted by θ g . Hence the determinstic parameter vector can be decomposed as θ = [θ d;θ g ]. The probability of failure can be expressed as follows:

Pf (θ ) =

∫g (x,θ )⩽0 f (x , θd) dx g

(1)

Importantly, the LSF often depends on a computationally intensive numerical model of the engineering system so that evaluation of the integral in Eq. (1) becomes a nontrivial task. To alleviate computational requirements, a variety of tailored approaches have been developed [1,2]. These include approximation methods such as the first/second order reliability method (FORM/SORM) [3], response surface

⁎

Corresponding author. E-mail address: [email protected] (I. Papaioannou).

https://doi.org/10.1016/j.strusafe.2018.05.003 Received 20 December 2017; Received in revised form 12 April 2018; Accepted 12 May 2018 0167-4730/ © 2018 Elsevier Ltd. All rights reserved.

approaches [4] and simulation techniques based on the Monte Carlo method. All these methods have their merits and disadvantages as discussed in the literature (e.g. [5–8]). In this contribution we focus on sampling-based methods. The Monte Carlo method is a simple and robust technique, that is able to handle any LSF, independent of its complexity. The efficiency of the standard Monte Carlo method does not depend on the dimension of the random variable space. The probability integral of Eq. (1) can be expressed as the expectation of the indicator function I (g (x , θ g ) ⩽ 0) , where I (g (x , θ g ) ⩽ 0) = 1 if g (x , θ g ) ⩽ 0 and I (g (x , θ g ) ⩽ 0) = 0 otherwise. Standard Monte Carlo estimates Pf (θ ) for a given parameter vector θ by generating ns independent samples {x (k ) , k = 1, …, ns} from the PDF f (x , θ d ) and taking the sample mean of I (g (x (k ) , θ g ) ⩽ 0) , i.e.

1 Pf = Ef̂ [I (g (x , θ g ) ⩽ 0)] = ns

ns

∑

I (g (x (k ) , θ g ) ⩽ 0)

k=1

(2)

The estimate of Eq. (2) is unbiased and has coefficient of variation:

δ Pf =

1−Pf (θ ) ns Pf (θ )

(3)

δ Pf is a measure of the statistical accuracy of Pf . According to Eq. (3), crude Monte Carlo requires approximately 10 k + 2 samples to achieve an accuracy of δ Pf = 10% for a probability in the order of 10−k . Hence, the

Structural Safety 75 (2018) 24–34

I. Papaioannou et al.

state function. Such methods are the directional sampling [28,29] and LS [13,14] methods. Both approaches require the solution of a line search problem for each sample to determine the intersection of the sampling direction or line with the limit-state surface. Application of these methods to reliability sensitivity is discussed in [29–32]. Alternative methods for solving the problem of Eq. (5) have been proposed in [33] and in [34]. The procedure introduced in [33] requires solving the limit-state equation g (x , θ g ) = 0 for one component of x . Estimation of the reliability sensitivity is achieved by conditional IS in terms of the remaining random variables. The method presented in [34] involves a linear approximation of the LSF in terms of the parameter vector and an approximation of the probability of failure in terms of the LSF. The approximations are constructed perturbing the limitstate parameters at the samples close to the failure surface and performing additional LSF evaluations. In [31], the authors introduced an approximation of the surface integral of Eq. (5) with a domain integral. This approach allows application of any simulation method for estimating Eq. (5). The same approximation was independently proposed in [35]. Therein, the domain integral approximation is estimated with SuS. In [36] this approach was combined with a polynomial chaos surrogate for efficient approximate reliability sensitivity analysis. In this contribution, we propose estimation of the domain integral approximation with SIS [18]. The basic idea of SIS is to gradually translate samples from the distribution of the random variables to samples from an approximately optimal importance sampling density. The transition of the samples is defined through the construction of a sequence of intermediate distributions. We demonstrate that this sequence is particularly suitable for estimating the domain integral approximation of Eq. (5). The proposed approach can also be applied for estimating Eq. (4) if estimation is performed at the independent standard normal space. The structure of the paper is as follows. In Section 2, we review the approximation of reliability sensitivity to limit-state parameters through a domain integral. Section 3 discusses estimation of the approximation with SIS. Section 4 focuses on application of the proposed approach to evaluation of reliability sensitivity of system problems with multiple failure modes. Section 5 tests the proposed method with a series of numerical examples. The paper closes with the conclusions in Section 6.

computational cost of crude Monte Carlo becomes intractable for low target failure probabilities. Several methods have been proposed that aim at enhancing the efficiency of the crude Monte Carlo method through decreasing the variance of the probability estimate. These include importance sampling (IS) and its adaptive variants [9–12], line sampling (LS) [13,14], subset simulation (SuS) [15,16] and sequential importance sampling (SIS) [17,18]. In many practical applications of reliability analysis, one is interested in understanding the influence of each component of the vectors X and θ on the probability of failure. The influence of the random variables X can be quantified by a variety of sensitivity analysis approaches, e.g. the FORM sensitivity indices [3] or other variance-based sensitivity measures (e.g. [19–21]). The sensitivity to deterministic parameters θ can be quantified through performing parameter studies for a range of parameter values or through evaluating the local partial derivative derivatives of Pf (θ ) to the components of θ . In this paper, we focus on efficient estimation of the latter. Approximations of the partial derivative of Pf (θ ) in terms of both distribution parameters θ d and limit-state parameters θ g are obtained as a byproduct of FORM/SORM solutions [1,22]. They involve the evaluation of first- or second-order derivatives of the isoprobabilistic transformation to the independent standard normal space and the derivative of the limit-state function g (.) to the respective parameters at the most probable failure point. The derivative of the probability of failure Pf (θ ) with respect to a distribution parameter θid is obtained through differentiating the integrand of Eq. (1) as follows:

∂Pf (θ ) ∂θid

=

∫g (x,θ )⩽0 g

∂f (x , θ d ) dx ∂θid

(4)

The expression in Eq. (4) is a (possibly high dimensional) domain integral and can be estimated using any reliability method including standard Monte Carlo methods. This has been acknowledged by Wu in [23], who estimated distribution parameter sensitivities with standard IS. Other variance reduction methods can also be applied, such as SuS [24,25] and adaptive IS based on surrogate models [26]. However, it should be noted that these approaches cannot be implemented in a transformed random variable space, which is commonly employed for estimating Pf (θ ) . This can be understood by considering that when transforming the random variable space to an equivalent standardized space (e.g. the independent standard normal space), distribution parameters become parameters of the transformed LSF, while the standardized distribution is parameter-free. As discussed next, reliability sensitivities to limit-state parameters are not expressed as domain integrals. We note that several sampling-based reliability methods benefit from such a transformation; e.g. it facilitates the choice of appropriate IS functions [2] and it can provide a basis for dimensionalindependent performance of sequential sampling methods [16,27,18]. A different approach was introduced in [19] for the estimation of distribution parameter sensitivities. The authors constructed a linear response surface using Monte Carlo samples and applied the FORM sensitivity results to this surrogate model. It is noted that this approach can also be applied to evaluate sensitivities to limit-state parameters, since the latter are also provided as byproducts of FORM. Consider the case where the LSF g (x , θ g ) is continuously differentiable and ∇x g (x , θ g ) ≠ 0 for all x and θ g on the surface {g (x , θ g ) = 0} . The derivative of Pf (θ ) with respect to a parameter of the LSF θig is then given by the following surface integral [22]:

∂Pf (θ ) ∂θig

=−

∫g (x,θ )=0 f (x , θd) ‖∇x g (x1 , θ g )‖ g

∂g (x , θ g ) ds (x ) ∂θig

2. Approximate reliability sensitivity analysis In this section, we derive an approximation of the surface integral that arises in reliability sensitivities to limit-state parameters through a domain integral. This approximation, originally introduced in [31], enables estimation of the surface integral with standard Monte Carlo methods. We focus on sensitivities to limit-state parameters, i.e. the case where θ = θ g . We note that if the reliability problem is recast in a standardized probability space through performing an isoprobabilistic transformation (e.g. [37]), distribution parameters become limit-state parameters. Therefore, in such settings the methods described here can be applied to the computation of reliability sensitivities to any deterministic parameter. 2.1. Approximation of the surface integral We derive the approximation of the surface integral of Eq (5) through considering the following definition of the probability of failure:

Pf (θ ) =

(5)

θ g)

= 0} . where ds (x ) denotes surface integration over the surface {g (x , Eq. (5) is a surface integral and cannot be estimated with classical Monte Carlo methods. However, one can estimate Eq. (5) through application of simulation methods that monitor the boundary of the limit-

∫

n

I (g (x , θ ) ⩽ 0) f (x ) dx

(6)

where I (g (x , θ ) ⩽ 0) is the indicator function of the failure domain and g (x , θ ) is a LSF that is continuously differentiable almost everywhere with respect to any component θi of θ . In the above expression, the dependence of parameters θ enters in the integrand. Therefore, one 25

Structural Safety 75 (2018) 24–34

I. Papaioannou et al.

proposed in [35] through considering a smooth approximation of the Dirac function, resulting in an approximation of the integration bound in Eq. (5). Therein, a number of smooth approximations of the Dirac function were considered and it was numerically shown that Eq. (11) results in the smaller mean-square error for an appropriate choice of the smoothing parameter. In the Appendix we show that under some regularity conditions Eq. σ → 0. A (10) indeed converges to the true surface integral of Eq (5) as ∼ further analysis for deriving error bounds for Eq. (10) has not yet been carried out and is left to future work. In the next section, we discuss σ in relation to the uncertainty of estimates appropriate choices of ∼ obtained by simulation methods.

could obtain an expression for the derivative of Pf to any component of θ through integrating the derivative of the integrand of Eq. (6). However, this is not possible as the indicator function is not differentiable. The key idea for deriving a domain integral approximation of Eq (5) is to substitute the indicator function in Eq. (6) with a smooth approximation and then take the derivative of the resulting approximation of the probability of failure. The indicator function I (g (x , θ ) ⩽ 0) can be expressed by the following limit

I (g (x , θ ) ⩽ 0) = lim Φ ⎛− σ→0 ⎝

g (x , θ ) ⎞ σ ⎠

(7)

where Φ(.) is the standard normal cumulative distribution function (CDF). It is noted that this is only one of several equivalent definitions of the indicator function through limits of smooth functions (e.g. see σ , with ∼ σ small enough, we can approximate [38]). Choosing σ = ∼ I (g (x , θ ) ⩽ 0) by the following expression

I (g (x , θ ) ⩽ 0) ≈ Φ ⎛− ⎝

g (x , θ ) ⎞ ∼ σ ⎠

σ 2.2. Choice of parameter ∼ The quality of the approximation of Eq. (10) depends on the choice σ : the bias of the approximation decreases with deof the parameter ∼ σ . However, a very small ∼ σ would lead to a large creasing value of ∼ uncertainty in the corresponding estimate of Eq. (10) obtained by any σ the domain Monte Carlo method. This is due to the fact that for small ∼ integral of Eq. (10) converges to a surface integral over the limit-state surface g (x , θ ) = 0 . Hence most of the samples values in Eq. (11) away σ should from this surface will be zero. Therefore, an optimal choice of ∼ achieve a compromise between the bias and the coefficient of variation of the corresponding Monte Carlo estimate. Here, we illustrate this with the help of a simple numerical example. Consider the limit-state function g (u, β ) = −u + β , where u is the outcome of a standard normal random variable and β is a deterministic parameter. It is noted that any linear function of normal random variables can be reduced to a similar form. The probability of failure equals Φ(−β ) and the sensitivity with respect to β is −φ (β ) . Let us assume that the probability of failure is estimated with crude Monte Carlo, with a target coefficient of variation of 5%. Following Eq. (3), the required number of samples is ns = 0.05−2 /Φ(−β ) = 400/Φ(−β ) . We use the same samples for estimating the reliability sensitivity. We study the behavior of the resulting Monte Carlo estimate of the reliability sensitivity for σ. different values of the parameter ∼ Fig. 2 plots the relative bias, i.e. the absolute difference between the approximation and the true value of the sensitivity divided by the absolute value of the latter, the coefficient of variation of the sensitivity estimate and the Euclidean norm of the two, for β = 3.5 and and difσ . The norm of the two measures is ferent values of the parameter ∼ equivalent to the mean-square error of the estimate of the sensitivity divided by the absolute value of the true sensitivity. The plots show that the minimum value of the norm, which corresponds to the optimal σ , would be acceptable for most problems, which supports the choice of ∼ applicability of the approximation of Eq. (10) in Monte Carlo settings. σ leads to slightly larger coefficient of variation The optimal value of ∼ than the one of the probability estimate. σ cannot be computed, In practical situations, the optimal value of ∼ since the true value of the probability sensitivity is not known and hence the relative bias of the approximation cannot be estimated. One σ and choose the minimum ∼ σ possible remedy is to vary the value of ∼ that leads to a coefficient of variation smaller than a target value, e.g. the same target value used for the probability of failure. It is noted that the same samples of the limit-state function and its derivative with respect to the parameter can be used for estimation of the coefficient of σ . Therefore, the additional computavariation for different values of ∼ tional cost for this adaptive approach is negligible. In [35] an empirical σ is proposed that depends on the number of failed samples in choice of ∼ a Monte Carlo analysis. However, this choice depends only on the magnitude of the probability of failure and does not account for possible variability of the derivatives of the limit-state function in the vicinity of the failure surface. Next, we employ an adaptive importance sampling approach to estimate the approximation of Eq. (10), which σ. uses an adaptive approach for choosing the parameter ∼

(8)

The approximation of the indicator function in Eq. (8) is illustrated σ . The concept of a smooth approxin Fig. 1 for different values of ∼ imation of the indicator function has been used in numerous fields of engineering and science. In reliability analysis, it has been used for deriving adaptive importance sampling approaches e.g. in [17,26,18], see also Section 3.1. Inserting Eq. (8) into Eq. (6), we obtain an ap∼ σ ) and exproximation of the probability of failure denoted by Pf (θ, ∼ pressed as:

∼ Pf (θ, ∼ σ) =

∫

n

Φ ⎛− ⎝

g (x , θ ) ⎞ f (x ) dx ∼ σ ⎠

(9)

Taking the derivative of the above expression with respect to a limit-state parameter θi , we arrive at the following approximation of the reliability sensitivity: ∼ ∂Pf (θ, ∼ σ) g (x , θ ) ∂g (x , θ ) 1 f (x ) dx = − n ∼φ⎛ ∼ ⎞  σ ∂θi (10) ⎝ σ ⎠ ∂θi

∫

wherein φ (.) is the standard normal PDF. The expression in Eq. (10) is a domain integral and can be estimated by any standard simulation method, e.g. crude Monte Carlo or importance sampling. For example, the crude Monte Carlo estimate of Eq. (10) reads:

1 ∼̂ Pf , θi = ns

ns

∑ k=1

(k ) (k ) ⎡− 1 φ ⎛⎜ g (x , θ ) ⎞⎟ ∂g (x , θ ) ⎤ ∼ ⎥ ⎢ ∼ σ ∂θi ⎠ ⎦ ⎣ σ ⎝

(11)

{x (k ) ,

where k = 1, …, ns} are independent samples from f (x ) . We note that the approximation of Eq. (10) was independently

Fig. 1. Illustration of the approximation of the indicator function for different σ. values of ∼ 26

Structural Safety 75 (2018) 24–34

I. Papaioannou et al.

Fig. 2. Relative bias and coefficient of variation of estimate of parameter sensitivity; Linear limit-state in terms of normal random variable and in terms of limit-state parameter.

3. Sequential importance sampling for reliability sensitivity

where ∞ = σ0 > … > σL > 0 . SIS samples the distributions {hj (x ), j = 0, …, L} in a step-wise manner and estimates each normalizing constant Pj by IS using as IS density the function hj − 1. Assume that samples {x (k ) , k = 1, …, ns} distributed according to hj − 1 are available. It is

We introduce a methodology for reliability sensitivity analysis based on an adaptive sampling method termed sequential importance sampling (SIS). We first summarize the basic principle of SIS for reliability analysis – for a more detailed description the reader is referred to [18]. Thereafter, we discuss the post processing step for estimation of reliability sensitivities.

Pj =

∫

n

ηj (x ) dx = Pj − 1

where wj (x ) =

ηj (x ) ηj − 1 (x )

=

3.1. Sequential importance sampling for reliability analysis

malizing constants Sj =

IS aims at decreasing the variance of the crude Monte Carlo estimate of the probability of failure through sampling from an alternative sampling density, the so-called IS density h (x ) . For notational convenience we drop the dependence of the probability of failure on the parameter vector θ in this section. The integral in Eq. (6) can be modified as follows:

Sj ̂ =

Pf =

∫

n

I (g (x ) ⩽ 0) w (x ) h (x ) dx = Eh [I (g (x ) ⩽ 0) w (x )]

(12)

where w (x ) = is the importance weight function. The IS estimate of Pf is given as follows: ns

∑ k=1

I (g (x (k ) ) ⩽ 0) w (x (k ) ) (13)

{x (k ) ,

where k = 1, …, ns} are independent samples from h (x ) . The probability estimate of Eq. (13) is unbiased provided that the support of h (x ) contains the failure domain {g (x ) ⩽ 0} . An appropriate choice of the IS function can lead to significantly smaller variance compared to the one of the crude Monte Carlo estimate. The theoretically optimal IS density is given by the following expression:

hopt (x ) =

1 I (g (x ) ⩽ 0) f (x ) Pf

1 ⎛ g (x ) ⎞ 1 Φ ⎜− η (x ) ⎟ f (x ) = Pj ⎝ σj ⎠ Pj j

wj (x ) hj − 1 (x ) dx = Pj − 1 Ehj − 1 [wj (x )]

Φ(−g (x ) / σj ) Φ(−g (x ) / σj − 1 ) Pj Pj − 1

(16)

. An estimate of the ratio of nor-

is given by:

Pj 1 = Eĥ j − 1 [wj (x )] = ns Pj − 1

P Pf = M ns

(14)

ns

∑

wj (x (k ) )

(17)

k=1

ns

∑ k=1

I (g (x (k ) ) ⩽ 0)

f (x (k ) ) P = M ns ηM (x (k ) )

ns

∑ k=1

I (g (x (k ) ) ⩽ 0) Φ(−g (x )/ σL )

(18)

wherein the estimate PL of the normalizing constant of hL is obtained as

The IS function of Eq. (14) leads to a variance of the IS estimate of zero. However, the optimal IS function cannot be used in practice since it requires knowledge of the failure domain and target probability Pf . The basic idea of SIS is to adaptively sample a sequence of distributions that gradually approximate the optimal IS density h opt (x ) . An approximation of h opt (x ) is obtained through inserting the smooth approximation of the indicator function introduced in Eq. (8) into Eq. (14). The resulting approximation of h opt (x ) becomes more accurate σ . We define the distribution with decreasing value of the parameter ∼ sequence {hj (x ), j = 0, …, L} , with

hj (x ) =

n

Sampling from the sequence of Eq. (15) can be performed through a resample-move scheme starting from samples from h 0 (x ) = f (x ) . In each step, samples from distribution hj − 1 are weighted according to the importance weight function wj (x ) . The resulting weighted samples are resampled to obtain uniformly weighted samples from hj through application of a resampling method (e.g. [39]). In addition, to decrease sample degeneracy the samples are moved to regions of high probability density of hj through Markov chain Monte Carlo (MCMC) sampling. For performing the latter step, tailored MCMC algorithms are described in [18]. These include an independent Metropolis–Hastings algorithm for application to general problems in low- to moderate-dimensional random variables spaces and a conditional sampling approach that is especially efficient for tackling high-dimensional problems. These algorithms are optimized for application in standard normal space. In the final sampling step, samples {x (k ) , k = 1, …, ns} from the distribution hL are used to estimate Pf with IS, as follows:

f (x ) h (x )

1 Pf = Eĥ [I (g (x ) ⩽ 0) w (x )] = ns

∫

PL =

L

∏

Sj ̂

j=1

(19)

3.2. Adaptive choice of parameter σj A crucial ingredient of the SIS approach is the choice of the parameters {σj, j = 1, …, L} of the intermediate distributions. To ensure effectiveness of the IS approach for estimating each ratio of normalizing constants, each parameter σj should be chosen such that distribution hj is not too different from hj − 1. However, σj should not be too close to σj − 1 because this would increase the total number of intermediate sampling

(15) 27

Structural Safety 75 (2018) 24–34

I. Papaioannou et al.

σ as small as possible such that the coefficient of failure Pf , we choose ∼ ∼ variation of Pf ,̂ θi is close to δ Pf . Assuming again that the samples at each sampling level are independent from samples at other levels, the coef∼ ficient of variation of P ̂ is given as:

steps and hence the computational cost of SIS. In [17] the parameters σj are chosen using a deterministic mathematical sequence. Here, we employ the approach proposed in [18] where σj is chosen adaptively, such that the sample coefficient of variation δ ŵ of the importance j

f , θi

weights {wj (x (k ) ), k = 1, …, ns} is equal to a target value δtarget . In this approach, at each sampling step of the SIS procedure, one solves the following optimization problem

σj = argmin (δ ŵ j (σ )−δtarget )2 σ ∈ (0, σj − 1 )

L

δ∼ P̂

f , θi

I (g (x (k ) ) ⩽ 0) Φ(−g (x (k ) ) / σL )

that δ ∼ P̂

f , θi

δS2 ̂

, k = 1, …, ns } is smaller than δtarget and L is set to the

f , θi

current step j. In [40] it is noted that this adaptive approach of choosing the parameters of the intermediate distributions based on the coefficient of variation of the weights is equivalent to requiring that the so-called effective sample size (ESS) takes a target value. The ESS, which measures the number of samples related to an independent sampling approach, is often used in adaptive sequential Monte Carlo approaches for Bayesian inference (e.g. [41]). We discuss now how this adaptive approach allows one to control the variability of the SIS estimate of the probability of failure. For simplicity, we make the assumption that at each sampling step j of the SIS approach the available samples from hj − 1 are independent. In addition, we assume that these samples are independent from samples at other sampling levels. We note that both assumptions could be satisfied if a sufficient burn-in period is applied in the MCMC step. The coefficient of variation of the SIS estimate in Eq. (18) is then obtained as follows:

∑

+

δS2 ̂ f

1 2 δw ns j

(26)

Sf , θi

Sj

where δtarget is the same target value used in the program of Eq. (20). We note that solving Eqs. (20) and (27) does not require additional LSF evaluations; therefore the additional computational cost is negligible. The upper bound in the allowable space of the program of Eq. (27) is set equal to σL , which is the parameter of the sampling density hL . This is done to ensure that the support of the integrand in the approximation of the sensitivity is included in the finite sample approx∼ imation of the density h used to estimate P ̂ . Using this upper bound

(21)

4. Reliability sensitivity of system problems We now discuss estimation of the sensitivity of the probability of failure for reliability problems with multiple modes of failure. Consider a system problem consisting of m component failure modes defined by corresponding continuously differentiable LSFs gi (x , θ ), i = 1, …, m . We define the generalized LSF gsys (x , θ ) as follows

(23)

That is, requiring that an estimate of δ wj takes a target value as done in Eq (20) allows us to bound the coefficient of variation of the SIS probability estimate. In [18] it is suggested to choose δtarget = 1.5, which corresponds to ε ≈ 0. 05 L + 1 for ns = 1000 .

gsys (x , θ ) = min [max{ai gi (x , θ )}] 1 ⩽ k ⩽ K i ∈ Ck

Having estimated the probability of failure with SIS, one can use the available samples from density hL to estimate the approximation of the reliability sensitivity of Eq. (10) with IS. The corresponding estimate reads: ns k=1

(k ) (k ) 1 ⎤ ⎡− 1 φ ⎛⎜ g (x , θ ) ⎞⎟ ∂g (x , θ ) ∼ (k ) )/ σ , θ ) ⎥ ⎢ ∼ σ θ σ ∂ Φ( g ( x − i L ⎝ ⎠ ⎦ ⎣

(28)

wherein Ck ⊆ {1, …, m} are index sets of component failure modes, whose joint failure results in system failure. This corresponds to the cutset representation of a system. The failure domain of the general system reliability problem can be expressed {gsys (x , θ ) ⩽ 0} . From Eq. (28) we can retrieve the component problem by setting m = 1, the parallel system problem by setting K = 1 and the series system problem by defining each set Ck by a single index. Similar definitions of the general system problem can be found in [3,15]. In the definition given here each component LSF is multiplied by a constant ai . This is done to avoid the situation where certain component LSFs dominate the behavior of gsys away from the surface gsys (x , θ ) = 0 . This is often the case when each component LSF is expressed in different units. A possible choice of the constants is ai = E[gi (x , θ )]−1 for some value of θ . Using the LSF of Eq. (28), the probability of failure of the general system problem can be estimated with SIS following the procedure discussed in Section 3.1. Each constant ai can be evaluated through

3.3. Estimation of the reliability sensitivity

∑

f , θi

might result in δ ŵ f , θi > δtarget , especially in cases where δtarget is set to a low value. In such cases, a small number of samples will have nonzero values in the sample approximation of the sensitivity. However, alσ > σL , which could lead to a δ ŵ f , θi ≈ δtarget , would also inlowing for ∼ crease considerably the bias in the sample approximation of the sensitivity because the samples from hL would be concentrated in a small area of the support of the integrand in Eq. (10).

(22)

ns ε L+1

(27)

σ ∈ (0, σL )

where δ wj is the coefficient of variation of wj (x ) . Bounding Eq. (21) by a small number ε 2 is equivalent to bounding the coefficient of variation of 2 the intermediate estimates by ε . This implies that for every L+1 {j = 1, …, L, f }

P ∼̂ Pf , θi = L ns

1 2 δw ns f , θi

∼ σ = argmin (δ ŵ f , θi (σ )−δtarget )2

wherein Sf ̂ = From Eq. (17) and because of the independence assumption of the samples from each hj − 1, we have:

δ wj ⩽

Comparing Eq. (25) with Eq. (21), one observes ≈ δ Pf if δSf̂ , θi is close to δSj,̂ j = 1, …, L . It is

solving

Pf . PL

j

(25)

∼̂ Pf , θi . PL

L

δS2 ̂ j

j=1

δS2 ̂ =

=

f , θi

j

where δ wf , θi is the coefficient of variation of the weight function wf , θi (x ) expressing the likelihood ratio of the integrand of Eq. (10) to the IS density hL (x ) . An approximation of δ wf , θi is obtained by taking the sample coefficient of variation δ ŵ f , θi of the term in brackets in Eq. (24). σ through To guarantee that δ ̂ ≈ δ ̂ for j = 1, …, L , we evaluate ∼

L

δ Pf =

δS2 ̂ + δS2 ̂

j=1

wherein Sf ̂ , θi =

(20)

The procedure is stopped when the coefficient of variation of the weights with respect to the optimal IS density

{

∑

=

(24)

As discussed in Section 2.2, an appropriate choice of the parameter ∼ σ in Eq. (24) should achieve a compromise between the bias and coefficient of variation of the estimate of the reliability sensitivity. The σ decreases. To enbias in the approximation of Eq. (10) decreases as ∼ sure consistency with the accuracy of the estimate of the probability of 28

Structural Safety 75 (2018) 24–34

I. Papaioannou et al.

Y

estimating E[gi (x , θ )]−1 using the samples from f (x ) at the first sampling level of SIS. We now look at estimation of the sensitivity of the probability of system failure to a component θj of the parameter vector θ . We consider again the approximation of Eq. (10). Application of this approximation to system problems would require an expression for the first derivative of gsys with respect to θj . The function gsys is not differentiable at all points [x , θ]: it is not defined for all [x , θ] where gi (x , θ ) = gj (x , θ ), i, j = 1, …, m , i ≠ j ; otherwise, it takes the following form:

∂gsys (x , θ ) ∂θj

∂garg =

min [max {ai gi (x , θ )}] (x ,

w L=100" Fig. 3. Cantilever beam subjected to biaxial bending according to [42].

Table 1 Uncertain parameters of Example 1.

θ)

1 ⩽ k ⩽ K i ∈ Ck

∂θj

Distribution

Mean

CV

X [lb] Y [lb] y [psi] E [psi]

Normal Normal Normal Normal

500 1000 40000

20% 10% 5% 5%

29 × 106

600 600 g11 (Y , X , y, w, t ) = y−⎛ 2 Y + 2 X ⎞ wt ⎠ ⎝ wt The second LSF restricts the maximum allowed displacements at the tip of the beam:

g12 (Y , X , E , w, t , D0) = D0−

5. Numerical examples

4L3 Y 2 X 2 ⎛ ⎞ +⎛ ⎞ 2 Ewt ⎝ t ⎠ ⎝ w2 ⎠

The two loads X and Y, the Young’s modulus E and the yield strength y of the beam are modeled by independent normal random variables, as shown in Table 1. We note that the normal distribution is strictly not an appropriate choice for modeling parameters with positive support, such as E and y; however, we use the same probabilistic model as in [42] for the sake of comparison. The parameters of the limit-state functions are the width w and height t of the cross-section of the beam and the allowable tip displacement D0 . The parameter sensitivities of the probability of failure are estimated at (w, t , D0) = (2.4, 3.9, 2.5) [in], which corresponds to the optimal reliability-based design for an allowable probability of failure of 3 × 10−3 . We compare the results obtained by SIS with the ones obtained by the LS method discussed in [31]. In the LS, the sampling direction is determined by an initial Monte Carlo run with 100 samples, while the samples size for the LS procedure is set to 500. In SIS, we use 1000 samples per level and set the target coefficient of variation of the weights for selection of the intermediate distributions to δtarget = 1.5. Table 2 shows the results for the LSF g11. The function g11 is linear in the normal random variables and therefore we evaluate the exact results for both the probability of failure and the sensitivities by running a FORM analysis. The mean estimates of the approximation of the sensitivities obtained by SIS are very close to the exact results, while the coefficient of variation (CV) of the estimates is consistent with the one of the probability estimate. The number of derivative evaluations in SIS indicates the number of unique samples among the samples in the final level of SIS obtained by the MCMC algorithm. The mean estimates obtained by LS show excellent agreement with the reference values. We note that LS requires the evaluation of the gradient of the LSF in terms of the random variables at each sample, in addition to the gradient in terms of the deterministic parameters, while SIS requires only the latter. The CVs of the LS estimates are moderately lower than the ones of the SIS estimates. In Table 3 the results for the LSF g12 are shown. This function is nonlinear in all three random variables. The reference solution in this case is evaluated by a LS simulation with 10 4 samples, preceded by a FORM analysis for determining the important direction pointing to the design point. The mean estimates of the sensitivities obtained by SIS are again very close to the exact results and the CVs are close, albeit slightly higher, to the CV of the probability estimate. The LS results are also

In this section, we investigate the performance of the proposed SIS method for reliability sensitivity analysis with numerical examples. The MCMC step required for sampling the intermediate distributions in SIS is performed with the independent Metropolis–Hastings algorithm proposed in [18], which uses as proposal density the Gaussian mixture model fitted with the weigthed samples. SIS is implemented in the independent standard normal space where all deterministic parameters are parameters of the transformed LSF. The first example examines the performance of the method for estimating sensitivities to limit state parameters while in the second and third example the method is applied to the estimation of sensitivities to distribution parameters. The first two examples consist of component reliability problems and the third example is a series system problem. For the second and third example, the reliability sensitivity results are compared in terms of the so-called partial elasticities (e.g. [22]). The partial elasticity e θi is defined as

∂Pf θi ∂θi Pf

Parameter

(29)

We note that evaluation of Eq. (10) requires the derivative of gsys to be continuous almost everywhere. This requirement holds under the assumption that all gi (x , θ ), i = 1, …, m , have continuous first derivatives. When applying SIS to estimate Eq. (10), the resulting estimate is a sample mean based on weighted samples from the continuous random vector X . In continuous probability spaces, the probability of occurrence of a sample point for which gi (x , θ ) = gj (x , θ ) is zero. Therefore, it is straightforward to apply the procedure described in Section 3.3 to estimate the reliability sensitivity of system failure with SIS; we only need to substitute the expression of Eq. (29) for the derivative of the LSF in Eq. (24).

e θi =

X

t

(30)

and expresses the approximate change in percentage of the probability of failure if θi is changed by 1%. The statistics of the probability and sensitivity estimates are computed with 500 independent simulation runs. For the first example, we compare the performance of the proposed SIS method with the LS method [31]. We note that, as discussed in Section 1, LS is one of a few sampling approaches that estimate the exact surface integral of Eq. (5). For the second and third example, we compare the results of SIS with the ones of the standard IS approach of Wu [23] using an IS density informed by a FORM analysis, which is common practice in reliability analysis [2]. 5.1. Cantilever beam The first example is a cantilever beam subjected to biaxial bending (Fig. 3). This example was studied in [42] in the context of reliabilitybased optimization. Two LSFs are considered. The first one represents yielding at the fixed end of the beam: 29

Structural Safety 75 (2018) 24–34

I. Papaioannou et al.

Table 2 Results for limit state function g11 of Example 1. The target coefficient of variation of the weights in SIS is chosen as δtarget = 1.5. Reference

LS (100 + 500 samples) Mean estimate

St. dev. (CV)

SIS (1000 samples/level) Mean estimate

St. dev. (CV)

Pf

3.03 × 10−3

3.03 × 10−3

2.18 × 10−4 (7.2%)

3.01 × 10−3

3.07 × 10−4 (10.2%)

∂Pf / ∂w

−5.76 × 10−2

−5.76 × 10−2

3.45 × 10−3 (6.0%)

−5.81 × 10−2

6.62 × 10−3 (11.4%)

∂Pf / ∂t

−3.53 × 10−2

−3.53 × 10−2 1600 500 500

2.12 × 10−3 (6.0%)

−3.57 × 10−2 4000 788

4.07 × 10−3 (11.4%)

g-eval ∇x g -eval ∇θ g -eval

very close to the exact solutions, however, requiring more LSF evaluations than for g11 because of the nonlinearity of the line search problem that needs to be solved for each sample. For this LSF, the two methods resulted in comparable CVs. 5.2. Roof truss The second example consists of the reliability assessment of a roof truss under a serviceability limit state (Fig. 4). This example was studied in [24] in the context of reliability sensitivity analysis. The LSF is defined as follows:

g2 (q, l, Ac , Ec , As , Es ) = 0.03−

Fig. 4. Roof truss according to [24].

ql 2 ⎛ 3.81 1.13 ⎞ + 2 ⎝ Ac Ec As Es ⎠ ⎜

Table 4 Uncertain parameters of Example 2.

⎟

where q is the distributed load applied on the roof truss, l is the roof span, Ac , Ec and As , Es are the cross sections and Young’s moduli of the concrete and steel beams, respectively. All these parameters are modeled by independent normal random variables, as shown in Table 4. Again, the normal distribution is strictly not an appropriate choice for modeling parameters l, Ac , Ec , As , Es due to their positive support; however, we use the same probabilistic model as in [24] for the sake of comparison. Table 5 shows the statistics of the estimated probability of failure and partial elasticities of the probability in terms of the means and standard deviations of the random variables. We compare the results obtained by SIS with the ones obtained by IS in the original space following the approach of [23]. As IS function, we employ an independent normal distribution centered at the design point obtained by a FORM analysis and with standard deviations equal to the standard deviations of the random variables. The reference solution is obtained by the aforementioned approach using 5 × 105 samples. In SIS, we set again the target coefficient of variation of the weights to δtarget = 1.5. The results show that the IS approach with 1000 samples is able to provide an accurate estimate of the probability of failure with a CV of about 5%. We note that in the number of LSF evaluations for the IS approach we do not count the LSF evaluations for solving the FORM optimization problem. The mean estimates of elasticities obtained by IS are accurate, however, the CV of the estimates varies considerably; it ranges from 3.2% to 45.5%. On the other hand, SIS with 1000 samples

Parameter

Distribution

Mean

CV

q [N/m] l [m]

Normal Normal Normal

20000 12

7% 1% 6%

As [m2 ]

Ac [m2 ] Es [N/m2 ] Ec [N/m2 ]

Normal

9.82 × 10−4 0.04

Normal

1 × 1011

6%

Normal

2 × 1010

6%

12%

per level results in slightly biased estimates, as the method evaluates an approximation of the sensitivities; however, the CV of the estimates has significantly lower range (from 6.8% to 13.6%). This is a result of the adaptive approach described in Section 3.3, which ensures that the σ of the approximation is chosen such that the CV of the parameter ∼ sensitivities is close to the CV of the probability of failure. It is noted that the IS approach in the original random variables space requires the evaluation of derivatives of the PDF of the random variables at each sample, which does not involve additional LSF evaluations. If the sampling is performed in standard normal space, the distribution parameters become parameters of the LSF. Hence for estimating the approximation of the sensitivities with SIS, we need to evaluate the gradient of the LSF for each sample at the final sampling level.

Table 3 Results for limit state function g12 of Example 1. The target coefficient of variation of the weights in SIS is chosen as δtarget = 1.5. Reference

LS (100 + 500 samples) Mean estimate

St. dev. (CV)

SIS (1000 samples/level) Mean estimate

St. dev. (CV)

Pf

2.52 × 10−4

2.52 × 10−4

3.55 × 10−5 (14.1%)

2.45 × 10−4

3.14 × 10−5 (12.8%)

∂Pf / ∂w

−8.84 × 10−3

−8.85 × 10−3

1.05 × 10−3 (11.9%)

−8.84 × 10−3

1.31 × 10−3 (14.8%)

∂Pf / ∂t

−2.95 × 10−3

−2.94 × 10−3

3.59 × 10−4 (12.2%)

−2.94 × 10−3

4.35 × 10−4 (14.8%)

∂Pf / ∂D0

−3.27 × 10−3

−3.27 × 10−3 2640 500 500

3.92 × 10−4 (12.0%)

−3.28 × 10−3 5000 724

4.85 × 10−4 (14.8%)

g-eval ∇x g -eval ∇θ g -eval

30

Structural Safety 75 (2018) 24–34

I. Papaioannou et al.

Table 5 Results for limit state function of Example 2. The target coefficient of variation of the weights in SIS is chosen as δtarget = 1.5. Note that the IS solution requires prior knowledge of the FORM design point. Reference

IS (1000 samples)

SIS (1000 samples/level)

Mean estimate

St. dev. (CV)

Mean estimate

St. dev. (CV)

9.38 × 10−3 23.60

9.39 × 10−3 23.60

5.26 × 10−4 (5.6%) 0.76 (3.2%)

9.24 × 10−3 24.68

8.87 × 10−4 (9.6%) 1.68 (6.8%)

e μl

51.56

51.12

6.08 (11.9%)

53.04

3.61 (6.8%)

e μ As

−19.49

−19.47

1.05 (5.4%)

−20.36

1.38 (6.8%)

Pf

e μq

e μ Ac

−9.11

−9.09

0.65 (7.2%)

−9.48

0.68 (7.2%)

e μ Es

−19.46

−19.54

1.02 (5.2%)

−20.36

1.38 (6.8%)

e μ Ec

−8.04

−8.06

1.07 (13.3%)

−8.39

0.59 (7.0%)

eσq

2.37

2.37

0.16 (6.9%)

2.49

0.18 (7.4%)

eσl

0.23

0.22

0.10 (45.5%)

0.27

0.037 (13.6%)

eσ As

1.29

1.28

0.16 (12.6%)

1.38

0.10 (7.2%)

eσ Ac

1.31

1.30

0.22 (16.5%)

1.47

0.13 (8.9%)

eσ Es

1.28

1.29

0.16(12.6%)

1.39

0.098 (7.1%)

eσ Ec

0.26

0.25

0.10 (41.8%)

0.30

0.037 (12.2)%

g-eval ∇θ g -eval

1000 (+FORM) -

3000 800

5.3. Elastoplastic frame

Table 6 Uncertain parameters of Example 3.

This example is an elastoplastic frame structure subjected to a horizontal load (Fig. 5). This example is taken from [43] and was studied in the context of reliability sensitivity analysis in [30]. The frame has four potential failure modes represented by the following LSFs:

Parameter

Distribution

Mean

CV

Mi (i = 1, 2, 3) [tm] S [tm]

Lognormal Lognormal

200 50

15% 40%

g31 (M1, M2, M3, S )=2M1 + 2M3−4.5S chosen in the vicinity of the FORM design points of the individual failure modes and the standard deviations were chosen such that the IS density has adequate spread to reach all four failure regions. We note that the choice of the standard deviations requires considerable knowledge of the structure of the joint failure domain. The reference solution was computed with IS using 5 × 105 samples. In SIS, we use 1000 samples per level and consider two cases for the target coefficient of variation of the weights, δtarget = 1.5 and 1.0 . The mean estimates for the IS approach with 2000 samples show good agreement with the reference values. The CV of the probability estimate is acceptable which shows that the chosen IS density is effective in evaluating the probability of system failure. However, as with Example 2 we observe inconsistent CVs of the partial elasticity estimates. A special case is the elasticity to the mean of the load S, e μS : The value of e μS is close to zero, which indicates that a small change in the mean of S will have negligible influence on the probability of failure. For this case, the CV is not an appropriate measure of the accuracy of the simulation. The mean estimates of the elasticities obtained with SIS and δtarget = 1.5 are close to the reference values, except from e μS for which a significant bias is observed. The latter is due to the fact that the σ requires that the CV of adaptive approach for selecting the parameter ∼ the weights is as close as possible to δtarget in the allowable range (0, σL) , where σL is the parameter of the final sampling density in SIS. Choosing σ would lead to a smaller bias but would also increase the CV a smaller ∼ of the estimate. Through decreasing δtarget to 1.0 , the value of σL decreases significantly and hence the bias in the estimate of the elasticities essentially disappears. As with Example 2, the CV of the elasticity estimates obtained with SIS is consistent with the one of the probability of failure, except from the special case of e μS .

g32 (M1, M2, M3, S )=2M1 + M2 + M3−4.5S g33 (M1, M2, M3, S )=M1 + M2 + 2M3−4.5S g34 (M1, M2, M3, S )=M1 + 2M2 + M3−4.5S The generalized LSF representing failure of the system is expressed as follows

g3,sys (M1, M2, M3, S ) = min {g3i (M1, M2, M3, S )} i = 1,2,3,4

The plastic moments Mi , i = 1, 2, 3, and the load S are modeled by independent lognormal random variables with means and CVs given in Table 6. The statistics of the estimated probability of failure and partial elasticities of the probability in terms of the means and standard deviations of the random variables are given in Table 7. We compare the performance of SIS with the one of the standard IS method in original space. For the latter, we use as IS function an independent normal distribution with means [185;185;185;160] and standard deviations set to twice the standard deviations of the random variables. The means were

S

6m M2 M1

M3

4.5m

6. Concluding remarks This paper presented a SIS method for estimating sensitivities of the probability of failure to deterministic parameters. SIS estimates the probability of failure through sampling a sequence of distributions that

Fig. 5. Elastoplastic frame according to [43]. 31

Structural Safety 75 (2018) 24–34

I. Papaioannou et al.

Table 7 Results for limit state function of Example 3. Note that the IS solution requires prior knowledge of the FORM design points. Reference

IS (2000 samples)

SIS (1000 samples/ level, δtarget = 1.5 )

SIS (1000 samples/ level, δtarget = 1.0 )

Mean estimate

St. dev. (CV)

Mean estimate

St. dev. (CV)

Mean estimate

St. dev. (CV)

5.33 × 10−4 −3.92

5.30 × 10−4 −3.89

3.97 × 10−5 (7.5%) 0.47 (12.0%)

5.19 × 10−4 −4.15

6.10 × 10−5 (11.5%) 0.29 (7.1%)

5.27 × 10−4 −4.05

4.53 × 10−5 (8.6%) 0.36 (9.0%)

eμM

−2.13

−2.20

0.51 (23.1%)

−2.59

0.18 (7.0%)

−2.25

0.28 (12.4%)

eμM

−3.92

−3.90

0.51 (13.0%)

−4.15

0.30 (7.2%)

−4.05

0.36 (9.0%)

e μS

−0.02

−0.02

0.086 (430.1%)

1.09

0.18 (18.1%)

0.01

0.063 (630.5%)

eσ M1

0.38

0.38

0.10 (27.4%)

0.57

0.060 (10.6%)

0.39

0.061 (15.7%)

eσ M2

0.30

0.31

0.10 (32.5%)

0.44

0.051 (11.7%)

0.31

0.052 (16.8%)

eσ M3

0.38

0.38

0.10 (26.7%)

0.57

0.059 (10.4%)

0.39

0.059 (15.2%)

eσS

8.95

8.96

0.15 (1.7%)

9.21

0.65 (7.1%)

9.16

0.76 (8.3%)

Pf

eμM

1 2

3

g-eval ∇θ g -eval

2000 (+FORM) –

4000 820

6000 686

n(y) d

g(x, ) = 0

n(y)

Fig. 6. Illustration of the integration domain in the approximation of Eq. (A.3).

able to produce accurate estimates of the sensitivities with consistent CVs with the ones obtained for the probability of failure. Finally, it should be noted that the proposed method requires evaluation of the gradient of the limit-state function to the deterministic parameters at each sample of the final sampling step. In problems involving limit-state functions expressed in terms of computationally expensive numerical models this task can involve significant additional computational cost. One possibility to reduce the computational cost is to reduce the number of gradient evaluations by accounting for the fact that the contribution of the samples away from the limit-state surface is insignificant. Another possibility is to combine the proposed approach with a surrogate model of the limit-state function based on the samples of the initial SIS run thus circumventing the gradient evaluations with the numerical model.

gradually approach the optimal IS density. Using the samples at the final sampling level, the proposed method estimates a domain integral approximation of the sensitivity of the probability of failure. The σ of the approximation is determined adaptively to ensure parameter ∼ consistency of the uncertainty of the sensitivity estimate with the one of the probability of failure. The method when implemented in the equivalent standard normal space can be applied for estimating sensitivities to any type of deterministic parameters, including parameters of the distribution of the random variables. We note that implementing the SIS approach in standard normal space has the advantage that tailored MCMC algorithms can be used for generating the samples at the intermediate levels. Moreover, the method is applicable to the estimation of reliability sensitivities of system problems with multiple failure modes. Numerical examples demonstrated that the proposed method is Appendix A

In this Appendix we show that under certain regularity conditions, the approximation of Eq. (10) converges to the true surface integral of Eq. (5) σ → 0 . For the sake of simplicity, we look at the case where the parameter vector θ reduces to a scalar θ representing a limit-state parameter. Eq. as ∼ (10) then reads: ∼ ∂Pf (θ , ∼ σ) g (x , θ) ∂g (x , θ) 1 f (x ) dx = − n ∼φ⎛ ∼ ⎞  σ ∂θ (A.1) ⎝ σ ⎠ ∂θ

∫

The true value of the derivative of Pf with respect to θ is:

dPf (θ) dθ

=−

∫g (x,θ)=0 f (x ) ‖∇x g (1x , θ)‖

∂g (x , θ) ds (x ) ∂θ

(A.2)

32

Structural Safety 75 (2018) 24–34

I. Papaioannou et al.

where ds (x ) denotes surface integration over the surface {g (x , θ) = 0} . We make the following assumptions: (a) The function g (x , θ) is continuously differentiable. (b) In a neighborhood of the surface g (x , θ) = 0 , one can make a one-to-one transformation x → (y, h) , where y is the coordinate on the surface and h is the distance from the point y along the surface normal n (y ) = ‖∇y g (y, θ)‖−1 ∇y g (y, θ) . Under assumption (b), we change variables in the integral of Eq. (A.1) by applying the mapping x → (y, h) . We further truncate the integration domain at a distance d to the surface g (x , θ) = 0 , leading to the following approximation: ∼ d 1 ∂Pf (θ , ∼ σ) g (y + hn (y ), θ) ⎞ ∂g (y + hn (y ), θ) −1 J (y + hn (y )) dhf (y ) ds (y ) ≈− φ⎛ ∼ g (y, θ ) = 0 −d ∼ σ ∂θ σ ⎝ ∂θ (A.3) ⎠

∫

∫

where J −1 (.) is the Jacobian of the inverse mapping (y, h) → x . The truncated integration domain in Eq. (A.3) is illustrated in Fig. 6. Since n (y ) is orthogonal to the limit-state surface at y one can always parametrize the surface such that J −1 (y, 0) = 1. Therefore, we have ∼ ∂Pf (θ , ∼ σ) ⎡ d 1 φ ⎛ g (y + hn (y ), θ) ⎞ dh⎤ ∂g (y, θ) f (y ) ds (y ) ≈− ∼ g (y, θ ) = 0 ⎢ −d ∼ σ ∂θ σ ⎝ (A.4) ⎠ ⎥ ⎣ ⎦ ∂θ

∫

∫

The integral over h can be approximated making a first-order Taylor expansion of g at y , which gives

g (y + hn (y ), θ) ⎞ dh ≈ ∼ σ ⎝ ⎠

∫−d ∼σ1 φ ⎛ d

g (y, θ) + h‖∇y g (y, θ)‖ ⎞ dh ∼ σ ⎝ ⎠

∫−d ∼σ1 φ ⎛ d

⎜

⎟

(A.5)

The above follows from the first order multivariate Taylor expansion:

g (y + hn (y ))−g (y, θ) ≈ hn (y )T ∇y g (y, θ) = h‖∇y g (y, θ)‖

(A.6)

Since y lies on the limit-state surface, one has g (y, θ) = 0 in Eq. (A.5). The resulting integral can be written as

h‖∇y g (y, θ)‖ ⎞ 1 ⎡Φ ⎛ d‖∇y g (y, θ)‖ ⎞−Φ ⎛− d‖∇y g (y, θ)‖ ⎞ ⎤ dh = ∼ ∼ ∼ ⎢ ⎥ σ ‖ g ( , θ )‖ σ σ y ∇ y ⎠⎦ ⎝ ⎠ ⎠ ⎝ ⎣ ⎝

∫−d ∼σ1 φ ⎛ d

⎜

⎟

⎜

⎟

⎜

⎟

σ → 0 . Inserting this into Eq. (A.3) shows that The term in brackets in Eq. (A.7) converges towards unity as ∼

(A.7) ∼ ∂Pf (θ, ∼ σ) ∂θ

converges to

dPf (θ ) dθ

σ → 0. as ∼

[20] Au S. Reliability-based design sensitivity by efficient simulation. Comput Struct 2005;83:1048–61. [21] Luyi L, Zhenzhou L, Jun F, Bintuan W. Moment-independent importance measure of basic variable and its state dependent parameter solution. Struct Saf 2012;38:40–7. [22] Breitung K. Asymptotic approximations of probability integrals. Lect Notes Math vol. 1592. Berlin: Springer; 1994. [23] Wu Y-T. Computational methods for efficient structural reliability and reliability sensitivity analysis. AIAA J 1994;32:1717–23. [24] Song S, Lu Z, Qiao H. Subset simulation for structural reliability sensitivity analysis. Reliab Eng Syst Safe 2009;94:658–65. [25] Jensen H, Mayorga F, Valdebenito M. Reliability sensitivity estimation of nonlinear structural systems under stochastic excitation: a simulation-based approach. Comput Methods Appl Mech Eng 2015;289:1–23. [26] Dubourg V, Sudret B. Meta-model-based importance sampling for reliability sensitivity analysis. Struct Saf 2014;49:27–36. [27] Au S-K, Patelli E. Rare event simulation in finite-infinite dimensional space. Reliab Eng Syst Saf 2016;148:67–77. [28] Bjerager P. Probability integration by directional simulation. J Eng Mech-ASCE 1988;114:1285–302. [29] Ditlevsen O, Bjerager P, Olesen R, Hasofer A. Directional simulation in Gaussian processes. Prob Eng Mech 1988;3:207–17. [30] Lu Z, Song S, Yue Z, Wang J. Reliability sensitivity method by line sampling. Struct Saf 2008;30:517–32. [31] Papaioannou I, Breitung K, Straub D. Reliability sensitivity analysis with Monte Carlo methods. Deodatis G, Ellingwood BR, Frangopol DM, editors. Proc. 11th International Conference on Structural Safety and Reliability (ICOSSAR) New York: Columbia University; 2013. [32] Valdebenito M, Jensen H, Hernández H, Mehrez L. Sensitivity estimation of failure probability applying line sampling. Reliab Eng Syst Saf 2018;171:99–111. [33] Royset JO, Polak E. Reliability-based optimal design using sample average approximations. Prob Eng Mech 2004;19:331–43. [34] Jensen HA, Valdebenito MA, Schuëller GI, Kusanovic DS. Reliability-based optimization of stochastic systems using line search. Comput Method Appl Mech Eng 2009;198:3915–24. [35] Lacaze S, Brevault L, Missoum S, Balesdent M. Probability of failure sensitivity with respect to decision variables. Struct Multidisciplinary Opt 2015;52:375–81. [36] Torii AJ, Lopez RH, Miguel LFF. Probability of failure sensitivity analysis using polynomial expansion. Prob Eng Mech 2017;48:76–84. [37] Rosenblatt M. Remarks on a multivariate transformation. Ann Math Stat 1952;23:470–2. [38] Spanier J, Oldham KB. An atlas of functions. Washington, DC: Hemisphere; 1987. [39] Doucet A, De Freitas N, Gordon NJ. Sequential Monte Carlo methods in practice. New York: Springer-Verlag; 2001.

References [1] Ditlevsen O, Madsen HO. Structural reliability methods. Chichester: Wiley; 1996. [2] Lemaire M. Structural reliability. Hoboken, NJ: Wiley; 2009. [3] Der Kiureghian A. First-and second-order reliability methods. In: Nikolaidis E, Ghiocel DM, Singhal S, editors. Engineering Design Reliability Handbook. Boca Raton, FL: CRC Press; 2005. [4] Sudret B. Meta-models for structural reliability and uncertainty quantification. Phoon KK, Beer M, Quek ST, Pang SD, editors. Proc. Fifth Asian-Pacific Symposium on Structural Reliability and its Applications (5APSSRA), Singapore, 2012 2012. [5] Rackwitz R. Reliability analysis – a review and some perspectives. Struct Saf 2001;23:365–95. [6] Valdebenito MA, Pradlwarter HJ, Schuëller GI. The role of the design point for calculating failure probabilities in view of dimensionality and structural nonlinearities. Struct Saf 2010;32:101–11. [7] Breitung K. Extrapolation, invariance, geometry and subset sampling. In: Caspeele R, Taerwe L, Proske D, editors. 14th International Probabilistic Workshop. Cham, CH: Springer; 2017. p. 33–43. [8] Breitung K, The geometry of limit state function graphs and subset simulation, arXiv preprint arXiv:1705.04453 (2017b). [9] Bucher C. Adaptive sampling - an iterative fast Monte Carlo procedure. Struct Saf 1988;5:119–26. [10] Au SK, Beck JL. A new adaptive importance sampling scheme. Struct Saf 1999;21:135–58. [11] Kurtz N, Song J. Cross-entropy-based adaptive importance sampling using Gaussian mixture. Struct Saf 2013;42:35–44. [12] Dubourg V, Sudret B, Deheeger F. Metamodel-based importance sampling for structural reliability analysis. Prob Eng Mech 2013;33:47–57. [13] Hohenbichler M, Rackwitz R. Improvements of second-order reliability estimates by importance sampling. J Eng Mech-ASCE 1988;114:2195–9. [14] Koutsourelakis PS, Pradlwarter HJ, Schuëller GI. Reliability of structures in high dimensions, part I: algorithms and applications. Prob Eng Mech 2004;19:409–17. [15] Au SK, Beck JL. Estimation of small failure probabilities in high dimensions by subset simulation. Prob Eng Mech 2001;15:263–77. [16] Papaioannou I, Betz W, Zwirglmaier K, Straub D. MCMC algorithms for subset simulation. Prob Eng Mech 2015;41:89–103. [17] Beaurepaire P, Jensen HA, Schuëller GI, Valdebenito MA. Reliability-based optimization using bridge importance sampling. Prob Eng Mech 2013;34:48–57. [18] Papaioannou I, Papadimitriou C, Straub D. Sequential importance sampling for structural reliability analysis. Struct Saf 2016;62:66–75. [19] Melchers RE, Ahammed M. A fast approximate method for parameter sensitivity estimation in Monte Carlo structural reliability. Comput Struct 2004;82:55–61.

33

Structural Safety 75 (2018) 24–34

I. Papaioannou et al. [40] Latz J, Papaioannou I, Ullmann E. Multilevel Sequential2 Monte Carlo for Bayesian inverse problems. J Comput Phys 2018;368:154–78. [41] Jasra A, Stephens DA, Doucet A, Tsagaris T. Inference for Lévy-driven stochastic volatility models via adaptive sequential Monte Carlo. Scandinavian J Stat 2011;38:1–22.

[42] Yang RJ, Gu L. Experience with approximate reliability-based optimization methods. Struct Multidisc Opt 2004;26:152–9. [43] Zhao Y-G, Ang AH. System reliability assessment by method of moments. J Struct Eng 2003;129:1341–9.

34