Scarce Sample-Based Reliability Estimation and Optimization Using Importance Sampling

Abstract

Importance sampling is a variance reduction technique that is used to improve the efficiency of Monte Carlo estimation. Importance sampling uses the trick of sampling from a distribution, which is located around the zone of interest of the primary distribution thereby reducing the number of realizations required for an estimate. In the context of reliability-based structural design, the limit state is usually separable and is of the form Capacity (C)–Response (R). The zone of interest for importance sampling is observed to be the region where these distributions overlap each other. However, often the distribution information of C and R themselves are not known, and one has only scarce realizations of them. In this work, we propose approximating the probability density function and the cumulative distribution function using kernel functions and employ these approximations to find the parameters of the importance sampling density (ISD) to eventually estimate the reliability. In the proposed approach, in addition to ISD parameters, the approximations also played a critical role in affecting the accuracy of the probability estimates. We assume an ISD which follows a normal distribution whose mean is defined by the most probable point (MPP) of failure, and the standard deviation is empirically chosen such that most of the importance sample realizations lie within the means of R and C. Since the probability estimate depends on the approximation, which in turn depends on the underlying samples, we use bootstrap to quantify the variation associated with the low failure probability estimate. The method is investigated with different tailed distributions of R and C. Based on the observations, a modified Hill estimator is utilized to address scenarios with heavy-tailed distributions where the distribution approximations perform poorly. The proposed approach is tested on benchmark reliability examples and along with surrogate modeling techniques is implemented on four reliability-based design optimization examples of which one is a multi-objective optimization problem.

1. Introduction

Reliability-based design optimization (RBDO) is a design approach that is used to generate reliable designs by accounting for uncertainties in the system variables such as material properties, loading conditions and geometry. Probabilistic approaches use probability distributions to model the uncertainties. In such approaches, the reliability of the system is measured as the probability of failure to satisfy a performance criterion. Mathematically, this involves calculating the hyper-volume of a multi-dimensional probability density function (PDF) under the failure region. This calculation becomes infeasible when the number of dimensions are large. Even in low dimensions, this calculation could become intractable due to complicated geometry of the failure region [1]. Attractive alternatives are analytical methods or sampling-based approaches.

Analytical methods, such as first-order and second-order reliability methods (FORM and SORM), transform the probability distributions to standard normal space and use linear or quadratic approximations of the performance function to estimate the failure probability at the most probable point (MPP) of failure [2,3,4]. In essence, analytical methods can estimate failure probability within a reasonable number of evaluations for linear or slightly non-linear performance functions. However, if the failure boundary is highly non-linear, analytical approaches are likely to lead to erroneous estimation. Additionally, other factors, such as island failure regions and multiple failure modes, limit their performance [5,6].

Sampling-based approaches such as Monte Carlo methods are highly effective for complex failure boundary problems and multi-modal failure problems [7]. In high-reliability applications where the failure probability is very low, MCS requires a very large number of model evaluations to obtain an accurate estimate. Most RBDO applications involve evaluating high-fidelity models which are computationally expensive, thus rendering MCS prohibitive for reliability estimation.

The computational burden of RBDO can be reduced by using surrogate-based methods wherein surrogate models which are cheaper to evaluate are constructed using limited high-fidelity model evaluations. Uncertainty from the random variables is then propagated through these surrogate models to obtain reliability estimates. Various surrogate modeling approaches, such as polynomial response surface (PRS), radial basis function (RBF), support vector machine (SVM) and Kriging among others have been adopted for reliability estimation [8]. Li and Xiu [9] proposed using cheaper to evaluate surrogates away from the limit state and high-fidelity model evaluations close to the limit state to improve accuracy and reduce cost. Dai et al. [10] proposed an SVM-based radial basis function to approximate the limit state function which is then used to estimate failure probability. Dubourg et al. [11] used error measure derived from Kriging to refine the surrogate during subset-simulation-based reliability estimation. Reliability estimation using surrogates may carry forward the bias from the surrogate approximation [12]. Surrogate modeling can also be used to build approximation models for reliability metrics instead of limit state function. Foschi et al. [13] used the combination of response surface, FORM and importance sampling to perform reliability estimation. Qu and Haftka [14] compared the accuracy of surrogates of different reliability metrics, such as failure probability, reliability index and probabilistic sufficiency factor (PSF) during RBDO. They conclude that inverse measure, such as PSF operating in performance space, performs better compared to a classical reliability metric. A survey of various surrogate-based RBDO frameworks is provided in [15]. The choice between surrogate for reliability estimation versus surrogate for limit state is a trade-off based on the number of design variables and the number of random variables as well as the cost of building the surrogate and cost of reliability estimate [16].

Employing variance reduction techniques to improve the efficiency and accuracy of Monte Carlo estimations is another way to alleviate the computational burden of RBDO. Several variance reduction techniques, such as importance sampling (IS) [17,18,19], subset simulation (SS) [1,20] and separable Monte Carlo (SMC) [21,22], are used to improve the accuracy of failure probability estimates while reducing the required number of high-fidelity model evaluations.

Importance sampling improves the accuracy of the estimate by drawing the sample realizations of the input random variables that have greater impact on the estimate, more often. In order to do that, an alternate sampling density known as importance sampling density (ISD) is chosen which enables sampling the important values more frequently. This introduces a bias in the estimator which is corrected by weighing the sample realizations. A good choice of ISD is proven to improve accuracy and in contrast, incorrect choice of ISD could lead to spurious estimates. Hence, the choice of ISD is critical, and several approaches have been explored to find the optimal ISD. Melchers [17] applied importance sampling for assessing reliability of parallel and series structural systems, where a multi-normal PDF centered at MPP was chosen as the ISD. Using the original distribution shifted to MPP, multivariate normal distribution located at MPP with various choices for the co-variance matrix ranging from the same as the original distribution to using only diagonal elements of the original co-variance matrix have been studied [23].

In cross-entropy-based methods [24], ISD is found by minimizing the Kullback–Leibler (KL) divergence between theoretical optimal ISD and a chosen family of distributions. Kurtz and Song [25] used a Gaussian mixture to obtain a non-parametric multi-modal PDF for ISD. It was observed that the coefficient of variation of the failure probability estimate converged as the number of Gaussian densities in the mixture increased. Cao and Choe [26] used an expectation-maximization (EM) algorithm to obtain a Gaussian mixture as the near optimal ISD. Cross-entropy information criterion (CIC) was used to select the number of Gaussian densities in the mixture. Geyer et al. [27] proposed a modified version of EM algorithm for updating the Gaussian-mixture-based ISD. For selecting the number of distributions in the mixture, density-based spatial clustering of applications with noise (DBSCAN) algorithm was used. In cross-entropy-based IS, one still requires the joint PDF of the original random variables to compute the near optimal ISD.

Kernel-based IS is another way to obtain the alternate sampling density wherein a kernel sampling density is constructed instead of choosing from a family of distributions [28,29]. Au and Beck [30] proposed Markov Chain Monte Carlo (MCMC) to distribute samples asymptotically according to the optimal ISD, and subsequently a kernel sampling density is constructed from these samples. Dai et al. [31] proposed a wavelet density estimation technique to construct the ISD from the MCMC samples. Botev et al. [32] proposed combining MCMC and IS to address the issues of biased sampling estimators that result from MCMC. Various adaptive importance sampling procedures have gained popularity recently. Dalbey and Swiler [33] proposed a Gaussian-process-based adaptive importance sampling where a Gaussian process surrogate is used to identify the likely regions of failure to adaptively improve the sampling density estimate. Zhao et al. [34] constructed a Kriging surrogate from the initial MCMC samples which is improved using an active learning process, and subsequently the surrogate is used for limit state evaluations for adaptive improvement of ISD. Wang and Song [35] used cross-entropy-based adaptive importance sampling for high-dimensional reliability analysis. Here, ISD is obtained by minimizing the KL divergence between a von Mises–Fisher mixture model and near optimal ISD.

In all the literature discussed above, while employing IS, determination of the appropriate ISD is considered to be the challenge. The probabilistic distributions and the corresponding parameters of the original random variables, however, are assumed to be known [22,36]. However, this is not necessarily always true. Here, we propose a framework to employ the importance sampling method for separable failure boundary of the form Capacity (C)–Response (R) using only scarce realizations of the R and C where no information about their distributions is known. We make use of tail-modeling techniques to approximate the cumulative distribution function (CDF) of the capacity and response. These approximations are used to locate a Gaussian ISD whose standard deviation is empirically chosen. Sample realizations drawn from the ISD are used to compute the failure probability. During this computation, kernel density estimates of PDF of response are utilized along with CDF approximation of capacity. Furthermore, bootstrap samples are generated from the original scarce samples of R and C. The proposed framework is applied on the bootstrap samples to obtain the confidence bounds for the reliability estimate. Reliability estimates obtained from scarce samples using the proposed approach are used to construct a surrogate of reliability index which is used for constraint evaluation during RBDO.

The rest of the paper is organized as follows. Section 2 presents the theoretical background of importance sampling in the context of reliability estimation. In Section 3, the methodology is discussed, and in later sections we present the results from the proposed method when applied to test cases and RBDO applications including a multi-objective optimization (MOO) application. Appendix A and Appendix B are used to provide short descriptions of kernel density estimation (KDE) and the third-order polynomial normal transformation (TPNT) technique.

2. Reliability Estimation Using Importance Sampling for Separable Limit States

In structural engineering, limit states are useful in prescribing performance requirements of a design. Thus, a limit state decomposes the design space of a system into safe and failure regions. Violation of a limit state is considered to be failure. and reliability is a measure of probability of such violations. For most structural problems, limit state can be expressed as the difference of Capacity $\left(C\right)$ and Response $\left(R\right)$, as presented in Equation (1), where R and C are functions of independent sets of random variables. This is referred to as a separable limit state [21,22,37], and we consider such limit states in this work.

$$G(C,R)=C-R$$

(1)

The system is said to fail when $C\le R$ and safe when $C>R$. When either capacity or response or both are functions of variables that are random, a probabilistic measure such as the probability of failure is as presented in Equation (2)

$$p_f=\underset{G\le 0}{\int \int }{f}_{CR}(c,r)dcdr$$

(2)

where $f_{CR}$ is the joint probability distribution function (PDF), and $G\le 0$ is the failure region. In the case of separable limit state, the joint PDF of capacity and response can be decomposed into the product of the marginal PDFs of capacity and response as presented in Equation (3)

$$p_f=\underset{-\infty }{\overset{\infty }{\int }}\underset{-\infty }{\overset{c\le r}{\int }}{f}_C\left(c\right)f_R\left(r\right)dcdr$$

(3)

$$p_f=\underset{-\infty }{\overset{\infty }{\int }}{F}_C\left(x\right)f_R\left(x\right)dx$$

(4)

where $f_C\left(c\right)$ and $f_R\left(r\right)$ are the marginal density functions of C and R, respectively. Equation (3) can also be written in a single integral form as presented in Equation (4). Here, $F_C\left(c\right)$ is the cumulative distribution function (CDF) of capacity. For low failure probabilities of order ${10}^{-4}$, to obtain a Monte Carlo estimate with $10\%$ coefficient of variation, the sample size required is $100\times \frac{1}{p_f}={10}^6$. Generating ${10}^6$ instances of expensive computer models is not feasible. To reduce such computational burden, importance sampling is used as in Equation (5).

$$p_{f_{IS}}=\underset{-\infty }{\overset{\infty }{\int }}\frac{F_C\left(x\right)f_R\left(x\right)}{h_x\left(x\right)}{h}_x\left(x\right)dx$$

(5)

where $h_x\left(x\right)$ is the alternate sampling density, and $p_{f_{IS}}$ is the failure probability computed using the importance sampling approach. It is the expectation of the integrand $\frac{F_C\left(x\right)f_R\left(x\right)}{h_x\left(x\right)}$ computed with respect to the ISD, $h_x$. Consequently, if $x_1,x_2,\dots ,x_N$ is an independent identically distributed (i.i.d.) random sample from $h_x$, then the expectation has an unbiased estimator:

$$\widehat{p_{f_{IS}}}=\sum _{i=1}^N\frac{F_C\left(x_i\right)f_R\left(x_i\right)}{h_x\left(x_i\right)},x_i\sim h_x$$

(6)

Alternatively, if the marginal distribution of the response $\left(f_R\left(r\right)\right)$ is integrated in Equation (3), then the importance sampling estimator of $p_f$ would be expressed as:

$$\widehat{p_{f_{IS}}}=\sum _{i=1}^N\frac{(1-F_R\left(x_i\right))f_C\left(x_i\right)}{h_x\left(x_i\right)},x_i\sim h_x$$

(7)

As mentioned in Section 1, in most of the literature, importance sampling is used when the distributions of R and C are known, and it becomes computationally expensive to draw samples in the region where tails of R and C overlap. Figure 1 presents a schematic of such a scenario where the limit state is of the separable form. Here, the region of interest lies in the tails of the distributions of R and C, and it contains the most probable point (MPP) of failure. Hence, it is only logical to locate the ISD at the MPP [3,17,22].

Figure 1 depicts a Gaussian ISD for which MPP serves as the mean of the distribution and is usually found by solving a constrained optimization problem. Once the mean of the ISD is determined, the spread of the distribution needs to be chosen. In order to have a finite variance for the importance sampling estimator, the ISD should not have a lighter tail than the original distribution [38]. In the literature [22,23], using the same co-variance matrix and sometimes using strictly the diagonal elements of the co-variance matrix of the original distribution is suggested. This has been shown to work for a range of possible limit states. However, sometimes a designer does not necessarily know the original distributions of input variables but has only a few realizations of R and C either through expensive computer simulations or physical experiments. Generally, response is considered to be the source of randomness since it is a measured quantity. However, capacity can also be a random variable. For instance, an example of capacity is yield strength which is a measured quantity. We know that it is random and hence modeled using probability distributions [39,40]. Other examples of capacity that are considered to be random include: member capacity under seismic loading [41], maximum flow capacity in measuring hydraulic reliability of water distribution system [42], and material fatigue properties [43,44].

In this work, we first characterize the distributions of R and C from limited samples and use the concept of importance sampling to estimate low failure probabilities. A Gaussian distribution is chosen as the ISD which is fully described by its two parameters, the mean and variance. Section 3 describes the procedure employed to identify the parameters of ISD. In the proposed framework, to demonstrate the methodology, distribution information of original input random variables is used only to generate the initial limited samples. However, the distribution information is not utilized anywhere afterwards.

3. Identifying Parameters of Gaussian ISD

Figure 2 presents an example scenario [23] of a separable limit state where the response and capacity follow normal distributions. It can be observed that the functional $F_C{f}_R$ integrated in Equation (4) is maximum at the point $(x\sim 11.55)$, near the point of intersection between the PDF of response and the CDF of capacity, i.e., $f_R\left(x\right)=F_C\left(x\right)$. In importance sampling, MPP is used to locate the ISD to maximize sampling of the failure probability content. Hence, it is only logical to define this point of intersection as the MPP and use it as the mean of Gaussian ISD.

This point can be anywhere within the bounds of capacity and response. In high-reliability scenarios, it is located at the tails of one of the distributions or both, so we need to extrapolate into the tails of R and C. The available PDF approximation methods capture the central part of the distribution better than the tails, whereas the main focus of the tail modeling techniques applied for CDF estimation is at the tails of the distribution. Hence, we can reduce the errors in finding the MPP by using the intersection of $1-F_R$ and $F_C$ instead of the intersection point of the curves $f_R$ and $F_C$. Here, $1-F_R$ is the complementary CDF of R. It is to be noted that the suitability of such a modification has been tested for only uni-modal type response. For multi-modal response, using a single MPP will only capture partial failure region. However, the approach can be extended by using a mixture of Gaussian distributions by using a summed failure probability integral as described in [45]. Algorithm 1 presents the bracketing procedure followed in the current work to find the MPP, while the same process is visually presented in Figure 3.

In Steps 3 and 6 of Algorithm 1, it is stated that the CDF approximations of capacity ($\widehat{F_C}$) and response ($\widehat{F_R}$) are obtained through a TPNT technique (details provided in Appendix B); however, it is an independent block in the proposed approach and hence can be replaced by any other suitable technique based on user discretion. In effect, the MPP finding problem reduces to finding the zeroes of the function $\widehat{F_C}\left(x\right)-(1-\widehat{F_R}\left(x\right))$. Hence, root finding algorithms can also be applied to find the MPP. However, it is advised to use derivative-free approaches to avoid numerical issues.

Algorithm 1 Finding ${\mu }_h$.

1: Obtain samples of response $\mathbf{r}=\{r_1,r_2,\dots ,r_M\}$ and capacity $\mathbf{c}=\{c_1,c_2,\dots ,c_N\}$ as in Figure 3a. 2: for$i=1,2,\dots ,M$do 3:     $\widehat{F_C}\left(r_i\right)$, $\widehat{F_R}\left(r_i\right)\leftarrow$ TPNT approximations of CDF of response and capacity at response sample $\mathbf{r}$. 4:     Obtain point A in Figure 3b ← maximum of ${\mathbf{r}}_{\mathbf{x}}=\left\{r_i\right|\widehat{F_C}\left(r_i\right)-(1-\widehat{F_R}\left(r_i\right))<0\}$. 5: for$i=1,2,\dots ,N$do 6:     $\widehat{F_C}\left(c_i\right)$, $\widehat{F_R}\left(c_i\right)\leftarrow$ TPNT approximations of CDF of response and capacity at capacity sample $\mathbf{c}$. 7:     Obtain point B in Figure 3c ← minimum of ${\mathbf{c}}_{\mathbf{x}}=\left\{c_i\right|\widehat{F_C}\left(c_i\right)-(1-\widehat{F_R}\left(c_i\right))>0\}$. 8: Generate a set of evenly spaced points $\mathbf{x}=\{x_i:i=1,2,\dots ,100\}$ between points A and B. 9: Obtain mean of ISD, ${\mu }_h\leftarrow$ maximum of $\mathbf{x}=\left\{x_i\right|\widehat{F_C}\left(x_i\right)-(1-\widehat{F_R}\left(x_i\right))<0\}$. (Alternatively, minimum of $\mathbf{x}=\left\{x_i\right|\widehat{F_C}\left(x_i\right)-(1-\widehat{F_R}\left(x_i\right))>0\}$ can also be used.)

In very low failure probability estimation, the region corresponding to the estimation is quite small, and it has been observed that aspects such as non-uniqueness of MPP and non-linearity of the limit states have little effect on the estimation [5]. Similarly, we find that errors in the MPP estimation have lesser effect on the reliability estimation compared to a poor choice of standard deviation.

The standard deviation of the ISD determines the importance ascribed to region around the MPP during sampling. The region of sampling itself could be bounded by the supports of capacity and response distributions. Different measures of spread could be applied based on the available knowledge. In this work, a 10% coefficient of variation (CoV) was used to calculate the standard deviation of the ISD to investigate the effect. Such a measure was found to be suitable for scarce samples of normally distributed capacity and response. Different measures of spread were investigated while testing the formulation on samples simulated to be belonging to a different distribution. The spread parameter obtained through Equation (8) was found to be appropriate as such a measure allows one to restrict 68% percent of importance sample realizations within the means of capacity and response.

$${\sigma }_h=\frac{{\overline{x}}_C-{\overline{x}}_R}{2}$$

(8)

where ${\overline{x}}_C$ is the sample mean of capacity, and ${\overline{x}}_R$ is the sample mean of the response.

4. Estimation of Reliability and Its Confidence Bounds

$$\widehat{p_f}=\sum _{i=1}^N\frac{\widehat{F_C}\left(x_i\right)\widehat{f_R}\left(x_i\right)}{h_x\left(x_i\right)},x_i\sim h_x:\mathcal{N}({\mu }_h,{\sigma }_h^2)$$

(9)

Using the now defined ISD, sample points are drawn at which the PDF of response and CDF of capacity are obtained to be used in the computation of failure probability as per Equation (6). However, instead of the actual values, approximations $\widehat{F_C}\left(x_i\right)$, $\widehat{f_R}\left(x_i\right)$ are used during the computation as presented in Equation (9). As mentioned earlier, TPNT is used to approximate CDF ($\widehat{F_C}$), whereas for PDF approximation ($\widehat{f_R}$) an adaptive kernel density estimation method is employed. Both of these methods are distribution-free methods; however, other suitable methods of approximation [46,47] can also be used. Despite the accuracy of the method chosen, errors from the approximations are bound to result in loss of accuracy in the failure probability estimate. Hence, it becomes necessary to quantify the confidence on the estimate. Here, confidence bounds on the estimate are computed using a non-parametric bootstrap method.

The underlying idea of a non-parametric bootstrap method is to recreate samples from the original sample by sampling with the replacement. The sample size of bootstrap samples must be the same as the original sample. From each of the bootstrap samples, a statistical quantity of interest (such as failure probability) can be estimated. By repeating the process many times (say B times), one can obtain a distribution around the estimate from the original sample. In the current work, the proposed approach is applied to the bootstrap samples of C and R to obtain confidence bounds on the reliability estimate. Since the original samples of C and R are themselves scarce, this process of bootstrapping is repeated for T ($=100$) original samples of C and R. Thus, the quantiles of $\widehat{\beta }$ obtained from T ($=100$) iterations are compared with mean and standard deviation of quantiles of ${\widehat{\beta }}_{boot}$ from the bootstrap samples. Figure 4 presents a schematic representation of the bootstrapping procedure followed in the current work.

In order to test the efficacy of the formulation, different tails for capacity and response were considered. To simulate samples of response and capacity belonging to different tail types, generalized extreme value (GEV) distribution is used which takes the shape parameter $\left(\xi \right)$ as an input along with location $\left(\theta \right)$ and scale $\left(\sigma \right)$ parameters. Shape parameter or tail index is a measure of heaviness of tails of a distribution. In GEV distributions, the shape parameter affects the lower-tail and upper-tail differently. In this study, the parameters are chosen such that response distributions are positively skewed, while capacity distributions are negatively skewed to better simulate the difficulties of sampling in tails. Individually, the shape parameters $({\xi }_r,{\xi }_c)$ correspond to three types of tail heaviness: heavy, medium and light. Here, tail heaviness is considered for the upper-tail of R and the lower-tail of C. The location parameter of the capacity distribution $\left({\theta }_c\right)$ is changed while keeping the response location $\left({\theta }_r\right)$ stationary so that each combination corresponds to a failure probability of ${10}^{-4}$. Nine study cases result because of the different combination of tails possible for both R and C. The parameters for the nine study cases are presented in

Table 1. GEV distributions parameters of R and C for nine study cases.

Response Tail	${\mathit{\xi }}_{\mathit{r}},{\mathit{\theta }}_{\mathit{r}},{\mathit{\sigma }}_{\mathit{r}}$	Capacity Tail	${\mathit{\xi }}_{\mathit{c}},{\mathit{\theta }}_{\mathit{c}},{\mathit{\sigma }}_{\mathit{c}}$
		Heavy	($-1.8,1,33.5$)
Heavy	($0.2,1,0$)	Medium	($-1,1,26.5$)
		Light	($-0.52,1,26.2$)
		Heavy	($-1.8,1,30$)
Medium	($0,1,0$)	Medium	($-1,1,10.7$)
		Light	($-0.52,1,9.5$)
		Heavy	($-1.8,1,30$)
Light	($-0.12,1,0$)	Medium	($-1,1,9.4$)
		Light	($-0.52,1,6.4$)

. The distribution parameters are only utilized to generate scarce samples of R and C. It is to be noted that the GEV distributions used here are continuous over the real line $\left(\mathbb{R}\right)$ and can have negative values. Though this does not reflect the real-world scenario, it does not deter in evaluating the performance of the proposed formulation.

To capture sampling variability, multiple iterations are carried out, and the procedure employed for identifying the variability in the estimates and bootstrap bounds is presented in Algorithm 2.

Algorithm 2 Confidence bounds using bootstrap.

1: for$j=1,2,\dots ,T$do 2:     ${\mathbf{r}}_{j\times M}\sim GEV({\xi }_r,{\theta }_r,{\sigma }_r)$, ${\mathbf{c}}_{j\times N}\sim GEV({\xi }_c,{\theta }_c,{\sigma }_c)\leftarrow$ Generate samples of response and capacity. 3:     ${\widehat{p_f}}_j\leftarrow$ Estimate failure probability by applying proposed approach to samples ${\mathbf{r}}_{j\times M}$ and ${\mathbf{c}}_{j\times N}$. 4:     for $i=1,2,\dots ,B$ do 5:         ${\mathbf{rb}}_{i\times M},{\mathbf{cb}}_{i\times N}\leftarrow$ Generate bootstrap samples from original samples ${\mathbf{r}}_{j\times M}$ and ${\mathbf{c}}_{j\times N}$. 6:         ${\widehat{p_{fb}}}_i\leftarrow$ Estimate failure probability by applying proposed approach to bootstrap samples ${\mathbf{rb}}_{i\times M}$ and ${\mathbf{cb}}_{i\times N}$. 7: 25th, 50th, 75th percentiles $\leftarrow {\left\{\widehat{p_f}\right\}}_{T\times 1}$. 8: ${25\mathrm{th}}_{1\times T}$, ${50\mathrm{th}}_{1\times T}$, ${75\mathrm{th}}_{1\times T}$ percentiles $\leftarrow {\left\{\widehat{p_{fb}}\right\}}_{B\times T}$.

$${\beta }_{IS}={\Phi }^{-1}(1-\widehat{p_f})$$

(10)

The failure probability estimates calculated through Algorithm 2 for each study case are converted into reliability indices using Equation (10). To facilitate comparison, the estimates are divided by the actual reliability index $({\beta }_a=3.71)$. For all nine study cases, the sample size of response $\left(M\right)$ and capacity $\left(N\right)$ are considered to be 50. The results of applying the formulation to each study case are presented in

Table 2. Both R and C unknown case: percentiles of $MM{\beta }_{IS}/{\beta }_a$.

		Heavy C		Medium C		Light C
	Percentile	Original Sample	Bootstrap Mean (Std)	Original Sample	Bootstrap Mean (Std)	Original Sample	Bootstrap Mean (Std)
Heavy R	$25\mathrm{th}$	0.96	1.10 (0.26)	1.21	1.42 (0.29)	1.79	2.08 (0.61)
	$50\mathrm{th}$	1.06	1.17 (0.29)	1.50	1.54 (0.32)	2.13	2.46 (0.89)
	$75\mathrm{th}$	1.22	1.36 (0.31)	1.76	1.77 (0.33)	2.59	3.08 (1.31)
Medium R	$25\mathrm{th}$	0.93	1.03 (0.20)	0.85	0.95 (0.18)	1.01	1.09 (0.21)
	$50\mathrm{th}$	1.02	1.09 (0.21)	0.97	1.01 (0.19)	1.13	1.20 (0.25)
	$75\mathrm{th}$	1.18	1.25 (0.23)	1.08	1.14 (0.22)	1.30	1.40 (0.32)
Light R	$25\mathrm{th}$	0.91	1.04 (0.23)	0.78	0.85 (0.14)	0.85	0.91 (0.12)
	$50\mathrm{th}$	1.01	1.10 (0.25)	0.88	0.91 (0.14)	0.94	0.98 (0.13)
	$75\mathrm{th}$	1.18	1.28 (0.26)	0.99	1.04 (0.17)	1.03	1.09 (0.18)

. Values under the ‘Original sample’ column represent the percentiles of ratio of the estimates from the original samples (from Step 7 of Algorithm 2) repeated for T (=100) iterations. Mean and standard deviation of the percentiles (from Step 8 of Algorithm 2) from the bootstrap repetitions (B = 100) are presented under the ‘Bootstrap’ column. An accurate estimate of the reliability index would be indicated with a value of one, and high precision is indicated by low variability between the percentiles and low standard deviation in the bootstrap percentiles. It is observed that the estimation is poorer in the case of heavy-tailed response and medium-tailed capacity and heavy-tailed response and light-tailed capacity. In both these cases, the failure region is situated further into the tail of the heavy-tailed response where the scarce sample-based PDF estimation is prone to high errors. In many instances, the probability density drops to zero prematurely which results in overestimation of reliability.

In certain reliability applications, either $f_R$ or $f_C$ is known, or it is easier to obtain samples from one of the distributions of R and C. The proposed method is tested for such scenarios where only one of the distributions is unknown, and during this exercise, the same combinations of tails for R and C are assumed.

Table 3. R unknown case: failure probability estimated as per Equation (6); (PDF approximation of R): percentiles of $MM{\beta }_{IS}/{\beta }_a$.

		Heavy C		Medium C		Light C
	Percentile	Original Sample	Bootstrap Mean (Std)	Original Sample	Bootstrap Mean (Std)	Original Sample	Bootstrap Mean (Std)
Heavy R	$25\mathrm{th}$	1.02	1.02 ($9\times {10}^{-3}$)	1.67	1.71 (0.15)	3.16	3.63 (0.66)
	$50\mathrm{th}$	1.03	1.03 ($7\times {10}^{-3}$)	1.78	1.72 (0.15)	3.64	3.67 (0.67)
	$75\mathrm{th}$	1.04	1.04 ($5\times {10}^{-3}$)	1.81	1.79 (0.06)	4.00	4.09 (0.30)
Medium R	$25\mathrm{th}$	0.99	0.99 ($2\times {10}^{-3}$)	1.02	1.03 (0.06)	1.25	1.32 (0.20)
	$50\mathrm{th}$	1.00	1.00 ($2\times {10}^{-3}$)	1.05	1.04 (0.05)	1.37	1.34 (0.20)
	$75\mathrm{th}$	1.00	1.00 ($2\times {10}^{-3}$)	1.08	1.07 (0.03)	1.45	1.51 (0.14)
Light R	$25\mathrm{th}$	0.99	0.99 ($2\times {10}^{-3}$)	0.99	0.99 (0.02)	0.98	1.03 (0.10)
	$50\mathrm{th}$	1.00	1.00 ($2\times {10}^{-3}$)	1.00	1.00 (0.02)	1.07	1.06 (0.10)
	$75\mathrm{th}$	1.00	1.00 ($2\times {10}^{-3}$)	1.02	1.01 (0.02)	1.12	1.13 (0.07)

and

Table 4. C unknown case: failure probability estimated as per Equation (6); (CDF approximation of C): percentiles of $MM{\beta }_{IS}/{\beta }_a$.

		Heavy C		Medium C		Light C
	Percentile	Original Sample	Bootstrap Mean (Std)	Original Sample	Bootstrap Mean (Std)	Original Sample	Bootstrap Mean (Std)
Heavy R	$25\mathrm{th}$	0.96	0.98 (0.10)	0.99	0.99 ($6\times {10}^{-3}$)	1.00	1.00 ($2\times {10}^{-3}$)
	$50\mathrm{th}$	1.03	1.00 (0.08)	1.00	1.00 ($3\times {10}^{-3}$)	1.00	1.00 ($2\times {10}^{-3}$)
	$75\mathrm{th}$	1.06	1.04 (0.05)	1.00	1.00 ($2\times {10}^{-3}$)	1.00	1.00 ($2\times {10}^{-3}$)
Medium R	$25\mathrm{th}$	0.93	1.13 (0.31)	0.87	0.92 (0.12)	0.96	0.96 (0.05)
	$50\mathrm{th}$	1.12	1.24 (0.35)	0.96	0.96 (0.11)	0.99	0.98 (0.04)
	$75\mathrm{th}$	1.42	1.52 (0.33)	1.03	1.02 (0.07)	1.00	1.00 (0.02)
Light R	$25\mathrm{th}$	0.91	1.05 (0.23)	0.82	0.90 (0.13)	0.87	0.89 (0.09)
	$50\mathrm{th}$	1.03	1.12 (0.25)	0.93	0.95 (0.14)	0.94	0.94 (0.09)
	$75\mathrm{th}$	1.20	1.29 (0.29)	1.03	1.06 (0.14)	1.00	0.99 (0.08)

present the results for R unknown and C unknown cases, respectively, using the proposed approach. From Table 3, it is observed that both the accuracy and precision of the estimation improved for most of the tail combinations compared to both R and C unknown cases. In the case of heavy-tailed response and medium-tailed and light-tailed capacity, and in addition, medium-tailed response and light-tailed capacity, the bias in the estimate has increased. This is again due to the erroneous PDF estimation of heavier-tailed response distribution using a scarce sample. Lesser bias in both R and C unknown cases is due to the interaction between the PDF and CDF approximation.

In Table 4, larger errors are observed for heavy-tailed capacity and medium-tailed response and heavy-tailed capacity and light-tailed response. This suggests that estimating heavy-tailed distribution further into the tails results in larger errors; this is akin to the observations made for the R unknown case. Furthermore, the largest errors observed in the C unknown case are smaller compared to the largest errors from the R unknown case which suggests that TPNT is a better approximation for CDF of heavy tailed distributions compared to adaptive KDE used for PDF approximation. To demonstrate that this is indeed the case, the alternate form of the importance sampling estimator presented in Equation (7) is used for the R unknown scenario, wherein CDF of the response distribution is approximated instead of the PDF.

Table 5. R unknown case: failure probability estimated as per Equation (7); (CDF approximation of R): percentiles of $MM{\beta }_{IS}/{\beta }_a$.

		Heavy C		Medium C		Light C
	Percentile	Original Sample	Bootstrap Mean (Std)	Original Sample	Bootstrap Mean (Std)	Original Sample	Bootstrap Mean (Std)
Heavy R	$25\mathrm{th}$	0.95	0.96 (0.09)	0.95	1.14 (0.29)	0.92	1.13 (0.31)
	$50\mathrm{th}$	1.02	0.99 (0.07)	1.13	1.22 (0.31)	1.11	1.23 (0.35)
	$75\mathrm{th}$	1.03	1.03 (0.02)	1.43	1.42 (0.30)	1.38	1.51 (0.41)
Medium R	$25\mathrm{th}$	0.99	0.99 (0.02)	0.84	0.89 (0.12)	0.77	0.89 (0.18)
	$50\mathrm{th}$	0.99	0.99 (0.01)	0.94	0.94 (0.11)	0.93	0.96 (0.20)
	$75\mathrm{th}$	1.00	1.00 ($2\times {10}^{-3}$)	1.02	1.01 (0.09)	1.08	1.11 (0.21)
Light R	$25\mathrm{th}$	0.99	0.99 ($2\times {10}^{-3}$)	0.91	0.93 (0.07)	0.79	0.84 (0.12)
	$50\mathrm{th}$	1.00	1.00 ($1\times {10}^{-3}$)	0.97	0.96 (0.06)	0.87	0.90 (0.12)
	$75\mathrm{th}$	1.00	1.00 ($1\times {10}^{-3}$)	1.00	0.99 (0.04)	0.99	0.98 (0.12)

presents the results, from which it is observed that the bias has reduced compared to the PDF approximation-based estimation. However, this also results in underestimation of reliability in the case of light-tailed response and light-tailed capacity. This is a trade-off between the choice of adaptive-kernel-based PDF approximation and TPNT-based CDF approximation. From these observations, it can be surmised that it would be beneficial to know the heaviness of R and C so that appropriate estimator can be chosen based on the efficacy of the approximation methods available.

An additional observation is that for one (either R or C) unknown case, the bootstrap standard deviation is larger when the bias in the estimates from the original samples are larger and smaller for more accurate cases. Thus, in practice where there is no actual value to compare with, bootstrap standard deviation can be used to discern the accuracy of the estimate. However, this does not track well in the case of both R and C unknown.

Tail-Index Estimation

Estimation of tail index is sometimes part of the distribution identification process and in turn CDF estimation where methods such as Hill estimator [48], Pickands estimator, Generalized Pareto fits and others are applied. Hill estimator considers k upper order statistics from a sample to evaluate the tail heaviness. Equation (11) presents the estimator for a positive sample of size n with the order statistics, $X_{1,n}\le X_{2,n}\le \dots X_{n,n}$.

$$H_{k,n}=\frac{1}{k}\sum _{i=0}^{k-1}log\frac{X_{n-i,n}}{X_{n-k,n}}$$

(11)

The estimate of tail heaviness is sensitive to the number of order statistics $\left(k\right)$ used and shifting of sample. Different modifications that address such sensitivity issues exist in the literature [49,50]. As we only require comparing between tails of R and C, especially when they are very different, obtaining exact estimates is of low priority. Hence, a modified Hill estimator is used where the sample is mean-shifted, and k is considered to be $10\%$ of the sample size. For upper-tail estimation, the largest k values are used and for lower-tail estimation, the absolute values of the smallest k values are used. The Modified Hill estimator is applied to assess the tail-heaviness of upper-tail of response and lower-tail of capacity, and the sample with the heavier tail is chosen for CDF approximation, and thereby the appropriate form of the estimator is chosen for failure probability estimation.

In both unknown scenarios, the heavy-tailed R and the light-tailed C which showed the largest deviation from the actual value was chosen to test the tail index estimate-based improvement on the proposed formulation. Results presented in

Table 6. Heavy R and light C: Percentiles of $MM{\beta }_{IS}/{\beta }_t$ from tail estimation-based improved formulation.

Percentile	Original Sample	Bootstrap Mean (Std)
25th	0.93	1.19 (0.39)
50th	1.12	1.32 (0.53)
75th	1.39	1.70 (1.02)

show the use of tail index estimation improved the estimates. It is to be noted that the modified Hill estimator was successful in contrasting between heaviness of tails of R and C 81 out of 100 iterations.

5. Reliability Estimation Examples

Based on the flowchart of the proposed approach in Figure 5, benchmark reliability estimation examples are tested in this section. In all examples, to account for sample variability 100 iterations are used with samples of size 50 for R and C in each iteration. For each iteration, bootstrap samples for R and C are generated 100 times.

5.1. Example 1: Concave Limit State 1

This is a concave limit state example taken from [51].

$$G=2.62-u_2-0.15{u_1}^2$$

(12)

Both $u_1$ and $u_2$ are standard normal variables, and the concave limit state is of the separable form. Here, R is taken as $0.15{u_1}^2$, and C is taken as $2.62-u_2$.

From tail-index estimates in 88 out of 100 iterations, response R is determined to be heavier-tailed than C which is consistent with the analytical form of R and C. The response distribution resulting from squaring a standard normal variable is an ${\chi }^2$-distribution with one degree of freedom. In this case, the ${\chi }^2$-distribution has a heavier tail than the standard normal distribution.

Table 7. Concave limit state-1: percentiles of $MM{\beta }_{IS}/{\beta }_t$ for ${\beta }_t=2.39$.

Percentile	Original Sample	Bootstrap Mean (Std)
25th	0.84	0.96 (0.20)
50th	0.96	1.03 (0.20)
75th	1.12	1.20 (0.22)

Note: ${\beta }_t$ is computed through MCS with 10⁸ samples.

provides the results obtained for a target reliability ${\beta }_t=2.39$. It is observed that the true value is contained within the percentile bounds (0.84–1.12) from the original sample estimates, and the bootstrap bounds obtained for each percentile also contain the true value. The variability of the estimates from the original sample is high, and the standard deviation computed from the bootstrap estimates is reflective of the variability presented in the original sample estimates.

5.2. Example 2: Concave Limit State 2

Example 2 is a concave limit state taken from [51].

$$G={u_1}^2-5{u_2}^2+45$$

(13)

Both $u_1$ and $u_2$ are standard normal variables, and the concave limit state is of the separable form where R is taken as $5{u_2}^2$ and C is taken as ${u_1}^2+45$.

Table 8. Concave limit state-2: percentiles of $MM{\beta }_{IS}/{\beta }_t$ for ${\beta }_t=2.82$.

Percentile	Original Sample	Bootstrap Mean (Std)
25th	0.82	0.89 (0.16)
50th	0.90	0.95 (0.16)
75th	1.02	1.08 (0.19)

Note: ${\beta }_t$ is computed through MCS with 10⁸ samples.

presents the numerical results obtained for this example by applying the proposed IS approach. From tail-index estimates, in 100 out of 100 iterations, response R is determined to be heavier-tailed than C. Here, both the response and capacity distributions resulting from squaring standard normal variables $u_1$ and $u_2$ are ${\chi }^2$-distributions. However, the upper tail of response is heavier because scaling a ${\chi }^2$-distribution results in an increase of tail heaviness.

The percentiles of $MM{\beta }_{IS}/{\beta }_t$ from the original samples show that the target value is contained within the first and third quartile bounds. For all the percentiles, it is observed that the bootstrap bounds also contain the target reliability.

5.3. Example 3: Roof Truss Example

This example is discussed as a case study for reliability estimation in [45]. The schematic in Figure 6 presents a roof truss subjected to a uniformly distributed load q which is transformed into the nodal load $P=\frac{ql}{4}$. Equation (14) presents the limit state constructed for perpendicular deflection ${\Delta }_C$ at node C as:

$$G\left(x\right)=0.03-{\Delta }_C$$

(14)

where

$${\Delta }_C=\frac{q{l}^2}{2}\left(\frac{3.81}{A_c{E}_c}+\frac{1.13}{A_s{E}_s}\right)$$

(15)

Here, $A_c$ and $A_s$ are areas of cross sections of the concrete reinforced bars and steel bars, respectively. Similarly, the variables $E_c$ and $E_s$ represent the moduli of elasticity, while l denotes the length. The input variables are considered to be mutually independent random variables, distributed normally. The parameters of their distributions are as presented in

Table 9. Mean and SD of random variables for roof truss example.

Random Variable	Mean (SD)
q (N/m)	20,000 (1600)
l (m)	12 (0.24)
$A_s\left({\mathrm{m}}^2\right)$	$9.82\times {10}^{-4}(5.89\times {10}^{-5})$
$A_c\left({\mathrm{m}}^2\right)$	0.04 (0.008)
$E_s\left({\mathrm{N}/\mathrm{m}}^2\right)$	$1.2\times {10}^{11}(8.4\times {10}^9)$
$E_c\left({\mathrm{N}/\mathrm{m}}^2\right)$	$3\times {10}^{10}(2.4\times {10}^9)$

.

The limit state is transformed into $G(C,R)=C-R$ form where $C=\frac{0.03}{q}$ and $R=\frac{{\Delta }_C}{q}$ becomes a function of five random variables that correspond to the geometry and material properties of the truss members.

Table 10. Roof truss example: percentiles of $MM{\beta }_{IS}/{\beta }_t$ for ${\beta }_t=3.4$.

Percentile	Original Sample	Bootstrap Mean (Std)
25th	0.84	0.94 (0.17)
50th	0.94	1.02 (0.20)
75th	1.07	1.18 (0.29)

Note: ${\beta }_t$ is computed through MCS with 10⁸ samples.

provides the estimates of reliability and the corresponding bootstrap bounds for a target reliability of ${\beta }_t=3.4$. It is observed that at each of the percentiles, the bootstrap bound captures the target reliability. From tail-index estimates, in 77 out of 100 iterations, response R is determined to be heavier-tailed than C.

5.4. Example 4: Propped Cantilever Beam Example

Figure 7 presents the schematic of the example which is taken from [52]. Equation (16) presents the original limit state for the maximum deflection of the beam ${\nu }_{max}$ against the maximum allowable deflection ${\nu }_{crit}$:

$$G={\nu }_{max}-{\nu }_{crit}$$

(16)

where deflection of the beam $\nu \left(x\right)$ is measured as per Equation (17), and for the considered loading condition, the maximum deflection $\left({\nu }_{max}\right)$ is obtained at $x=0.5528L$.

$$\nu \left(x\right)=\frac{q_0{x}^2}{120LEI}(4L^3-8L^{2}x-5L{x}^2-x^3)$$

(17)

where

$$I=\frac{b_f{d}^3-(b_f-t_w){(d-2t_f)}^3}{12}$$

(18)

The limit state is modified to convert into $G(C,R)=C-R$ form, where $C({\nu }_{crit},E)={\nu }_{crit}E$, and R is a function of remaining seven random variables. Mean and standard deviation (SD) of the normally distributed random variables is presented in

Table 11. Mean and SD of random variables for propped cantilever beam example.

Random Variable	Mean (SD)
$q_0$ (kN/m)	20 (2)
L (m)	6 (0.3)
E (GPa)	210 (10)
d (cm)	25 (0.5)
$b_f$ (cm)	25 (0.5)
$t_w$ (cm)	2 (0.2)
$t_f$ (cm)	2 (0.2)

Note: All random variables follow normal distribution.

. The critical displacement ${\nu }_{crit}$ is changed between 4–5 mm to generate three target reliability situations.

Table 12. Propped cantilever beam: percentiles of $MM{\beta }_{IS}/{\beta }_t$ for different ${\nu }_{crit}$ and ${\beta }_t$.

	${\mathit{\nu }}_{\mathit{crit}}=4.0\mathbf{mm},{\mathit{\beta }}_{\mathit{t}}=2.98$		${\mathit{\nu }}_{\mathit{crit}}=4.5\mathbf{mm},{\mathit{\beta }}_{\mathit{t}}=3.50$		${\mathit{\nu }}_{\mathit{crit}}=5.0\mathbf{mm},{\mathit{\beta }}_{\mathit{t}}=3.97$
Percentile	Original	Bootstrap Mean (Std)	Original	Bootstrap Mean (Std)	Original	Bootstrap Mean (Std)
25th	0.84	1.03 (0.30)	0.82	0.99 (0.27)	0.81	0.99 (0.30)
50th	0.99	1.13 (0.36)	0.99	1.09 (0.32)	0.96	1.12 (0.39)
75th	1.29	1.31 (0.42)	1.24	1.29 (0.44)	1.26	1.35 (0.52)

Note: ${\beta }_t$ values are computed through MCS with 10⁸ samples.

presents the results obtained from the proposed IS formulation for different target reliability indices $\left({\beta }_t\right)$. As the target reliability increases, the variability from both the original sample estimates and bootstrap sample estimates remain similar.

6. Application to RBDO Examples

The proposed importance sampling approach is demonstrated on two benchmark and two real-world RBDO examples. Algorithm 3 delineates the steps followed to perform RBDO using the proposed approach of reliability estimation. For the benchmark examples, efficiency between the proposed approach and a crude Monte Carlo approach is compared using the total number of limit state evaluations which is computed as $N_{doe}\times N_{is}$ for the proposed approach and $N_{doe}\times N_{mcs}$ for the MCS-based approach. For the real-world examples, MCS is used only to validate the optima obtained through the proposed approach.

Algorithm 3 RBDO using proposed importance sampling approach.

1: $\mathbf{d}=\{d_1,d_2,\dots ,d_{N_{doe}}\}\leftarrow$ Generate a design of experiment (DoE) in design variable space. 2: for  $i=1,2,\dots ,N_{doe}$do 3:     Propagate uncertainty at DoE points. 4:     $\mathbf{r},\mathbf{c}\leftarrow$ Obtain response and capacity samples through simulations or physical experiments. In the examples, analytical functions are utilized. 5:     ${\widehat{p_f}}_i\leftarrow$ Apply proposed importance sampling formulation to samples $\mathbf{r},\mathbf{c}$. 6:     ${\beta }_{I{S}_i}\leftarrow$ Transform failure probability to reliability index using (10). 7: $\widehat{{\beta }_{\mathbf{IS}}}\leftarrow$ Construct surrogate model for reliability index using $\left\{{\beta }_{I{S}_i}\right\}$ and $\mathbf{d}$. 8: $d^{*}\leftarrow$ Optimization using $\widehat{{\beta }_{\mathbf{IS}}}$ for evaluation of reliability constraint.

6.1. Cantilever Beam Example

This is a standard RBDO example taken from [53]. In this example, the weight of a cantilever beam (shown in Figure 8) is minimized while considering two failure modes: bending stress does not exceed yield strength (19) and the tip displacement does not exceed the allowable displacement limit of $2.5in$ (20).

$$G_s={\sigma }_y-\frac{6L}{wt}\left(\frac{X}{w}+\frac{Y}{t}\right)$$

(19)

$$G_d=D_0-\frac{4L^3}{Ewt}\sqrt{{\left(\frac{Y}{t^2}\right)}^2+{\left(\frac{X}{w^2}\right)}^2}$$

(20)

The length of the beam $\left(L\right)$ and the density are held constant, while the width $\left(w\right)$ and thickness $\left(t\right)$ of the beam are allowed to change which transforms the objective function from weight to area of cross section, $A=wt$. The horizontal $\left(X\right)$ and vertical $\left(Y\right)$ loads along with modulus of elasticity $\left(E\right)$ are random variables whose uncertainty characteristics are presented in

Table 13. Random variables for cantilever beam example.

Random Variable	Mean (SD)
$X\left(lb\right)$	500 (100)
$Y\left(lb\right)$	${10}^3$ (100)
${\sigma }_y\left(psi\right)$	$4\times {10}^4$ ($2\times {10}^3$)
$E\left(psi\right)$	$2.9\times {10}^7$ ($1.45\times {10}^6$)

Note: All random variables follow normal distribution.

.

RBDO of the cantilever beam requires that the failure probability of both the limits states does not exceed $1.35\times {10}^{-3}$ which translates to a target reliability index ${\beta }_t\ge 3$. Thus, the optimization is formulated as:

$$\begin{array}{cc}\hfill \underset{w,t}{min}& A=wt\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& g1={\Phi }^{-1}(1-Pr(G_s\le 0))\ge 3\hfill \\ \hfill & g2={\Phi }^{-1}(1-Pr(G_d\le 0))\ge 3\hfill \end{array}$$

(21)

As stated in Algorithm 3, a DoE of size $(N_{doe}=40)$ is created in $d=(w,t)$ space using latin hypercube sampling (LHS). Additionally, four corner points are added to the DoE. At each design point $\left(d_i\right)$, samples of response and capacity for both limit states are obtained by simulating uncertainty as per Table 13, considering a sample size of $N_{is}=50$. For limit state $G_s$, ${\sigma }_y$ is considered as capacity, while the rest of the equation is considered as a response. Similarly, for limit state $G_d$, $D_0\times E$ is taken as capacity by regrouping the random variables in the limit state function to convert it into a separable form [54]. Using these samples, failure probability estimates for both limit states are obtained through the proposed approach. A surrogate model for reliability index $\beta$ is constructed, and this surrogate is used for constraint evaluation during optimization.

In a similar fashion, MCS-based failure probability estimates are also obtained at the design points using samples of size $N_{mcs}={10}^6$. A surrogate model using the MCS estimates is constructed for a reliability index, and optimization is carried out. At design points that are closer to the lower and upper bounds, failure probability estimates could be either zeros or ones. These estimates when converted into reliability indices become $-\infty$ and ∞; hence, these are modified during surrogate construction. These singularities are more common in the case of MCS estimation compared to IS estimation.

The optima obtained from both MCS and IS are compared in

Table 14. RBDO results of cantilever beam example.

Surrogate Model for $\mathit{\beta }$ Constraints	Reliability Estimation	Optima $\left({\mathit{d}}^{*}\right)$		Objective Function Value	$\widehat{\mathit{\beta }}$ at ${\mathit{d}}^{*}$		${\mathit{\beta }}_{\mathit{MCS}}$ at ${\mathit{d}}^{*}$
		$\mathit{w}$ (in)	$\mathit{t}$ (in)	$\mathit{A}$ (in${}^{\mathbf{2}}$)	${\mathit{g}}_{\mathbf{1}}$	${\mathit{g}}_{\mathbf{2}}$	${\mathit{g}}_{\mathbf{1}}$	${\mathit{g}}_{\mathbf{2}}$
WAS	IS	2.59	3.74	9.69	3.00	3.64	3.25	3.69
	MCS	2.59	3.66	9.50	3.00	3.44	2.95	3.39

. As the cost of the knowledge of random variables is not quantified, computational cost is only compared with MCS. Additionally, different surrogate choices, such as PRS, RBF, Kriging and weighted average surrogate (WAS), were considered for $\widehat{\mathbf{fi}}$ construction, and the optima obtained by using WAS model are presented. The accuracy of the surrogate for both constraints is evaluated using a generalized mean square error (GMSE) metric. GMSE for constraint $g_1$ using IS-based estimates is reported as 0.91 (mean of reliability estimates = 2.21), whereas the error from the MCS-based surrogate is 0.09 (mean of reliability estimates = 1.97). For constraint $g_2$, GMSE (vs. mean) is reported as 0.70 (vs. 1.98) for the IS surrogate and 0.12 (vs. 1.41) for the MCS surrogate.

To validate the results, at the optima, a reliability index is calculated using MCS with sample size of ${10}^7$ which is also presented in the Table 14. It is observed that the surrogate from IS estimates has less global accuracy. However, reliability indices obtained at the optima using MCS $({\beta }_{MCS}$ at $d^{*})$ suggest that the surrogate is reasonably approximated at the constraint boundaries. The optima obtained from IS and MCS are very close, but the computational savings are hugely in favour of IS $\left(50\mathrm{vs}.{10}^6\right)$.

6.2. Bracket Structure Example

This is originally from [55] where a bracket structure is subjected to a tip load $\left(P\right)$ in addition to its own weight due to gravity (g) as presented in Figure 9. The weight of the bracket structure is optimized while considering two failure modes:

(i)

Maximum bending stress of beam CD at point B $\left({\sigma }_B\right)$ does not exceed its yield strength $\left(f_y\right)$,

(ii)

Maximum axial load on beam AB $\left(F_{AB}\right)$ does not exceed the Euler critical buckling load $\left(F_{buckling}\right)$.

Equations (22) and (23) present the limit states as:

$$\begin{array}{cc}\hfill & G_{CD}=f_y-{\sigma }_B,\mathrm{where}\hfill \end{array}$$

(22a)

$$\begin{array}{cc}\hfill & {\sigma }_B=\frac{6M_B}{w_{CD}{t}^2}\hfill \end{array}$$

(22b)

$$\begin{array}{cc}\hfill & M_B=\frac{PL}{3}+\frac{\rho g{w}_{CD}t{L}^2}{18}\hfill \end{array}$$

(22c)

$$\begin{array}{cc}\hfill & G_{AB}=F_{buckling}-F_{AB},\mathrm{where}\hfill \end{array}$$

(23a)

$$\begin{array}{cc}\hfill & F_{buckling}=\frac{{\pi }^{2}Et{w}_{AB}^{3}9sin{\theta }^2}{48L^2}\hfill \end{array}$$

(23b)

$$\begin{array}{cc}\hfill & F_{AB}=\frac{1}{cos\theta }\left(\frac{3P}{2}+\frac{3\rho g{w}_{CD}tL}{4}\right)\hfill \end{array}$$

(23c)

During RBDO of the bracket structure, the target reliability index for both limit states is $\beta \ge 2$, and the design parameters are means of the geometrical parameters of the structure: width of CD $\left({\mu }_{w_{CD}}\right)$, width of AB $\left({\mu }_{w_{AB}}\right)$ and thickness of AB and CD $\left({\mu }_t\right)$ which are bounded between $50\mathrm{mm}$ and $300\mathrm{mm}$. The uncertainty characteristics of the random variables is presented in

Table 15. Random variables for bracket structure example.

Type	Variable	Distribution	Mean	C.o.V
Random	P (kN)	Gumbel	100	15%
	E (GPa)	Gumbel	200	8%
	$f_y$ (MPa)	Lognormal	225	8%
	$\rho (\mathrm{kg}·{\mathrm{m}}^{-3})$	Weibull	7860	10%
	L (m)	Gaussian	5	5%
Design	$w_{AB}$ (mm)	Gaussian	${\mu }_{w_{AB}}$	5%
	$w_{CD}$ (mm)	Gaussian	${\mu }_{w_{CD}}$	5%
	t (mm)	Gaussian	${\mu }_t$	5%

. Thus, the RBDO is formulated as:

$$\begin{array}{cc}\hfill \underset{w_{CD},w_{AB},t}{min}& C={\mu }_{\rho }{\mu }_t{\mu }_L\left(\frac{4\sqrt{3}}{9}{\mu }_{w_{AB}}+{\mu }_{w_{CD}}\right)\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& g1={\Phi }^{-1}(1-Pr(G_{CD}\le 0))\ge 2\hfill \\ \hfill & g2={\Phi }^{-1}(1-Pr(G_{AB}\le 0))\ge 2\hfill \\ \hfill & 50\le {\mu }_{w_{CD}},{\mu }_{w_{AB}},{\mu }_t\le 300\left(\mathrm{in}\mathrm{mm}\right)\hfill \end{array}$$

(24)

The steps enumerated in Algorithm 3 are followed using a DoE of size $(N_{doe}=60)$, and scarce sample sizes of R and C during the IS approach are considered as $N_{is}=75$, while $N_{mcs}={10}^6$ realizations are used in MCS. The RBDO results are presented in

Table 16. RBDO results of bracket structure example.

Surrogate Model for $\mathit{\beta }$ Constraints	Reliability Estimation	Optima $\left({\mathit{d}}^{*}\right)$			Objective Function Value	$\widehat{\mathit{\beta }}$ at ${\mathit{d}}^{*}$		${\mathit{\beta }}_{\mathit{MCS}}$ at ${\mathit{d}}^{*}$
		${\mathit{w}}_{\mathit{AB}}\left(\mathrm{mm}\right)$	${\mathit{w}}_{\mathit{CD}}\left(\mathrm{mm}\right)$	$\mathit{t}\left(\mathrm{mm}\right)$	$\mathit{Weight}$ (kg)	${\mathit{g}}_{\mathbf{1}}$	${\mathit{g}}_{\mathbf{2}}$	${\mathit{g}}_{\mathbf{1}}$	${\mathit{g}}_{\mathbf{2}}$
WAS	IS	58	89	300	1576	2.00	2.00	2.59	2.87
	MCS	62	77	300	1474	2.00	2.00	2.02	3.55

which has the same format as the first example. GMSE (vs. mean) for constraint $g_1$ using IS-based estimates is reported as 0.58 (1.37), whereas the error from the MCS-based surrogate is 0.24 (1.02). For constraint $g_2$, GMSE (vs. mean) is reported as 0.97 (vs. 3.72) for the IS surrogate and 0.30 (vs. 3.08) for the MCS surrogate. It is observed that the IS estimate-based surrogate has less global accuracy; however, the MCS-based reliability indices computed at the optima show that it approximates reasonably well near the constraint boundaries. It is to be noted that in both engineering examples, we assume that both R and C are unknown. Any knowledge of the uncertainty characteristics of either R or C will only improve the accuracy of the reliability estimates.

6.3. Torque Arm Example

This example presents RBDO of torque arm where the mass of the component is to be minimized adhering to a probabilistic constraint on the allowable stress. Unlike previous examples, there is no analytical expression available for limit state evaluation in this case. It is a shape optimization problem originally from Bennett and Botkin [56]. Researchers have used it as a benchmark example for reliability estimation [52,57]. Rahman and Wei [58] perform RBDO where a constant allowable stress limit is considered.

In the current study, seven design variables (see Figure 10) as per [59] are considered for altering the shape of the torque arm. Figure 11 presents the base design of the torque arm around which the optimization is to be performed. Here, a horizontal load $(F_x=-2789$ N) and a vertical load $(F_y=5066$ N) are applied at the right hole while the left hole is fixed. The torque arm has modulus of elasticity $E=207\times {10}^5 \mathrm{N}·{\mathrm{cm}}^{-2}$, density $\rho =7.850\times {10}^{-3} \mathrm{kg}·{\mathrm{cm}}^{-3}$ and Poisson’s ratio $\nu =0.3$.

Equation (25) presents the limit state equation for the torque arm.

$$G={\sigma }_{max}(d_i,F_x,F_y)-{\sigma }_{all}$$

(25)

where ${\sigma }_{max}$ is the maximum von Mises stress developed in the torque arm, and ${\sigma }_{all}$ is the allowable stress limit. Since there is no analytical expression that relates the design variables with the response ${\sigma }_{max}$, finite-element analysis is used to compute the stresses developed in the torque arm for a given loading condition. A MATLAB finite-element toolbox developed by CALFEM [60] is used for this purpose. The thickness of the finite-elements in the mesh is considered to be 0.3 cm. At the end of each finite-element analysis, the maximum of the stresses is used as ${\sigma }_{max}$ in Equation (25).

For RBDO, the design variables $(d_1$ to $d_7)$, loads $(F_x$ and $F_y)$ and the allowable stress $\left({\sigma }_{all}\right)$ are considered to be uncertain. Each design point $\left(d_i\right)$ is considered to be normally distributed about itself with a coefficient of variation of 10%. The uncertainty characteristics of the remaining random variables is presented in

Table 17. Mean and SD of random variables in torque arm example.

Random Variable	Distribution Type	Mean; SD
$F_x$ (N)	Normal	$-2789;278.9$
$F_y$ (N)	Normal	$5066;506.6$
${\sigma }_{all}$ (MPa)	Lognormal	$800;80$

.

Equation (26) presents the RBDO formulation of the torque arm where a target reliability index of three is considered.

$$\begin{array}{cc}\hfill \underset{\{d_1,\dots ,d_7\}}{min}& Mass\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& g={\mathsf{\Phi }}^{-1}(1-Pr(G\le 0))\ge 3\hfill \end{array}$$

(26)

As per Algorithm 3, a DoE of size $(N_{doe}=200)$ is generated where reliability is estimated using an IS approach. Sample sizes of $R={\sigma }_{max}$ and $C={\sigma }_{all}$ during IS approach are considered as $N_{is}=100$. An RBF-based surrogate is constructed using IS-based reliability estimates. The error (GMSE vs. range) for the RBF surrogate of $\beta$ is found to be 6% which indicates a good fit. This surrogate is used for constraint (g) evaluation during optimization.

The stress contour of the torque arm design obtained as a result of the optimization is presented in Figure 12. Here, mean values for design parameters and mean loading condition are considered. The maximum von Mises stress is observed to be 523.92 MPa for the optimum design. The optimal mass of the torque arm is 0.801 kg. It can be observed that the mass is distributed to meet the target reliability. An MCS (with ${10}^5$ sample size) is used to validate the reliability of RBDO optima obtained from the IS approach, and it is observed that ${\beta }_{MCS}=3.00$. For the Monte Carlo simulation of ${10}^5$ simulations, it took 8.7 h using parallel computing toolbox of MATLAB on a system with the following specifications: Intel Xeon 10 Core 2.20 GHz 64 bit processor, with 32GB RAM.

Table 18. Design variable bounds and optimum design of torque arm (in cm).

DV	${\mathit{d}}_{\mathit{L}}$	${\mathit{d}}_{\mathit{U}}$	${\mathit{d}}^{*}$ (Optimum)
$d_1$	1.80	3.20	2.15
$d_2$	1.25	1.60	1.28
$d_3$	1.20	4.60	1.59
$d_4$	−0.10	0.40	−0.09
$d_5$	−0.30	0.30	0.30
$d_6$	−0.90	0.80	0.30
$d_7$	0.40	1.80	0.54

presents the design variable bounds and the optimum design obtained from RBDO.

Though the errors in the individual reliability estimations are quantified via bootstrap, they were not propagated into the surrogate model during RBDO. Using the surrogate model for reliability index evaluation instead of direct evaluation enabled the smoothing of noise from the IS-based estimation which leads to better convergence during optimization. From MCS validation at the optima, it is observed that the proposed approach results in heavier but safer designs. Further analysis is required to understand the error propagation from different stages of the proposed approach.

6.4. Car Side-Impact Problem—A Multi-Objective Reliability-Based Design Optimization (MORBDO) Example

We demonstrate the proposed methodology on an MORBDO example taken from [61]. In this example, the objective is to minimize the weight $\left(f_1\right)$ of a car as well as the average rib deflection $\left(f_2\right)$ during a crash. A car is subjected to a side-impact based on European Enhanced Vehicle-Safety Committee (EEVC) procedures. The effect of the side-impact on a dummy in terms of head injury criteria, load in the abdomen, pubic symphysis force, viscous criterion, and rib deflections at the upper, middle and lower rib locations are considered as constraints. The MORBDO formulation is made up of seven uncertain design variables ($x_1,\dots ,x_7$) and four random variables ($p_1,\dots ,p_4$). Equation (27) presents the optimization formulation of the car side-impact problem:

$$\begin{array}{cc}\hfill \underset{{\mu }_{\mathbf{x}}}{min}& f_1=f({\mu }_{\mathbf{x}},{\mu }_{\mathbf{p}})\hfill \\ \hfill \underset{{\mu }_{\mathbf{x}}}{min}& f_2={\displaystyle \frac{g_2({\mu }_{\mathbf{x}},{\mu }_{\mathbf{p}})+g_3({\mu }_{\mathbf{x}},{\mu }_{\mathbf{p}})+g_4({\mu }_{\mathbf{x}},{\mu }_{\mathbf{p}})}{3}}\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& {\mathsf{\Phi }}^{-1}\left(1-\mathbb{P}\mathbf{r}\left(g_i\left(\mathbf{x},\mathbf{p}\right)\right)\le b_i\right)\ge {\beta }_t,\hfill \\ \hfill & i=1,\dots ,10.\hfill \end{array}$$

(27)

where

$$\begin{array}{cc}\hfill g_1\left(\mathbf{x},\mathbf{p}\right)\equiv & \mathrm{Abdomen}\mathrm{load}\le 1\mathrm{kN}\hfill \\ \hfill g_2\left(\mathbf{x},\mathbf{p}\right)\equiv & \mathrm{Upper}\mathrm{rib}\mathrm{deflection}\le 32\mathrm{mm}\hfill \\ \hfill g_3\left(\mathbf{x},\mathbf{p}\right)\equiv & \mathrm{Middle}\mathrm{rib}\mathrm{deflection}\le 32\mathrm{mm}\hfill \\ \hfill g_4\left(\mathbf{x},\mathbf{p}\right)\equiv & \mathrm{Lower}\mathrm{rib}\mathrm{deflection}\le 32\mathrm{mm}\hfill \\ \hfill g_5\left(\mathbf{x},\mathbf{p}\right)\equiv & \mathrm{Upper}\mathrm{viscous}\mathrm{criteria}\le 0.32\mathrm{m}/\mathrm{s}\hfill \\ \hfill g_6\left(\mathbf{x},\mathbf{p}\right)\equiv & \mathrm{Middle}\mathrm{viscous}\mathrm{criteria}\le 0.32\mathrm{m}/\mathrm{s}\hfill \\ \hfill g_7\left(\mathbf{x},\mathbf{p}\right)\equiv & \mathrm{Lower}\mathrm{viscous}\mathrm{criteria}\le 0.32\mathrm{m}/\mathrm{s}\hfill \\ \hfill g_8\left(\mathbf{x},\mathbf{p}\right)\equiv & \mathrm{Pubic}\mathrm{symphysis}\mathrm{force}\le 4.0\mathrm{kN}\hfill \\ \hfill g_9\left(\mathbf{x},\mathbf{p}\right)\equiv & \mathrm{Velocity}\mathrm{of}\mathrm{B}-\mathrm{pillar}\mathrm{at}\mathrm{middle}\mathrm{point}\hfill \\ \hfill & \le 10\mathrm{mm}/\mathrm{ms}\hfill \\ \hfill g_{10}\left(\mathbf{x},\mathbf{p}\right)\equiv & \mathrm{Velocity}\mathrm{of}\mathrm{front}\mathrm{door}\mathrm{at}\mathrm{B}-\mathrm{pillar}\hfill \\ \hfill & \le 15.7\mathrm{mm}/\mathrm{ms};\hfill \end{array}$$

$$\begin{array}{cc}\hfill 1.0\le x_1& \le 1.5,0.45\le x_2\le 1.0,0.5\le x_3\le 1.5,\hfill \\ \hfill 0.5\le x_4& \le 1.5,0.875\le x_5\le 2.625,0.4\le x_6\le 1.2,\hfill \\ \hfill 0.4\le x_7& \le 1.2,{\mu }_{p_1}=0.345,{\mu }_{p_2}=0.192,{\mu }_{p_3}={\mu }_{p_4}=0.\hfill \end{array}$$

Analytical expressions of the objective functions and constraints as well as the physical descriptions of the design variables ($x_1-x_7$), and random variables ($p_1-p_4$) are presented in Appendix C.

In order to solve the MORBDO problem as per Algorithm 3, a DoE of size $(N_{doe}=200)$ is generated within the design bounds where reliability of the ten constraints ($g_1,\dots ,g_{10}$) is estimated using the IS approach. In this example, the number of design variables is seven; hence, $2^7=128$ corner points have been sampled. Next, we performed a space-filling sampling using LHS. In each dimension, we sampled 10 points (in total $7\times 10=70$) using LHS design as per a thumb rule in DoE sampling [62]. Without loss of generalization, we added two to make a round sampling number. Sample sizes of response and capacity during IS are considered as $N_{is}=50$. PRS surrogate for reliability indices of the ten constraints is constructed using the IS-based estimates. In order to improve the surrogate accuracy, an additional DoE of size 200 is generated within adjusted design bounds. The final surrogate errors (GMSE vs. range) for all constraints is found to range between $3.76\%$ and $8.38\%$. Using these surrogates for constraint evaluation, the MOO problem presented in Equation (27) is solved using NSGA-II [63].

The Pareto optimal front corresponding to the different target reliability indices (${\beta }^t=[2.5,3.0,3.5]$) along with the deterministic Pareto optimal front is presented in Figure 13. The NSGA-II algorithm is applied for the four instances using the population size of 200 for 100 generations. The GA uses the following parameters: number of offspring is 50, probability ($p_c$) of simulated binary crossover (SBX) is 0.9, crossover parameter (${\eta }_c$) is 20, probability of mutation ($p_m$) is 0.9, and the mutation parameter (${\eta }_m$) is 50.

It is evident from Figure 13 that as the targeted reliability increases, the respective reliable Pareto optimal front shifts inside the feasible criterion space and away from the deterministic Pareto optimal front to ensure more reliable solutions. The Pareto solutions were further validated using MCS ($N_{mcs}={10}^6$), and all the solutions were found to meet the target reliability.

7. Conclusions

A scarce sample-based importance sampling approach to estimate reliability is proposed when there is little or no information about the uncertainty characteristics of the random variables involved. The proposed formulation was tested on different tail heaviness for R and C distributions. In the case of one of the distributions being heavy-tailed, a tail-index estimate-based improvement to choose between PDF and CDF approximation was employed and shown to improve the accuracy (50th percentile by 1.9 times). Confidence bounds on the reliability estimate obtained through the bootstrap procedure have been shown to be indicative of the accuracy of the estimate. The proposed IS approach has been applied for reliability estimation and RBDO examples and found to be effective in terms of computational savings ($50\times 44=2200$ for cantilever beam and $75\times 60=4500$ for bracket structure example) as compared to MCS where the sample size for each reliability estimate was ${10}^6$. The approach is demonstrated on a non-analytical RBDO example which yielded a design that met the target reliability index (validated by MCS). The proposed approach has also been demonstrated on the car side-impact problem which is a multi-objective reliability-based design optimization example.

While the tail-index estimate-based alternative reduces the errors from the approximation, establishing the superiority of CDF approximation versus PDF approximation was only achieved through post-analysis of results for the specific methods used in this work. Future work could include incorporating region-wise best methods for both PDF and CDF approximation, using a suite of methods assessed through cross-validation errors. Active learning approaches could be employed during optimization to further reduce the number of design points, thereby reducing the number of reliability estimations required to obtain the optima.