Probability estimation with importance sampling simulation on cantilever beam example

In this example we estimate a failure probability with the importance sampling simulation algorithm provided by the ImportanceSamplingExperiment class.

The main steps of this method are:

  • run a FORM analysis,

  • create an importance distribution based on the results of the FORM results,

  • run a sampling-based probability estimate algorithm.

Description of the problem

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

One end of the cantilever beam 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{FL^3}{3EI}

Failure occurs when the beam deviation is too large:

d \ge 30 (cm)

Four independent random variables are considered:

  • E: Young’s modulus [Pa]

  • F: load [N]

  • L: length [m]

  • I: section [m^4]

Stochastic model (simplified model, no units):

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

  • F ~ LogNormal(30000, 9000, 15000)

  • L ~ Uniform(250, 260)

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

Define the cantilever beam model

[1]:

from __future__ import print_function
import openturns as ot

[2]:

# 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]

[3]:

# 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)

[4]:

# Create the joint probability distribution
distribution = ot.ComposedDistribution(marginals, copula)

[5]:

# create the model
model = ot.SymbolicFunction(['E', 'F', 'L', 'I'], ['F*L^3/(3*E*I)'])

Define the event

We create the event whose probability we want to estimate.

[6]:

vect = ot.RandomVector(distribution)
G = ot.CompositeRandomVector(model, vect)
event = ot.ThresholdEvent(G, ot.Greater(), 30.0)

Run a FORM analysis

[7]:

# Define a solver
optimAlgo = ot.Cobyla()
optimAlgo.setMaximumEvaluationNumber(1000)
optimAlgo.setMaximumAbsoluteError(1.0e-10)
optimAlgo.setMaximumRelativeError(1.0e-10)
optimAlgo.setMaximumResidualError(1.0e-10)
optimAlgo.setMaximumConstraintError(1.0e-10)

[8]:

# Run FORM
algo = ot.FORM(optimAlgo, event, distribution.getMean())
algo.run()
result = algo.getResult()

Configure an importance sampling algorithm

The getStandardSpaceDesignPoint method returns the design point in the U-space.

[9]:

standardSpaceDesignPoint = result.getStandardSpaceDesignPoint()
standardSpaceDesignPoint

[9]:

[-0.602386,2.31056,0.355794,-0.533677]

The key point is to define the importance distribution in the U-space. To define it, we use a multivariate standard Gaussian centered around the design point in the U-space.

[10]:

dimension = distribution.getDimension()
dimension

[10]:

4

[11]:

myImportance = ot.Normal(dimension)
myImportance.setMean(standardSpaceDesignPoint)
myImportance

[11]:

Normal(mu = [-0.602386,2.31056,0.355794,-0.533677], sigma = [1,1,1,1], R = [[ 1 0 0 0 ]
[ 0 1 0 0 ]
[ 0 0 1 0 ]
[ 0 0 0 1 ]])

Create the design of experiment corresponding to importance sampling. This generates a WeightedExperiment with weights fitting to the importance distribution.

[12]:

experiment = ot.ImportanceSamplingExperiment(myImportance)
type(experiment)

[12]:

openturns.weightedexperiment.ImportanceSamplingExperiment

Create the standard event corresponding to the event. This pushes the original problem to the U-space, with Gaussian independent marginals.

[13]:

standardEvent = ot.StandardEvent(event)

Run the importance sampling simulation

We then create the simulation algorithm.

[14]:

algo = ot.ProbabilitySimulationAlgorithm(standardEvent, experiment)
algo.setMaximumCoefficientOfVariation(0.1)
algo.setMaximumOuterSampling(40000)
algo.run()

[15]:

# retrieve results
result = algo.getResult()
probability = result.getProbabilityEstimate()
probability

[15]:

0.0056136313619849386

In order to compute the confidence interval, we use the getConfidenceLength method, which returns the length of the interval. In order to compute the bounds of the interval, we divide this length by 2.

[16]:

alpha = 0.05

[17]:

pflen = result.getConfidenceLength(1-alpha)
print("%.2f%% confidence interval = [%f,%f]" % ((1-alpha)*100,probability-pflen/2,probability+pflen/2))

95.00% confidence interval = [0.004516,0.006711]