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:
Therefore, the state limit G used here is:
Two independent random variables R and S are considered:
R (the strength):
S (the load):
- F ~ Normal(75e3, 5e3) [N]
- R ~ LogNormal(300, 30) [N]
This two dimensional stochastic problem can be solved by calculating directly the failure probability:
If R and S are independant, then:
The numerical application gives:
from __future__ import print_function import openturns as ot
# 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)
# create the model model = ot.SymbolicFunction(['R', 'F'], ['R-F/(_pi*100.0)'])
# 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)
# create a Monte Carlo algorithm algo = ot.MonteCarlo(event) algo.setMaximumCoefficientOfVariation(0.05) algo.setMaximumOuterSampling(int(1e5)) algo.run()
# retrieve results result = algo.getResult() probability = result.getProbabilityEstimate() print('Pf=', probability)