Probability estimation: directional simulation

In this basic example we are going to estimate a failure probability with the directional simulation.

The directional simulation algorithm operates in the standard space and is parametrized by two concepts:

  1. a Root Strategy to evaluate the intersections of each direction with the limit state function and take into account the contribution of the direction to the event probability. The available corresponding classes are:

    • RiskyAndFast

    • MediumSafe

    • SafeAndSlow

  2. a Sampling Strategy to chose directions in the standard space. The available corresponding classes are:

    • RandomDirection

    • OrthogonalDirection

Let’s consider the following analytical example of a cantilever beam, of Young’s modulus E, length L, section modulus I.

One end is built in a wall and we apply a concentrated bending load F at the other end of the beam, resulting in a deviation:

d = \frac{F*L^3}{3*E*I}

It is considered that failure occurs when the beam deviation is too important:

d \ge 30 (cm)

Four independent random variables are considered:

  • E (the Young’s modulus) [Pa]

  • F (the load) [N]

  • L (the length) # [m]

  • I (the section) # [m^4]

Stochastic model (simplified model, no units):

  • E ~ Beta(0.93, 3.2, 2.8e7, 4.8e7)

  • F ~ LogNormal(30000, 9000, 15000)

  • L ~ Uniform(250, 260)

  • I ~ Beta(2.5, 4.0, 3.1e2, 4.5e2)

from __future__ import print_function
import openturns as ot
# Create the marginal distributions of the parameters
dist_E = ot.Beta(0.93, 2.27, 2.8e7, 4.8e7)
dist_F = ot.LogNormalMuSigma(30000, 9000, 15000).getDistribution()
dist_L = ot.Uniform(250, 260)
dist_I = ot.Beta(2.5, 1.5, 3.1e2, 4.5e2)
marginals = [dist_E, dist_F, dist_L, dist_I]
# Create the Copula
RS = ot.CorrelationMatrix(4)
RS[2, 3] = -0.2
# Evaluate the correlation matrix of the Normal copula from RS
R = ot.NormalCopula.GetCorrelationFromSpearmanCorrelation(RS)
# Create the Normal copula parametrized by R
copula = ot.NormalCopula(R)
# Create the joint probability distribution
distribution = ot.ComposedDistribution(marginals, copula)
# create the model
model = ot.SymbolicFunction(['E', 'F', 'L', 'I'], ['F*L^3/(3*E*I)'])
# create the event we want to estimate the probability
vect = ot.RandomVector(distribution)
G = ot.CompositeRandomVector(model, vect)
event = ot.ThresholdEvent(G, ot.Greater(), 30.0)
# root finding algorithm
solver = ot.Brent()
rootStrategy = ot.MediumSafe(solver)
# direction sampling algorithm
samplingStrategy = ot.OrthogonalDirection()
# Create a simulation algorithm
algo = ot.DirectionalSampling(event, rootStrategy, samplingStrategy)
# retrieve results
result = algo.getResult()
probability = result.getProbabilityEstimate()
print('Pf=', probability)
Pf= 0.006797156019437116