Estimate a probability with Latin Hypercube Sampling

In this example we show how to use the LHS algorithm to estimate the probability of an event.

Introduction

We consider a simple beam stressed by a traction load F at both sides.

The geometry is supposed to be deterministic; the diameter D is equal to:

D=0.02 \textrm{ (m)}.

By definition, the yield stress is the load divided by the surface. Since the surface is \pi D^2/4, the stress is:

S = \frac{F}{\pi D^2/4}.

Failure occurs when the beam plastifies, i.e. when the axial stress gets larger than the yield stress:

R - \frac{F}{\pi D^2/4} \leq 0

where R is the strength.

Therefore, the limit state function G is:

G(R,F) = R - \frac{F}{\pi D^2/4},

for any R,F\in\mathbb{R}.

The value of the parameter D is such that:

D^2/4 = 10^{-4},

which leads to the equation:

G(R,F) = R - \frac{F}{10^{-4} \pi}.

We consider the following distribution functions.

Variable

Distribution

R

LogNormal(\mu_R=3\times 10^6, \sigma_R=3\times 10^5) [Pa]

F

Normal(\mu_F=750, \sigma_F=50) [N]

where \mu_R=E(R) and \sigma_R^2=V(R) are the mean and the variance of R.

The failure probability is:

P_f = \text{Prob}(G(R,F) \leq 0).

The exact P_f is

P_f = 0.02920.

[1]:

from __future__ import print_function
import openturns as ot

Create the joint distribution of the parameters.

[2]:

distribution_R = ot.LogNormalMuSigma(300.0, 30.0, 0.0).getDistribution()
distribution_F = ot.Normal(75e3, 5e3)
marginals = [distribution_R, distribution_F]
distribution = ot.ComposedDistribution(marginals)

Create the model.

[3]:

model = ot.SymbolicFunction(['R', 'F'], ['R-F/(pi_*100.0)'])

Create the event whose probability we want to estimate.

[4]:

vect = ot.RandomVector(distribution)
G = ot.CompositeRandomVector(model, vect)
event = ot.ThresholdEvent(G, ot.Less(), 0.0)

Create a LHS algorithm.

[5]:

algo = ot.LHS(event)
algo.setMaximumCoefficientOfVariation(0.05)
algo.setMaximumOuterSampling(int(1e5))
algo.run()

Retrieve results.

[6]:

result = algo.getResult()
probability = result.getProbabilityEstimate()
print('Pf=', probability)

Pf= 0.029342988609791055