Simulation methods#

Simulation estimates the event probability directly from sampled evaluations, including when a local design-point approximation is inadequate. The starting point is the probability integral in Equation (2). Direct Monte Carlo is the basic estimator; importance sampling, line sampling and subset simulation seek to reduce the cost of rare-event estimation. [Nowak2000] [Faber2009]

See Estimating failure probability by simulation for the practical workflow and FORM, SORM and design-point diagnostics for the local approximations being compared.

The principle of simulation methods is to carry out random sampling in the physical (or standardized) space. For each of the samples the limit state function is evaluated to figure out, whether the configuration is desired or undesired. The probability of failure \(p_f\) is estimated by the number of undesired configurations, respected to the total numbers of samples. [Lemaire2010]

For this analysis Equation (2) can be rewritten as

(1)#\[ p_f = P(g(\mathbf{X})\leq 0) = \iiint\limits_{g(\mathbf{X})\leq 0} I(g(\mathbf{X})\leq 0) f_{\mathbf{X}}(\mathbf{x}) d \mathbf{x}\]

where \(I\) is an indicator function that is equals to 1 if \(g(\mathbf{X}) \leq 0\) and otherwise 0. Equation (1) can be interpreted as expected value of the indicator function. Therefore, the probability of failure can be estimated such as [Malioka2009]

(2)#\[ \tilde{p}_f = \text{Ex}\left[I(g(\mathbf{X})\leq 0)\right] = \frac{1}{n}\sum_{i=1}^{n} I(g(\mathbf{X})\leq 0)\]

Crude Monte Carlo simulation#

The Crude Monte Carlo simulation (CMC) is the most simple form and corresponds to a direct application of Equation (2). A large number \(n\) of samples are simulated for the set of random variables \(\mathbf{X}\). All samples that lead to a failure are counted \(n_f\) and after all simulations the probability of failure \(p_f\) may be estimated by [Faber2009]

(3)#\[ \tilde{p}_f = \frac{n_f}{n}\]

Theoretically, an infinite number of simulations will provide an exact probability of failure. However, time and the power of computers are limited; therefore, a suitable amount of simulations \(n\) are required to achieve an acceptable level of accuracy. One possibility to reach such a level is to limit the coefficient of variation CoV for the probability of failure. [Lemaire2010]

(4)#\[ \text{CoV} = \sqrt{\frac{1-p_f}{n p_f}} \approx \frac{1}{\sqrt{n p_f}} \qquad \text{for} \quad p_f \to 0\]

Importance sampling#

To decrease the number of simulations and the coefficient of variation, other methods can be performed. One commonly applied method is the Importance Sampling simulation method (IS). Here the prior information about the failure surface is added to Equation (1)

(5)#\[ p_f = P(g(\mathbf{X})\leq 0) = \iiint\limits_{g(\mathbf{X})\leq 0} I(g(\mathbf{X})\leq 0) \frac{f_{\mathbf{X}}(\mathbf{x})}{h_{\mathbf{X}}(\mathbf{x})} h_{\mathbf{X}}(\mathbf{x}) d \mathbf{x}\]

where \(h_{X} (\mathbf{X})\) is the importance sampling probability density function of \(\mathbf{X}\). Consequently Equation (2) is extended to [Faber2009]

(6)#\[ \tilde{p}_f = \text{Ex}\left[I(g(\mathbf{X})\leq 0) \frac{f_{\mathbf{X}}(\mathbf{x})}{h_{\mathbf{X}}(\mathbf{x})}\right] = \frac{1}{n}\sum_{i=1}^{n} I(g(\mathbf{X})\leq 0)\frac{f_{\mathbf{X}}(\mathbf{x})}{h_{\mathbf{X}}(\mathbf{x})}\]

The key to this approach is to choose \(h_{X} (\mathbf{X})\) so that samples are obtained more frequently from the failure domain. For this reason, often a FORM (or SORM) analysis is performed to find a prior design point. [Baker2010]

Joint-density contours and a shaded failure region in physical coordinates.

Representation of a physical space with a set \(\mathbf{X}\) of any two random variables. The shaded area denotes the failure domain and g(mathbf{X}) = 0 the failure surface.

Direct Monte Carlo samples, with failures beyond the limit-state boundary.

For the CMC method every dot corresponds to one configuration of the random variables \(\mathbf{X}\). Dots in shaded areas lead to a failure.

Importance samples concentrated around the design point on the failure boundary.

The IS simulation method uses a distribution centered on the design point \(\mathbf{x}^*\), is obtained from a FORM (or SORM) analysis. More dots in the failure domain can be observed.

Line sampling#

Line Sampling (LS) is a variance-reduction technique that exploits the important direction \(\boldsymbol{\alpha}\) identified by FORM to reduce the n-dimensional sampling problem to a family of one-dimensional problems [Koutsourelakis2004].

The important direction \(\boldsymbol{\alpha}\) is the unit vector from the origin in standard normal space toward the most probable failure point. For each of \(N\) random samples \(\mathbf{u}_i\) drawn from \(\mathcal{N}(\mathbf{0}, \mathbf{I})\), the component along \(\boldsymbol{\alpha}\) is projected out to obtain the foot-point

\[\mathbf{v}_i = \mathbf{u}_i - \left(\mathbf{u}_i^T \boldsymbol{\alpha}\right) \boldsymbol{\alpha}\]

which lies in the \((n-1)\)-dimensional hyperplane perpendicular to \(\boldsymbol{\alpha}\). A root-finding step then locates the scalar \(c_i\) such that

\[g\!\left(\mathbf{v}_i + c_i\,\boldsymbol{\alpha}\right) = 0\]

The failure probability is estimated as the average of the one-dimensional conditional failure probabilities along each line:

(7)#\[\hat{p}_f = \frac{1}{N} \sum_{i=1}^{N} \Phi(-c_i)\]

where \(\Phi\) is the standard normal CDF. Each term \(\Phi(-c_i)\) is the probability that a point drawn from \(\mathcal{N}(0,1)\) along the \(i\)-th line lies in the failure domain.

The variance of the estimator is

\[\widehat{\operatorname{Var}}\!\left[\hat{p}_f\right] = \frac{1}{N}\,\operatorname{Var}\!\left[\Phi(-c_i)\right]\]

giving a coefficient of variation

\[\text{CoV} = \frac{\operatorname{Std}\!\left[\Phi(-c_i)\right]}{\sqrt{N}\,\hat{p}_f}\]

Line Sampling is particularly efficient when the failure surface is nearly planar near the design point, because all \(c_i\) are then close to \(\beta_{\text{FORM}}\) and \(\operatorname{Var}[\Phi(-c_i)]\) is small.

Subset simulation#

Subset Simulation (SS) is an adaptive simulation method that decomposes the rare failure event \(F = \{g(\mathbf{u}) \le 0\}\) into a sequence of more frequent nested intermediate events [AuBeck2001]:

\[F_1 \supset F_2 \supset \cdots \supset F_m = F\]

where \(F_j = \{g(\mathbf{u}) \le y_j\}\) and the thresholds satisfy \(y_1 > y_2 > \cdots > y_m = 0\). By the chain rule of probability,

(8)#\[p_f = P(F_1) \prod_{j=2}^{m} P(F_j \mid F_{j-1})\]

Each conditional probability is targeted at a user-specified level \(p_0\) (typically 0.1), making every factor in the product relatively large and easy to estimate.

Algorithm

  1. Level 0 — Generate \(N\) samples from \(\mathcal{N}(\mathbf{0}, \mathbf{I})\) and evaluate the LSF. Choose \(y_1\) as the \(p_0\)-th quantile of the LSF values, so that \(N p_0\) samples satisfy \(g \le y_1\). If \(y_1 \le 0\), the failure probability is estimated directly as \(\hat{p}_f = N_{\text{fail}} / N\).

  2. Levels \(j \ge 1\) — Use the \(N p_0\) samples satisfying \(g \le y_{j-1}\) as seeds for Modified Metropolis–Hastings (MMH) chains. Generate \(N\) new samples distributed approximately as \(\mathcal{N}(\mathbf{0}, \mathbf{I})\) conditioned on \(g \le y_{j-1}\). Set \(y_j\) as the \(p_0\)-th quantile of the new LSF values. Stop when \(y_j \le 0\).

  3. Final level — Count the actual failures (\(g \le 0\)) in the last conditional sample: \(\hat{p}_m = N_{\text{fail}} / N\).

  4. Estimate — \(\hat{p}_f = \hat{p}_1 \hat{p}_2 \cdots \hat{p}_m\)

Modified Metropolis–Hastings (MMH)

To generate samples from \(\mathcal{N}(\mathbf{0}, \mathbf{I})\) conditioned on \(g(\mathbf{u}) \le y_j\), the MMH algorithm applies Metropolis updates component-wise. For each component \(d\):

\[\xi_d \sim u_d + \sigma\, \mathcal{U}(-1, 1)\]

with acceptance probability

\[\alpha_d = \min\!\left(1,\; e^{-(\xi_d^2 - u_d^2)/2}\right)\]

After assembling all accepted components into a candidate \(\mathbf{u}'\), the entire vector is accepted only if \(g(\mathbf{u}') \le y_j\); otherwise the current state is retained. This ensures the stationary distribution is \(\mathcal{N}(\mathbf{0}, \mathbf{I}) \mid g(\mathbf{u}) \le y_j\).

Coefficient of variation

Ignoring correlations within the Markov chains (a lower bound on the true variance), the CoV of the estimator is approximated by [AuBeck2001]

\[\delta^2(\hat{p}_f) \approx \sum_{j=1}^{m} \frac{1 - \hat{p}_j}{N\,\hat{p}_j}\]

Subset Simulation is particularly effective for small failure probabilities (roughly \(p_f < 10^{-3}\)), where crude Monte Carlo would require an impractically large number of samples. A benchmark comparison of simulation methods on high-dimensional problems is given in [Schueller2007].

Use this method: Estimating failure probability by simulation · Simulation methods for reliability analysis · Reliability algorithms and results

For coordinate conventions, see Notation and probability spaces.