.. 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 :ref:`Go to the end ` 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-18 .. code-block:: Python import openturns as ot import openturns.viewer as otv import numpy as np .. GENERATED FROM PYTHON SOURCE LINES 19-29 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 31-36 .. code-block:: Python distX1 = ot.Exponential(1.0) distX2 = ot.Normal() distX = ot.ComposedDistribution([distX1, distX2]) .. GENERATED FROM PYTHON SOURCE LINES 37-38 We can draw the bidimensional PDF of the distribution `distX` over :math:`[0,-10] \times [10,10]`: .. GENERATED FROM PYTHON SOURCE LINES 38-47 .. code-block:: Python 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("lower right") 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 48-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:: Python 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-69 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 71-72 We first define RandomVector objects and the failure event associated to the output random variable. .. GENERATED FROM PYTHON SOURCE LINES 72-78 .. code-block:: Python vectorX = ot.RandomVector(distX) vectorY = ot.CompositeRandomVector(f, vectorX) s = 10.0 event = ot.ThresholdEvent(vectorY, ot.Greater(), s) .. GENERATED FROM PYTHON SOURCE LINES 79-86 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 88-89 We shall represent this curve using a :class:`~openturns.Contour` object. .. GENERATED FROM PYTHON SOURCE LINES 89-98 .. code-block:: Python 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 99-105 .. code-block:: Python 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 106-108 We can superimpose the event boundary with the 2D-PDF ot the input variables : .. GENERATED FROM PYTHON SOURCE LINES 108-114 .. code-block:: Python mycontour.setColor("black") mycontour.setLabels(["event"]) graphPDF.add(mycontour) graphPDF.setLegendPosition("lower right") 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 115-119 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 121-137 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 mentioned, 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 139-142 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 142-144 .. code-block:: Python print("Is Elliptical ? ", distX.isElliptical()) .. rst-class:: sphx-glr-script-out .. code-block:: none Is Elliptical ? False .. GENERATED FROM PYTHON SOURCE LINES 145-153 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 153-156 .. code-block:: Python transformation = distX.getIsoProbabilisticTransformation() inverseTransformation = distX.getInverseIsoProbabilisticTransformation() .. GENERATED FROM PYTHON SOURCE LINES 157-161 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 163-168 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 170-172 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 172-174 .. code-block:: Python xi = vectorX.getRealization() .. GENERATED FROM PYTHON SOURCE LINES 175-179 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 179-181 .. code-block:: Python ui = [distX1.computeCDF(xi[0]), xi[1]] .. GENERATED FROM PYTHON SOURCE LINES 182-184 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 184-186 .. code-block:: Python zi = [-ot.Normal().computeInverseSurvivalFunction(ui[0])[0], ui[1]] .. GENERATED FROM PYTHON SOURCE LINES 187-188 The sought transform then maps a point in the physical space to the standard space : .. GENERATED FROM PYTHON SOURCE LINES 188-191 .. code-block:: Python print(xi, "->", ui, "->", zi) .. rst-class:: sphx-glr-script-out .. code-block:: none [0.0423321,1.51134] -> [0.04144857740421459, 1.5113438501790912] -> [-1.7341179636821182, 1.5113438501790912] .. GENERATED FROM PYTHON SOURCE LINES 192-194 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 194-197 .. code-block:: Python transformX1 = distX1.getIsoProbabilisticTransformation() inverseTransformX1 = distX1.getInverseIsoProbabilisticTransformation() .. GENERATED FROM PYTHON SOURCE LINES 198-204 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 204-211 .. code-block:: Python 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 .. code-block:: none zi = [-1.7341179636821182, 1.5113438501790912] zi1D = [-1.7341179636821182, 1.5113438501790912] zi2D = [-1.73412,1.51134] .. GENERATED FROM PYTHON SOURCE LINES 212-215 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 215-231 .. code-block:: Python 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, r"Boundary of the event $\partial \mathcal{D}$") curveCX.setLineStyle("solid") curveCX.setColor("blue") graphStandardSpace.add(curveCX) .. GENERATED FROM PYTHON SOURCE LINES 232-233 We add the origin to the previous graph. .. GENERATED FROM PYTHON SOURCE LINES 233-248 .. code-block:: Python cloud = ot.Cloud([0.0], [0.0]) cloud.setColor("black") cloud.setPointStyle("fcircle") cloud.setLegend("origin") graphStandardSpace.add(cloud) graphStandardSpace.setGrid(True) graphStandardSpace.setLegendPosition("lower right") # 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 249-257 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 259-260 We configure the Cobyla solver that we use for the optimization : .. GENERATED FROM PYTHON SOURCE LINES 260-267 .. code-block:: Python 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 268-270 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 270-272 .. code-block:: Python algoFORM = ot.FORM(solver, event, distX.getMean()) .. GENERATED FROM PYTHON SOURCE LINES 273-274 We are ready to run the algorithm and store the result : .. GENERATED FROM PYTHON SOURCE LINES 274-277 .. code-block:: Python algoFORM.run() result = algoFORM.getResult() .. GENERATED FROM PYTHON SOURCE LINES 278-280 The design point can be retrieved in both physical and standard space with respectively the `getPhysicalSpaceDesignPoint` and `getStandardSpaceDesignPoint` methods. .. GENERATED FROM PYTHON SOURCE LINES 280-286 .. code-block:: Python 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 .. 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 287-288 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 288-292 .. code-block:: Python betaHL = result.getHasoferReliabilityIndex() print("Hasofer index : ", betaHL) .. rst-class:: sphx-glr-script-out .. code-block:: none Hasofer index : 3.176696193499833 .. GENERATED FROM PYTHON SOURCE LINES 293-294 We visualize it on the previous graph. .. GENERATED FROM PYTHON SOURCE LINES 294-311 .. code-block:: Python 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("lower right") 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 312-323 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 325-332 .. code-block:: Python 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 .. code-block:: none abscissa of the design point u0 = 2.413844253971037 value of the failure boundary at u0 = 2.065335164471684 value of the gradient of the failure boundary at u0 = -1.1706609709179403 .. GENERATED FROM PYTHON SOURCE LINES 333-339 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 339-347 .. code-block:: Python 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 348-356 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 356-359 .. code-block:: Python pf = ot.Normal().computeSurvivalFunction(betaHL) print("FORM : Pf = ", pf) .. rst-class:: sphx-glr-script-out .. code-block:: none FORM : Pf = 0.0007448149708283228 .. GENERATED FROM PYTHON SOURCE LINES 360-361 This proability of failure is the one computed in the FORMResult and obtained with the `getEventProbability` method : .. GENERATED FROM PYTHON SOURCE LINES 361-365 .. code-block:: Python pf = result.getEventProbability() print("Probability of failure (FORM) Pf = ", pf) .. rst-class:: sphx-glr-script-out .. code-block:: none Probability of failure (FORM) Pf = 0.0007448149708283228 .. GENERATED FROM PYTHON SOURCE LINES 366-374 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 376-378 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 378-384 .. code-block:: Python 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 .. code-block:: none abscissa of the design point u0 = 2.413844253971037 value of the hessian of the failure boundary at u0 = 0.9401058369722994 .. GENERATED FROM PYTHON SOURCE LINES 385-391 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 391-404 .. code-block:: Python 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 405-415 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 oscillating parabola of concern we use the gradient and hessian previously computed : .. GENERATED FROM PYTHON SOURCE LINES 415-419 .. code-block:: Python curvature = (d2u0[0, 0, 0]) / (1 + (du0[0, 0]) ** 2) ** (3 / 2) print("Curvature (analytic formula) = ", curvature) .. rst-class:: sphx-glr-script-out .. code-block:: none Curvature (analytic formula) = 0.2575913913877353 .. GENERATED FROM PYTHON SOURCE LINES 420-421 We build the SORM algorithm and run it : .. GENERATED FROM PYTHON SOURCE LINES 421-424 .. code-block:: Python algoSORM = ot.SORM(solver, event, distX.getMean()) algoSORM.run() .. GENERATED FROM PYTHON SOURCE LINES 425-426 The SORM result is obtained with the `getResult` method : .. GENERATED FROM PYTHON SOURCE LINES 426-428 .. code-block:: Python resultSORM = algoSORM.getResult() .. GENERATED FROM PYTHON SOURCE LINES 429-431 The principal curvatures of the osculating paraboloid at the design point is obtained by the `getSortedCurvatures` method : .. GENERATED FROM PYTHON SOURCE LINES 431-434 .. code-block:: Python print("Curvature (estimated) = ", resultSORM.getSortedCurvatures()[1]) .. rst-class:: sphx-glr-script-out .. code-block:: none Curvature (estimated) = 0.25761034541451805 .. GENERATED FROM PYTHON SOURCE LINES 435-444 Once the curvature is obtained there are several ways of approximating the failure probability :math:`P_f`. The library 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 444-448 .. code-block:: Python 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 .. code-block:: none SORM : Pf = 0.0005523531956154747 .. GENERATED FROM PYTHON SOURCE LINES 449-450 We can compare with the different estimators implemented in OpenTURNS : .. GENERATED FROM PYTHON SOURCE LINES 450-458 .. code-block:: Python 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 .. code-block:: none Probability of failure (SORM Breintung) Pf = 0.0005523440504786004 Probability of failure (SORM Hohenbichler) Pf = 0.0005420328660300074 Probability of failure (SORM Tvedt) Pf = 0.0005381057564255441 .. GENERATED FROM PYTHON SOURCE LINES 459-460 Display all figures .. GENERATED FROM PYTHON SOURCE LINES 460-462 .. code-block:: Python otv.View.ShowAll() .. GENERATED FROM PYTHON SOURCE LINES 463-464 Reset default settings .. GENERATED FROM PYTHON SOURCE LINES 464-465 .. code-block:: Python ot.ResourceMap.Reload() .. _sphx_glr_download_auto_reliability_sensitivity_reliability_plot_form_explained.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_form_explained.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_form_explained.py `