Estimate a probability with Monte CarloΒΆ

In this basic example we are going to estimate a probability by means of a simulation algorithm.

This model is a simple beam, restrained at one side and stressed by a traction load F at the other side.

The geometry is supposed to be deterministic: the diameter D is fixed to D=20 mm.

It is considered that failure occurs when the beam plastifies, i.e. when the axial stress gets bigger than the yield stress:

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

Therefore, the state limit G used here is:

G = \sigma_e - \frac{F}{\pi-D^2/4}

Two independent random variables R and S are considered:

  • R (the strength):

    R = \sigma_e

  • S (the load):

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

Stochastic model:

  • F ~ Normal(75e3, 5e3) [N]
  • R ~ LogNormal(300, 30) [N]

Theoretical results:

This two dimensional stochastic problem can be solved by calculating directly the failure probability:

P_f = \text{Prob}(R-S \leq 0) = \int_{r-s \leq 0} f_{R, S}(r, s)drds

If R and S are independant, then:


The numerical application gives:

P_f = 0.0292

In [10]:
from __future__ import print_function
import openturns as ot
In [11]:
# create the joint distribution of the parameters
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)
In [12]:
# create the model
model = ot.SymbolicFunction(['R', 'F'], ['R-F/(_pi*100.0)'])
In [13]:
# create the event we want to estimate the probability
vect = ot.RandomVector(distribution)
G = ot.RandomVector(model, vect)
event = ot.Event(G, ot.Less(), 0.0)
In [14]:
# create a Monte Carlo algorithm
algo = ot.MonteCarlo(event)
In [15]:
# retrieve results
result = algo.getResult()
probability = result.getProbabilityEstimate()
print('Pf=', probability)
Pf= 0.02858192505510653