Estimate a probability with Monte Carlo¶
In this example we estimate a probability by means of a simulation algorithm, the Monte-Carlo algorithm. To do this, we need the classes
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:
By definition, the yield stress is the load divided by the surface. Since the surface is , the stress is:
Failure occurs when the beam plastifies, i.e. when the axial stress gets larger than the yield stress:
where is the strength.
Therefore, the limit state function is:
for any .
The value of the parameter is such that:
which leads to the equation:
We consider the following distribution functions.
LogNormal(, ) [Pa]
Normal(, ) [N]
where and are the mean and the variance of .
The failure probability is:
The exact is
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 whose probability we want to estimate.
vect = ot.RandomVector(distribution) G = ot.CompositeRandomVector(model, vect) event = ot.ThresholdEvent(G, ot.Less(), 0.0)
Create a Monte Carlo algorithm.
experiment = ot.MonteCarloExperiment() algo = ot.ProbabilitySimulationAlgorithm(event, experiment) algo.setMaximumCoefficientOfVariation(0.05) algo.setMaximumOuterSampling(int(1e5)) algo.run()
result = algo.getResult() probability = result.getProbabilityEstimate() print('Pf=', probability)