.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_reliability_sensitivity/reliability/plot_form_explained.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_form_explained.py: An illustrated example of a FORM probability estimate ===================================================== .. GENERATED FROM PYTHON SOURCE LINES 7-14 Abstract -------- In this example we illustrate the different steps of a FORM/SORM analysis on a simple example. We focus on the different steps and compare them with an analytic computation whenever possible. .. GENERATED FROM PYTHON SOURCE LINES 14-19 .. code-block:: default import openturns as ot import openturns.viewer as otv import numpy as np from matplotlib import pylab as plt .. GENERATED FROM PYTHON SOURCE LINES 20-30 Position of the problem ----------------------- We consider a bivariate random vector :math:`X = (X_1, X_2)` with the following independent marginals : - an exponential distribution with parameter :math:`\lambda=1`, :math:`X_1 \sim \mathcal{E}(1.0)` ; - a standard unit gaussian :math:`X_2 \sim \mathcal{N}(0,1)`. The support of the input vector is :math:`[0, +\infty[ \times \mathbb{R}` .. GENERATED FROM PYTHON SOURCE LINES 32-37 .. code-block:: default distX1 = ot.Exponential(1.0) distX2 = ot.Normal() distX = ot.ComposedDistribution([distX1, distX2]) .. GENERATED FROM PYTHON SOURCE LINES 38-39 We can draw the bidimensional PDF of the distribution `distX` over :math:`[0,-10] \times [10,10]` : .. GENERATED FROM PYTHON SOURCE LINES 39-48 .. code-block:: default ot.ResourceMap_SetAsUnsignedInteger("Contour-DefaultLevelsNumber", 8) graphPDF = distX.drawPDF([0, -10], [10, 10]) graphPDF.setTitle(r'2D-PDF of the input variables $(X_1, X_2)$') graphPDF.setXTitle(r'$x_1$') graphPDF.setYTitle(r'$x_2$') graphPDF.setLegendPosition("bottomright") view = otv.View(graphPDF) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_001.png :alt: 2D-PDF of the input variables $(X_1, X_2)$ :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 49-51 We consider the model :math:`f : (x_1, x_2) \mapsto x_1 x_2` which maps the random input vector :math:`X` to the output variable :math:`Y=f(X) \in \mathbb{R}`. We also draw the isolines of the model `f`. .. GENERATED FROM PYTHON SOURCE LINES 51-59 .. code-block:: default f = ot.SymbolicFunction(['x1', 'x2'], ['x1 * x2']) graphModel = f.draw([0.0, -10.0], [10.0, 10.0]) graphModel.setXTitle(r'$x_1$') graphModel.setXTitle(r'$x_2$') graphModel.setTitle(r'Isolines of the model : $Y = f(X)$') view = otv.View(graphModel) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_002.png :alt: Isolines of the model : $Y = f(X)$ :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 60-68 We want to estimate the probability :math:`P_f` of the output variable to be greater than a prescribed threshold :math:`s=10` : this is the failure event. This probability is simply expressed as an integral : .. math:: P_f = \int_{\mathcal{D}} \mathbf{1}_{\mathcal{D}}(x) df_{X_1,X_2}(x) where :math:`\mathcal{D} = \{ (x_1, x_2) \in [0,+\infty[ \times \mathbb{R} / x_1 x_2 \geq s \}` is the failure domain. In the general case the probability density function :math:`f_{X_1,X_2}` and the domain of integration :math:`\mathcal{D}` are difficult to handle. .. GENERATED FROM PYTHON SOURCE LINES 70-71 We first define RandomVector objects and the failure event associated to the ouput random variable. .. GENERATED FROM PYTHON SOURCE LINES 71-77 .. code-block:: default vectorX = ot.RandomVector(distX) vectorY = ot.CompositeRandomVector(f, vectorX) s = 10.0 event = ot.ThresholdEvent(vectorY, ot.Greater(), s) .. GENERATED FROM PYTHON SOURCE LINES 78-85 This event can easily be represented with a 1D curve as it is a branch of an hyperbole. If :math:`y = x_1 x_2 = 10.0`, then the boundary of the domain of failure is the curve : .. math:: h : x_1 \mapsto \frac{10.0}{x_1} .. GENERATED FROM PYTHON SOURCE LINES 87-88 We shall represent this curve using a :class:`~openturns.Contour` object. .. GENERATED FROM PYTHON SOURCE LINES 88-99 .. code-block:: default nx, ny = 15, 15 xx = ot.Box([nx], ot.Interval([0.0], [10.0])).generate() yy = ot.Box([ny], ot.Interval([-10.0], [10.0])).generate() inputData = ot.Box([nx, ny], ot.Interval( [0.0, -10.0], [10.0, 10.0])).generate() outputData = f(inputData) mycontour = ot.Contour(xx, yy, outputData, [10.0], ["10.0"]) myGraph = ot.Graph("Representation of the failure domain", r"$X_1$", r"$X_2$", True, "") myGraph.add(mycontour) .. GENERATED FROM PYTHON SOURCE LINES 100-106 .. code-block:: default texts = [r" Event : $\mathcal{D} = \{Y \geq 10.0\}$"] myText = ot.Text([[4.0, 4.0]], texts) myText.setTextSize(1) myGraph.add(myText) view = otv.View(myGraph) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_003.png :alt: Representation of the failure domain :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 107-109 We can superimpose the event boundary with the 2D-PDF ot the input variables : .. GENERATED FROM PYTHON SOURCE LINES 109-115 .. code-block:: default mycontour.setColor("black") mycontour.setLabels(["event"]) graphPDF.add(mycontour) graphPDF.setLegendPosition("bottomright") view = otv.View(graphPDF) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_004.png :alt: 2D-PDF of the input variables $(X_1, X_2)$ :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 116-118 From the previous figure we observe that in the failure domain the PDF takes small (and even very small) values. Consequently the probability of the failure, the integral :math:`P_f` is also expected to be small. The FORM/SORM methods estimate this kind of integral. .. GENERATED FROM PYTHON SOURCE LINES 120-136 The FORM approximation ---------------------- The basic steps of the FORM (or SORM) algorithm are : - an isoprobabilistic transform ; - finding the design point : that is the nearest point wrt the origin in the standard space ; - estimating the probability integral. As mentionned, both the density function and the domain of integration are complex in general. The first step of the FORM method makes the density easier to work with and the second step tackles the domain of integration problem. Variable transform ^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 138-141 OpenTURNS has several isoprobabilistic transforms and the FORM/SORM classes implement the Generalized Nataf and Rosenblatt transforms. In this case the `distX` distribution is not elliptical so the default method is the Rosenblatt transform. .. GENERATED FROM PYTHON SOURCE LINES 141-143 .. code-block:: default print("Is Elliptical ? ", distX.isElliptical()) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Is Elliptical ? False .. GENERATED FROM PYTHON SOURCE LINES 144-152 We seek an isoprobabilistic transform :math:`T` such as .. math:: T : X \mapsto Z where each component of :math:`Z` is a standard unit gaussian. The isoprobabilistic transform and its inverse are methods of the distribution `distX` : .. GENERATED FROM PYTHON SOURCE LINES 152-155 .. code-block:: default transformation = distX.getIsoProbabilisticTransformation() inverseTransformation = distX.getInverseIsoProbabilisticTransformation() .. GENERATED FROM PYTHON SOURCE LINES 156-160 The main goal of this step is to work with a simpler probability density function of the input variables as they will be standard gaussian unit and uncorrelated. The domain of integration will still be complicated but will be handled with a well chosen approximate. .. GENERATED FROM PYTHON SOURCE LINES 162-167 We detail the Rosenblatt transform in this simple case. In this example we consider independent variables so the transform is simpler, we only have to perform the transformation on each variable. For the second one is already a standard unit gaussian we transform the first variable only. .. GENERATED FROM PYTHON SOURCE LINES 169-171 We draw a realization of the random input vector. This point is said to be in the physical space. We shall focus on the first component. .. GENERATED FROM PYTHON SOURCE LINES 171-173 .. code-block:: default xi = vectorX.getRealization() .. GENERATED FROM PYTHON SOURCE LINES 174-178 The first step of the Rosenblatt transform is to build a random variable :math:`u` with a uniform law in ]0,1[. This is done through an evaluation of the CDF of `distX1` at the given point in the physical space. Once again, please note that the second component is left unchanged. .. GENERATED FROM PYTHON SOURCE LINES 178-180 .. code-block:: default ui = [distX1.computeCDF(xi[0]), xi[1]] .. GENERATED FROM PYTHON SOURCE LINES 181-183 The second step is to build a standard unit gaussian from a uniform variable. This is done by a simple call to the probit function. The point `zi` is said to be in the standard space. .. GENERATED FROM PYTHON SOURCE LINES 183-185 .. code-block:: default zi = [-ot.Normal().computeInverseSurvivalFunction(ui[0])[0], ui[1]] .. GENERATED FROM PYTHON SOURCE LINES 186-187 The sought transform then maps a point in the physical space to the standard space : .. GENERATED FROM PYTHON SOURCE LINES 187-190 .. code-block:: default print(xi, "->", ui, "->", zi) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none [0.291882,-0.55673] -> [0.25314355280016465, -0.5567298885851408] -> [-0.6646301103295218, -0.5567298885851408] .. GENERATED FROM PYTHON SOURCE LINES 191-193 We also build the isoprobabilistic transform :math:`T_1` and its inverse :math:`T_1^{-1}` for the first marginal : .. GENERATED FROM PYTHON SOURCE LINES 193-196 .. code-block:: default transformX1 = distX1.getIsoProbabilisticTransformation() inverseTransformX1 = distX1.getInverseIsoProbabilisticTransformation() .. GENERATED FROM PYTHON SOURCE LINES 197-203 We can check the result of our experiment against : - the 2D-transform :math:`T` ; - the 1D-transform :math:`T_1` and the second component unchanged ; and observe the results are the same. .. GENERATED FROM PYTHON SOURCE LINES 203-210 .. code-block:: default zi1D = [transformX1([xi[0]])[0], xi[1]] zi2D = transformation(xi) print("zi = ", zi) print("zi1D = ", zi1D) print("zi2D = ", zi2D) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none zi = [-0.6646301103295218, -0.5567298885851408] zi1D = [-0.6646301103295218, -0.5567298885851408] zi2D = [-0.66463,-0.55673] .. GENERATED FROM PYTHON SOURCE LINES 211-214 We can represent the boundary of the event in the standard space : that is a composition of the hyperbole :math:`h : x \mapsto 10/x` and the inverse transform :math:`T_1^{-1}` defined by :math:`inverseTransformX1`. .. GENERATED FROM PYTHON SOURCE LINES 214-228 .. code-block:: default failureBoundaryPhysicalSpace = ot.SymbolicFunction(['x'], ['10.0 / x']) failureBoundaryStandardSpace = ot.ComposedFunction( failureBoundaryPhysicalSpace, inverseTransformX1) x = np.linspace(1.1, 5.0, 100) cx = np.array([failureBoundaryStandardSpace([xi])[0] for xi in x]) graphStandardSpace = ot.Graph( 'Failure event in the standard space', r'$u_1$', r'$u_2$', True, '') curveCX = ot.Curve(x, cx, 'Boundary of the event $\partial \mathcal{D}$') curveCX.setLineStyle("solid") curveCX.setColor("blue") graphStandardSpace.add(curveCX) .. GENERATED FROM PYTHON SOURCE LINES 229-230 We add the origin to the previous graph. .. GENERATED FROM PYTHON SOURCE LINES 230-245 .. code-block:: default cloud = ot.Cloud([0.0], [0.0]) cloud.setColor("black") cloud.setPointStyle("fcircle") cloud.setLegend("origin") graphStandardSpace.add(cloud) graphStandardSpace.setGrid(True) graphStandardSpace.setLegendPosition("bottomright") # Some annotation texts = [r"Event : $\mathcal{D} = \{Y \geq 10.0\}$"] myText = ot.Text([[3.0, 4.0]], texts) myText.setTextSize(1) graphStandardSpace.add(myText) view = otv.View(graphStandardSpace) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_005.png :alt: Failure event in the standard space :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 246-254 The design point ^^^^^^^^^^^^^^^^ The FORM and SORM methods assume that the failure probability integral has its support in the vicinity of the closest point of the domain to the origin. Then the second step of the method is to find this point, *the design point*, through a minimization problem under constraints. .. GENERATED FROM PYTHON SOURCE LINES 256-257 We configure the Cobyla solver that we use for the optimization : .. GENERATED FROM PYTHON SOURCE LINES 257-264 .. code-block:: default solver = ot.Cobyla() solver.setMaximumIterationNumber(10000) solver.setMaximumAbsoluteError(1.0e-3) solver.setMaximumRelativeError(1.0e-3) solver.setMaximumResidualError(1.0e-3) solver.setMaximumConstraintError(1.0e-3) .. GENERATED FROM PYTHON SOURCE LINES 265-267 We build the FORM algorithm with its basic constructor. The starting point for the optimization algorithm is the mean of the input variables. .. GENERATED FROM PYTHON SOURCE LINES 267-269 .. code-block:: default algoFORM = ot.FORM(solver, event, distX.getMean()) .. GENERATED FROM PYTHON SOURCE LINES 270-271 We are ready to run the algorithm and store the result : .. GENERATED FROM PYTHON SOURCE LINES 271-274 .. code-block:: default algoFORM.run() result = algoFORM.getResult() .. GENERATED FROM PYTHON SOURCE LINES 275-277 The design point can be retrieved in both physical and standard space with respectively the `getPhysicalSpaceDesignPoint` and `getStandardSpaceDesignPoint` methods. .. GENERATED FROM PYTHON SOURCE LINES 277-283 .. code-block:: default designPointPhysicalSpace = result.getPhysicalSpaceDesignPoint() designPointStandardSpace = result.getStandardSpaceDesignPoint() print("Design point in physical space : ", designPointPhysicalSpace) print("Design point in standard space : ", designPointStandardSpace) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Design point in physical space : [4.84183,2.06513] Design point in standard space : [2.41384,2.06513] .. GENERATED FROM PYTHON SOURCE LINES 284-285 We can get the Hasofer index with the `getHasoferReliabilityIndex` method which is the distance of the design point to the origin : .. GENERATED FROM PYTHON SOURCE LINES 285-289 .. code-block:: default betaHL = result.getHasoferReliabilityIndex() print("Hasofer index : ", betaHL) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Hasofer index : 3.176696193499824 .. GENERATED FROM PYTHON SOURCE LINES 290-291 We visualize it on the previous graph. .. GENERATED FROM PYTHON SOURCE LINES 291-305 .. code-block:: default cloud = ot.Cloud([designPointStandardSpace[0]], [designPointStandardSpace[1]]) cloud.setColor("red") cloud.setPointStyle("fcircle") cloud.setLegend("design point") graphStandardSpace.add(cloud) graphStandardSpace.setGrid(True) graphStandardSpace.setLegendPosition("bottomright") cc = ot.Curve([0.0, designPointStandardSpace[0]], [ 0.0, designPointStandardSpace[1]], r'$\beta_{HL}$ distance') cc.setLineStyle("dashed") cc.setColor("black") graphStandardSpace.add(cc) view = otv.View(graphStandardSpace) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_006.png :alt: Failure event in the standard space :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 306-317 Estimating the failure probability integral ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The last step of the FORM algorithm is to replace the domain of integration by the half-space at the design point. In this simple example the half-space is delimited by the tangent at the design point in the standard space. The expression of the failure domain boundary in the standard space is the composition of the hyperbole :math:`h:x \mapsto 10/x` and the inverse transform on the first variable. We can compute the gradient (here the first derivative of a 1D function :math:`h(u_0)` ) at any given point with the getGradient method : .. GENERATED FROM PYTHON SOURCE LINES 319-327 .. code-block:: default u0 = [designPointStandardSpace[0]] du0 = failureBoundaryStandardSpace.getGradient().gradient(u0) print("abscissa of the design point u0 = ", u0[0]) print("value of the failure boundary at u0 = ", failureBoundaryStandardSpace(u0)[0]) print("value of the gradient of the failure boundary at u0 = ", du0[0, 0]) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none abscissa of the design point u0 = 2.4138442539794656 value of the failure boundary at u0 = 2.065335164461817 value of the gradient of the failure boundary at u0 = -1.1706609709100166 .. GENERATED FROM PYTHON SOURCE LINES 328-334 In the standard space the equation of the tangent :math:`\Pi_{u_0}(x)` is given by .. math:: \Pi_{u_0}(x) = (h \circ T^{-1}) (u_0) + \frac{d}{dx} (h \circ T^{-1}) (u_0) (x-u_0) .. GENERATED FROM PYTHON SOURCE LINES 334-342 .. code-block:: default x = np.linspace(1.1, 5.0, 100) hyperplane = failureBoundaryStandardSpace(u0)[0] + du0[0, 0] * (x-u0) curveHyperplane = ot.Curve(x, hyperplane, r'$\Pi_{u_0}$ (FORM)') curveHyperplane.setLineStyle("dashed") curveHyperplane.setColor("green") graphStandardSpace.add(curveHyperplane) view = otv.View(graphStandardSpace) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_007.png :alt: Failure event in the standard space :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 343-351 In the standard space the PDF of the input variables is rotationally invariant so .. math:: P_f \approx E(\beta_{HL}), where :math:`E(.)` is the survival function of the standard unit gaussian. .. GENERATED FROM PYTHON SOURCE LINES 351-354 .. code-block:: default pf = ot.Normal().computeSurvivalFunction(betaHL) print("FORM : Pf = ", pf) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none FORM : Pf = 0.0007448149708283453 .. GENERATED FROM PYTHON SOURCE LINES 355-356 This proability of failure is the one computed in the FORMResult and obtained with the `getEventProbability` method : .. GENERATED FROM PYTHON SOURCE LINES 356-360 .. code-block:: default pf = result.getEventProbability() print("Probability of failure (FORM) Pf = ", pf) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Probability of failure (FORM) Pf = 0.0007448149708283453 .. GENERATED FROM PYTHON SOURCE LINES 361-369 The SORM approximation ---------------------- The SORM approximate uses an osculating paraboloid instead of the half-space delimited by the tangent at the design point. In this case it is a simple parabola we can obtain through Taylor expansion at the design point. However, in the general case one has to manipulate the gradient and the hessian in the standard space which is cumbersome. .. GENERATED FROM PYTHON SOURCE LINES 371-373 We need the value of the second derivative of the failure boundary function at the design point in the standard space : .. GENERATED FROM PYTHON SOURCE LINES 373-379 .. code-block:: default u0 = [designPointStandardSpace[0]] d2u0 = failureBoundaryStandardSpace.getHessian().hessian(u0) print("abscissa of the design point u0 = ", u0[0]) print("value of the hessian of the failure boundary at u0 = ", d2u0[0, 0, 0]) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none abscissa of the design point u0 = 2.4138442539794656 value of the hessian of the failure boundary at u0 = 0.9401058369642105 .. GENERATED FROM PYTHON SOURCE LINES 380-386 In the standard space the equation of the osculating parabola :math:`\mathcal{P}_{u_0}(x)` at :math:`u_0` is given by .. math:: \mathcal{P}_{u_0}(x) = h \circ T^{-1} (u_0) + \frac{d}{dx} (h \circ T^{-1})(u_0) (x-u_0) + \frac{1}{2} \frac{d^2}{dx^2} (h \circ T^{-1})(u_0) (x-u_0)^2 .. GENERATED FROM PYTHON SOURCE LINES 386-396 .. code-block:: default x = np.linspace(1.1, 5.0, 100) parabola = failureBoundaryStandardSpace( u0)[0] + du0[0, 0] * (x-u0) + 0.5 * d2u0[0, 0, 0] * (x-u0)**2 curveParabola = ot.Curve(x, parabola, r'$\mathcal{P}_{u_0}$ (SORM)') curveParabola.setLineStyle("dashed") curveParabola.setColor("orange") graphStandardSpace.add(curveParabola) view = otv.View(graphStandardSpace) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_008.png :alt: Failure event in the standard space :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_008.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 397-407 The next step is to estimate the principal curvatures of the osculating paraboloid. For any regular function :math:`g` the curvature :math:`\kappa(x_0)` at the point :math:`x_0` in cartesian coordinates reads as .. math:: \kappa(x_0) = \frac{g''(x_0)}{(1+[g'(x_0)]^2)^{3/2}}. For the oscilating parabola of concern we use the gradient and hessian previously computed : .. GENERATED FROM PYTHON SOURCE LINES 407-411 .. code-block:: default curvature = (d2u0[0, 0, 0]) / (1 + (du0[0, 0]) ** 2)**(3/2) print("Curvature (analytic formula) = ", curvature) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Curvature (analytic formula) = 0.2575913913885428 .. GENERATED FROM PYTHON SOURCE LINES 412-413 We build the SORM algorithm and run it : .. GENERATED FROM PYTHON SOURCE LINES 413-416 .. code-block:: default algoSORM = ot.SORM(solver, event, distX.getMean()) algoSORM.run() .. GENERATED FROM PYTHON SOURCE LINES 417-418 The SORM result is obtained with the `getResult` method : .. GENERATED FROM PYTHON SOURCE LINES 418-420 .. code-block:: default resultSORM = algoSORM.getResult() .. GENERATED FROM PYTHON SOURCE LINES 421-423 The principal curvatures of the osculating paraboloid at the design point is obtained by the `getSortedCurvatures` method : .. GENERATED FROM PYTHON SOURCE LINES 423-426 .. code-block:: default print("Curvature (estimated) = ", resultSORM.getSortedCurvatures()[1]) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Curvature (estimated) = 0.25761034541532546 .. GENERATED FROM PYTHON SOURCE LINES 427-435 Once the curvature is obtained there are several ways of approximating the failure probability :math:`P_f`. OpenTURNS implements the Breitung, Hohenbichler and Tvedt estimates. For instance, the Breitung approximation gives .. math:: P_f \approx E(\beta_{HL}) \frac{1}{\sqrt{1+\beta_{HL}\kappa}} .. GENERATED FROM PYTHON SOURCE LINES 435-439 .. code-block:: default coeff = (1.0 + betaHL*curvature) ** (-0.5) pf = (1.0 - ot.Normal().computeCDF(betaHL)) * coeff print("SORM : Pf = ", pf) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none SORM : Pf = 0.0005523531956150853 .. GENERATED FROM PYTHON SOURCE LINES 440-441 We can compare with the different estimators implemented in OpenTURNS : .. GENERATED FROM PYTHON SOURCE LINES 441-449 .. code-block:: default pfBreitung = resultSORM.getEventProbabilityBreitung() pfHohenbichler = resultSORM.getEventProbabilityHohenbichler() pfTvedt = resultSORM.getEventProbabilityTvedt() print("Probability of failure (SORM Breintung) Pf = ", pfBreitung) print("Probability of failure (SORM Hohenbichler) Pf = ", pfHohenbichler) print("Probability of failure (SORM Tvedt) Pf = ", pfTvedt) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Probability of failure (SORM Breintung) Pf = 0.0005523440504782278 Probability of failure (SORM Hohenbichler) Pf = 0.0005420328660296243 Probability of failure (SORM Tvedt) Pf = 0.0005381057564251503 .. GENERATED FROM PYTHON SOURCE LINES 450-451 Display all figures .. GENERATED FROM PYTHON SOURCE LINES 451-452 .. code-block:: default otv.View.ShowAll() .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 1.305 seconds) .. _sphx_glr_download_auto_reliability_sensitivity_reliability_plot_form_explained.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_form_explained.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_form_explained.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_