.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_reliability_sensitivity/reliability/plot_axial_stressed_beam.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note Click :ref:`here ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_reliability_sensitivity_reliability_plot_axial_stressed_beam.py: Axial stressed beam : comparing different methods to estimate a probability =========================================================================== .. GENERATED FROM PYTHON SOURCE LINES 6-13 In this example, we compare four methods to estimate the probability in the :ref:`axial stressed beam ` example : * Monte-Carlo simulation, * FORM, * directional sampling, * importance sampling with FORM design point: FORM-IS. .. GENERATED FROM PYTHON SOURCE LINES 15-17 Define the model ---------------- .. GENERATED FROM PYTHON SOURCE LINES 19-28 .. code-block:: default from __future__ import print_function import numpy as np from openturns.usecases import stressed_beam as stressed_beam import openturns as ot import openturns.viewer as viewer from matplotlib import pylab as plt ot.Log.Show(ot.Log.NONE) .. GENERATED FROM PYTHON SOURCE LINES 29-30 We load the model from the usecases module : .. GENERATED FROM PYTHON SOURCE LINES 30-32 .. code-block:: default sm = stressed_beam.AxialStressedBeam() .. GENERATED FROM PYTHON SOURCE LINES 33-34 The limit state function is defined in the `model` field of the data class : .. GENERATED FROM PYTHON SOURCE LINES 34-36 .. code-block:: default limitStateFunction = sm.model .. GENERATED FROM PYTHON SOURCE LINES 37-39 The probabilistic model of the axial stressed beam is defined in the data class. We get the first marginal and draw it : .. GENERATED FROM PYTHON SOURCE LINES 39-43 .. code-block:: default R_dist = sm.distribution_R graph = R_dist.drawPDF() view = viewer.View(graph) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_001.png :alt: plot axial stressed beam :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 44-45 We get the second marginal and draw it : .. GENERATED FROM PYTHON SOURCE LINES 47-51 .. code-block:: default F_dist = sm.distribution_F graph = F_dist.drawPDF() view = viewer.View(graph) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_002.png :alt: plot axial stressed beam :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 52-53 These independent marginals define the joint distribution of the input parameters : .. GENERATED FROM PYTHON SOURCE LINES 53-56 .. code-block:: default myDistribution = sm.distribution .. GENERATED FROM PYTHON SOURCE LINES 57-58 We create a `RandomVector` from the `Distribution`, then a composite random vector. Finally, we create a `ThresholdEvent` from this `RandomVector`. .. GENERATED FROM PYTHON SOURCE LINES 60-65 .. code-block:: default inputRandomVector = ot.RandomVector(myDistribution) outputRandomVector = ot.CompositeRandomVector( limitStateFunction, inputRandomVector) myEvent = ot.ThresholdEvent(outputRandomVector, ot.Less(), 0.0) .. GENERATED FROM PYTHON SOURCE LINES 66-68 Using Monte Carlo simulations ----------------------------- .. GENERATED FROM PYTHON SOURCE LINES 70-79 .. code-block:: default cv = 0.05 NbSim = 100000 experiment = ot.MonteCarloExperiment() algoMC = ot.ProbabilitySimulationAlgorithm(myEvent, experiment) algoMC.setMaximumOuterSampling(NbSim) algoMC.setBlockSize(1) algoMC.setMaximumCoefficientOfVariation(cv) .. GENERATED FROM PYTHON SOURCE LINES 80-81 For statistics about the algorithm .. GENERATED FROM PYTHON SOURCE LINES 81-83 .. code-block:: default initialNumberOfCall = limitStateFunction.getEvaluationCallsNumber() .. GENERATED FROM PYTHON SOURCE LINES 84-85 Perform the analysis. .. GENERATED FROM PYTHON SOURCE LINES 87-89 .. code-block:: default algoMC.run() .. GENERATED FROM PYTHON SOURCE LINES 90-99 .. code-block:: default result = algoMC.getResult() probabilityMonteCarlo = result.getProbabilityEstimate() numberOfFunctionEvaluationsMonteCarlo = limitStateFunction.getEvaluationCallsNumber() - \ initialNumberOfCall print('Number of calls to the limit state =', numberOfFunctionEvaluationsMonteCarlo) print('Pf = ', probabilityMonteCarlo) print('CV =', result.getCoefficientOfVariation()) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Number of calls to the limit state = 12693 Pf = 0.030568029622626627 CV = 0.04998535791736659 .. GENERATED FROM PYTHON SOURCE LINES 100-104 .. code-block:: default graph = algoMC.drawProbabilityConvergence() graph.setLogScale(ot.GraphImplementation.LOGX) view = viewer.View(graph) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_003.png :alt: ProbabilitySimulationAlgorithm convergence graph at level 0.95 :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 105-107 Using FORM analysis ------------------- .. GENERATED FROM PYTHON SOURCE LINES 109-110 We create a NearestPoint algorithm .. GENERATED FROM PYTHON SOURCE LINES 110-119 .. code-block:: default myCobyla = ot.Cobyla() # Resolution options: eps = 1e-3 myCobyla.setMaximumEvaluationNumber(100) myCobyla.setMaximumAbsoluteError(eps) myCobyla.setMaximumRelativeError(eps) myCobyla.setMaximumResidualError(eps) myCobyla.setMaximumConstraintError(eps) .. GENERATED FROM PYTHON SOURCE LINES 120-121 For statistics about the algorithm .. GENERATED FROM PYTHON SOURCE LINES 121-123 .. code-block:: default initialNumberOfCall = limitStateFunction.getEvaluationCallsNumber() .. GENERATED FROM PYTHON SOURCE LINES 124-125 We create a FORM algorithm. The first parameter is a NearestPointAlgorithm. The second parameter is an event. The third parameter is a starting point for the design point research. .. GENERATED FROM PYTHON SOURCE LINES 127-129 .. code-block:: default algoFORM = ot.FORM(myCobyla, myEvent, myDistribution.getMean()) .. GENERATED FROM PYTHON SOURCE LINES 130-131 Perform the analysis. .. GENERATED FROM PYTHON SOURCE LINES 133-135 .. code-block:: default algoFORM.run() .. GENERATED FROM PYTHON SOURCE LINES 136-143 .. code-block:: default resultFORM = algoFORM.getResult() numberOfFunctionEvaluationsFORM = limitStateFunction.getEvaluationCallsNumber() - \ initialNumberOfCall probabilityFORM = resultFORM.getEventProbability() print('Number of calls to the limit state =', numberOfFunctionEvaluationsFORM) print('Pf =', probabilityFORM) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Number of calls to the limit state = 98 Pf = 0.02998278558231473 .. GENERATED FROM PYTHON SOURCE LINES 144-147 .. code-block:: default graph = resultFORM.drawImportanceFactors() view = viewer.View(graph) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_004.png :alt: Importance Factors from Design Point - Unnamed :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 148-150 Using Directional sampling -------------------------- .. GENERATED FROM PYTHON SOURCE LINES 152-153 Resolution options: .. GENERATED FROM PYTHON SOURCE LINES 153-161 .. code-block:: default cv = 0.05 NbSim = 10000 algoDS = ot.DirectionalSampling(myEvent) algoDS.setMaximumOuterSampling(NbSim) algoDS.setBlockSize(1) algoDS.setMaximumCoefficientOfVariation(cv) .. GENERATED FROM PYTHON SOURCE LINES 162-163 For statistics about the algorithm .. GENERATED FROM PYTHON SOURCE LINES 163-165 .. code-block:: default initialNumberOfCall = limitStateFunction.getEvaluationCallsNumber() .. GENERATED FROM PYTHON SOURCE LINES 166-167 Perform the analysis. .. GENERATED FROM PYTHON SOURCE LINES 169-171 .. code-block:: default algoDS.run() .. GENERATED FROM PYTHON SOURCE LINES 172-181 .. code-block:: default result = algoDS.getResult() probabilityDirectionalSampling = result.getProbabilityEstimate() numberOfFunctionEvaluationsDirectionalSampling = limitStateFunction.getEvaluationCallsNumber() - \ initialNumberOfCall print('Number of calls to the limit state =', numberOfFunctionEvaluationsDirectionalSampling) print('Pf = ', probabilityDirectionalSampling) print('CV =', result.getCoefficientOfVariation()) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Number of calls to the limit state = 8926 Pf = 0.029480392428094374 CV = 0.049931236082334296 .. GENERATED FROM PYTHON SOURCE LINES 182-186 .. code-block:: default graph = algoDS.drawProbabilityConvergence() graph.setLogScale(ot.GraphImplementation.LOGX) view = viewer.View(graph) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_005.png :alt: DirectionalSampling convergence graph at level 0.95 :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 187-189 Using importance sampling with FORM design point: FORM-IS --------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 191-192 The `getStandardSpaceDesignPoint` method returns the design point in the U-space. .. GENERATED FROM PYTHON SOURCE LINES 194-197 .. code-block:: default standardSpaceDesignPoint = resultFORM.getStandardSpaceDesignPoint() standardSpaceDesignPoint .. raw:: html

[-1.59355,0.999463]



.. GENERATED FROM PYTHON SOURCE LINES 198-199 The key point is to define the importance distribution in the U-space. To define it, we use a multivariate standard Gaussian and configure it so that the center is equal to the design point in the U-space. .. GENERATED FROM PYTHON SOURCE LINES 201-204 .. code-block:: default dimension = myDistribution.getDimension() dimension .. rst-class:: sphx-glr-script-out Out: .. code-block:: none 2 .. GENERATED FROM PYTHON SOURCE LINES 205-209 .. code-block:: default myImportance = ot.Normal(dimension) myImportance.setMean(standardSpaceDesignPoint) myImportance .. raw:: html

Normal(mu = [-1.59355,0.999463], sigma = [1,1], R = [[ 1 0 ]
[ 0 1 ]])



.. GENERATED FROM PYTHON SOURCE LINES 210-211 Create the design of experiment corresponding to importance sampling. This generates a `WeightedExperiment` with weights corresponding to the importance distribution. .. GENERATED FROM PYTHON SOURCE LINES 213-215 .. code-block:: default experiment = ot.ImportanceSamplingExperiment(myImportance) .. GENERATED FROM PYTHON SOURCE LINES 216-217 Create the standard event corresponding to the event. This transforms the original problem into the U-space, with Gaussian independent marginals. .. GENERATED FROM PYTHON SOURCE LINES 219-221 .. code-block:: default standardEvent = ot.StandardEvent(myEvent) .. GENERATED FROM PYTHON SOURCE LINES 222-223 We then create the simulation algorithm. .. GENERATED FROM PYTHON SOURCE LINES 225-229 .. code-block:: default algo = ot.ProbabilitySimulationAlgorithm(standardEvent, experiment) algo.setMaximumCoefficientOfVariation(cv) algo.setMaximumOuterSampling(40000) .. GENERATED FROM PYTHON SOURCE LINES 230-231 For statistics about the algorithm .. GENERATED FROM PYTHON SOURCE LINES 231-233 .. code-block:: default initialNumberOfCall = limitStateFunction.getEvaluationCallsNumber() .. GENERATED FROM PYTHON SOURCE LINES 234-236 .. code-block:: default algo.run() .. GENERATED FROM PYTHON SOURCE LINES 237-238 retrieve results .. GENERATED FROM PYTHON SOURCE LINES 238-247 .. code-block:: default result = algo.getResult() probabilityFORMIS = result.getProbabilityEstimate() numberOfFunctionEvaluationsFORMIS = limitStateFunction.getEvaluationCallsNumber() - \ initialNumberOfCall print('Number of calls to the limit state =', numberOfFunctionEvaluationsFORMIS) print('Pf = ', probabilityFORMIS) print('CV =', result.getCoefficientOfVariation()) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Number of calls to the limit state = 976 Pf = 0.027514306856638082 CV = 0.049927184185779384 .. GENERATED FROM PYTHON SOURCE LINES 248-250 Conclusion ---------- .. GENERATED FROM PYTHON SOURCE LINES 252-253 We now compare the different methods in terms of accuracy and speed. .. GENERATED FROM PYTHON SOURCE LINES 258-259 The following function computes the number of correct base-10 digits in the computed result compared to the exact result. .. GENERATED FROM PYTHON SOURCE LINES 261-266 .. code-block:: default def computeLogRelativeError(exact, computed): logRelativeError = -np.log10(abs(exact - computed) / abs(exact)) return logRelativeError .. GENERATED FROM PYTHON SOURCE LINES 267-268 The following function prints the results. .. GENERATED FROM PYTHON SOURCE LINES 270-284 .. code-block:: default def printMethodSummary(name, computedProbability, numberOfFunctionEvaluations): print("---") print(name, ":") print('Number of calls to the limit state =', numberOfFunctionEvaluations) print('Pf = ', computedProbability) exactProbability = 0.02919819462483051 logRelativeError = computeLogRelativeError( exactProbability, computedProbability) print("Number of correct digits=%.3f" % (logRelativeError)) performance = logRelativeError/numberOfFunctionEvaluations print("Performance=%.2e (correct digits/evaluation)" % (performance)) return .. GENERATED FROM PYTHON SOURCE LINES 285-293 .. code-block:: default printMethodSummary("Monte-Carlo", probabilityMonteCarlo, numberOfFunctionEvaluationsMonteCarlo) printMethodSummary("FORM", probabilityFORM, numberOfFunctionEvaluationsFORM) printMethodSummary("DirectionalSampling", probabilityDirectionalSampling, numberOfFunctionEvaluationsDirectionalSampling) printMethodSummary("FORM-IS", probabilityFORMIS, numberOfFunctionEvaluationsFORMIS) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none --- Monte-Carlo : Number of calls to the limit state = 12693 Pf = 0.030568029622626627 Number of correct digits=1.329 Performance=1.05e-04 (correct digits/evaluation) --- FORM : Number of calls to the limit state = 98 Pf = 0.02998278558231473 Number of correct digits=1.571 Performance=1.60e-02 (correct digits/evaluation) --- DirectionalSampling : Number of calls to the limit state = 8926 Pf = 0.029480392428094374 Number of correct digits=2.015 Performance=2.26e-04 (correct digits/evaluation) --- FORM-IS : Number of calls to the limit state = 976 Pf = 0.027514306856638082 Number of correct digits=1.239 Performance=1.27e-03 (correct digits/evaluation) .. GENERATED FROM PYTHON SOURCE LINES 294-300 We see that all three methods produce the correct probability, but not with the same accuracy. In this case, we have found the correct order of magnitude of the probability, i.e. between one and two correct digits. There is, however, a significant difference in computational performance (measured here by the number of function evaluations). * The fastest method is the FORM method, which produces more than 1 correct digit with less than 98 function evaluations with a performance equal to :math:`1.60 \times 10^{-2}` (correct digits/evaluation). A practical limitation is that the FORM method does not produce a confidence interval: there is no guarantee that the computed probability is correct. * The slowest method is Monte-Carlo simulation, which produces more than 1 correct digit with 12806 function evaluations. This is associated with a very slow performance equal to :math:`1.11 \times 10^{-4}` (correct digits/evaluation). The interesting point with the Monte-Carlo simulation is that the method produces a confidence interval. * The DirectionalSampling method is somewhat in-between the two previous methods. * The FORM-IS method produces 2 correct digits and has a small number of function evaluations. It has an intermediate performance equal to :math:`2.37\times 10^{-3}` (correct digits/evaluation). It combines the best of the both worlds: it has the small number of function evaluation of FORM computation and the confidence interval of Monte-Carlo simulation. .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 0.894 seconds) .. _sphx_glr_download_auto_reliability_sensitivity_reliability_plot_axial_stressed_beam.py: .. only :: html .. container:: sphx-glr-footer :class: sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_axial_stressed_beam.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_axial_stressed_beam.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_