.. 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-15 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. See :ref:`FORM ` and :ref:`SORM ` and to get more theoretical details. .. GENERATED FROM PYTHON SOURCE LINES 15-19 .. code-block:: Python import openturns as ot import openturns.viewer as otv import numpy as np .. GENERATED FROM PYTHON SOURCE LINES 20-29 Context ------- We consider a bivariate random vector :math:`\inputRV = (X_1, X_2)` with the following independent components that follow: - the exponential distribution with parameter :math:`\lambda=1`, :math:`X_1 \sim \mathcal{E}(1.0)` ; - the standard unit normal distribution :math:`X_2 \sim \mathcal{N}(0,1)`. The support of the input vector is :math:`[0, +\infty[ \times \Rset` .. GENERATED FROM PYTHON SOURCE LINES 29-34 .. code-block:: Python dist_X1 = ot.Exponential(1.0) dist_X2 = ot.Normal() dist_X = ot.JointDistribution([dist_X1, dist_X2]) .. GENERATED FROM PYTHON SOURCE LINES 35-36 We can draw the isolines of the PDF of the distribution `dist_X`: .. GENERATED FROM PYTHON SOURCE LINES 36-48 .. code-block:: Python ot.ResourceMap.SetAsUnsignedInteger("Contour-DefaultLevelsNumber", 8) graph_PDF = dist_X.drawPDF([0.0, -10.0], [20.0, 10.0]) graph_PDF.setTitle(r"2D-PDF of the input variables $(X_1, X_2)$") graph_PDF.setXTitle(r"$x_1$") graph_PDF.setYTitle(r"$x_2$") graph_PDF.setLegendPosition("lower right") contours = graph_PDF.getDrawable(0).getImplementation() contours.setColorMapNorm("log") graph_PDF.setDrawable(contours, 0) view = otv.View(graph_PDF, square_axes=True) .. 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 .. rst-class:: sphx-glr-script-out .. code-block:: none /home/devel/project/build/python/src/site-packages/openturns/viewer.py:729: UserWarning: Attempt to set non-positive ylim on a log-scaled axis will be ignored. colorbar = self._fig.colorbar( .. GENERATED FROM PYTHON SOURCE LINES 49-56 We consider the model from :math:`\Rset^2` into :math:`\Rset` defined by: .. math:: \model : (x_1, x_2) \mapsto x_1 x_2 We start by drawing the isolines of the model :math:`\model`. .. GENERATED FROM PYTHON SOURCE LINES 56-64 .. code-block:: Python g = ot.SymbolicFunction(["x1", "x2"], ["x1 * x2"]) graph_model = g.draw([0.0, -10.0], [20.0, 10.0]) graph_model.setXTitle(r"$x_1$") graph_model.setYTitle(r"$x_2$") graph_model.setTitle(r"Isolines of the model : $g$") view = otv.View(graph_model, square_axes=True) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_002.png :alt: Isolines of the model : $g$ :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 65-88 We consider the univariate output variable : .. math:: Y = \model(\inputRV) 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 for a continuous random vector :math:`\inputRV` as: .. math:: :label: PfDef P_f = \Prob{Y \geq s} = \int_{\set{D}} \mathbf{1}_{\set{D}}(x) \pdf d\vect{x} where: .. math:: \set{D} = \{ (x_1, x_2) \in [0,+\infty[ \times \Rset \, | \, \model(x_1, x_2) \geq s \} is the failure domain and :math:`\inputMeasure` is the probability density function (PDF) of :math:`\inputRV`. .. GENERATED FROM PYTHON SOURCE LINES 90-92 We first define random vectors and the failure event associated to the output random variable. .. GENERATED FROM PYTHON SOURCE LINES 92-98 .. code-block:: Python vector_X = ot.RandomVector(dist_X) vector_Y = ot.CompositeRandomVector(g, vector_X) s = 10.0 event = ot.ThresholdEvent(vector_Y, ot.Greater(), s) .. GENERATED FROM PYTHON SOURCE LINES 99-115 The boundary of the failure domain can easily be represented as it is a branch of an hyperbole: the boundary is the graph of the function defined from :math:`\Rset` into :math:`\Rset` by: .. math:: :label: defH h : x_1 \mapsto x_2 = \frac{s}{x_1} The boundary of the failure domain is also the isoline of the model :math:`\model` associated to the level :math:`s`: .. math:: \partial \set{D} = \{(x_1, x_2)\, |\, \model(x_1, x_2) = s \} We can draw it with the `draw` method of the function :math:`\model`. .. GENERATED FROM PYTHON SOURCE LINES 117-124 .. code-block:: Python nb_points = 101 graph_g = g.draw([0.0, -10.0], [20.0, 10.0], [nb_points] * 2) draw_boundary = graph_g.getDrawable(0) draw_boundary.setLevels([s]) draw_boundary.setLegend(r"Boundary $\partial \mathcal{D}$") graph_g.setDrawables([draw_boundary]) .. GENERATED FROM PYTHON SOURCE LINES 125-131 .. code-block:: Python texts = [r" $\mathcal{D} = \{(x_1, x_2)\, |\, g(x_1, x_2) \geq 10 \}$"] text_graph = ot.Text([[10.0, 3.0]], texts) text_graph.setTextSize(1) text_graph.setColor("black") graph_g.add(text_graph) .. GENERATED FROM PYTHON SOURCE LINES 132-139 .. code-block:: Python graph_g.setTitle("Failure domain in the physical space") graph_g.setXTitle(r"$x_1$") graph_g.setYTitle(r"$x_2$") graph_g.setLegendPosition("topright") view = otv.View(graph_g, square_axes=True) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_003.png :alt: Failure domain in the physical space :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 140-141 We can superimpose the event boundary with the bivariate PDF insolines of the input distribution: .. GENERATED FROM PYTHON SOURCE LINES 141-146 .. code-block:: Python draw_boundary.setColor("black") graph_PDF.add(draw_boundary) graph_PDF.setLegendPosition("lower right") view = otv.View(graph_PDF, square_axes=True) .. 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 147-152 From the previous figure, we observe that in the failure domain, the PDF takes small (and even very small) values. Consequently the failure probability :math:`P_f` is also expected to be small. The FORM/SORM methods estimate the failure probability. .. GENERATED FROM PYTHON SOURCE LINES 154-163 The FORM/SORM approximations ---------------------------- The basic steps of the FORM and SORM algorithms are: - use an isoprobabilistic transformation to map the input random vector into the standard space; - find the design point which is the nearest point to the origin in the standard space; - estimate the probability. .. GENERATED FROM PYTHON SOURCE LINES 165-186 Isoprobabilistic transformation ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The interest of the isoprobabilistic transformation is the rotational invariance of the distribution in the standard space. This property reduces the dimension of the problem to 1. See :ref:`Isoprobabilistic transformation ` for more theoretical details. OpenTURNS has several isoprobabilistic transformations, depending on the distribution of the input random vector: - the Nataf transformation is used if the input distribution has a normal copula, - the Generalized Nataf transformation is used if the input distribution has an elliptical copula, - the Rosenblatt transformation is used in any other cases. The Nataf and Rosenblatt transformations map the input random vector into a vector that follows a normal distribution with zero mean and unit variance. The Generalized Nataf transformation maps the input random vector into a vector that follows the standard elliptical distribution associated to the elliptical copula of the input distribution. In this example, the input distribution is not elliptical so the isoprobabilistic transformation is the Rosenblatt transformation. .. GENERATED FROM PYTHON SOURCE LINES 186-188 .. code-block:: Python print("Is Elliptical ? ", dist_X.isElliptical()) .. rst-class:: sphx-glr-script-out .. code-block:: none Is Elliptical ? False .. GENERATED FROM PYTHON SOURCE LINES 189-211 The Rosenblatt transformation :math:`T` is defined by: .. math:: :label: defT T : \vect{x} \mapsto \vect{z} such that the random vector :math:`\standardRV = T(\inputRV)` follows a bivariate normal distribution with zero mean and unit variance. It follows that the components :math:`Z_1` and :math:`Z_2` are independent. We detail the Rosenblatt transform in this simple case where the input random vector :math:`\inputRV` has independent components. Then, the Rosenblatt transform is defined by: .. math:: z_i = \Phi^{-1} \circ F_i(x_i) where :math:`F_i` is the cumulative distribution function (CDF) of :math:`X_i` and :math:`\Phi` the CDF of the univariate normal distribution with zero mean and unit variance. Note that in this example, :math:`\Phi^{-1} \circ F_2 = I_d` as :math:`F_2 = \Phi`. The isoprobabilistic transform and its inverse are methods of the distribution: .. GENERATED FROM PYTHON SOURCE LINES 211-214 .. code-block:: Python transformation = dist_X.getIsoProbabilisticTransformation() inverse_transformation = dist_X.getInverseIsoProbabilisticTransformation() .. GENERATED FROM PYTHON SOURCE LINES 215-217 Let us detail this transformation, step by step. We draw a realization of the random input vector. This point is said to be in the physical space. .. GENERATED FROM PYTHON SOURCE LINES 217-219 .. code-block:: Python xi = vector_X.getRealization() .. GENERATED FROM PYTHON SOURCE LINES 220-222 We build `zi` the mapping of `xi` through the Rosenblatt transformation. The point `zi` is said to be in the standard space. Note that the second component remained unchanged. .. GENERATED FROM PYTHON SOURCE LINES 222-226 .. code-block:: Python ui = [dist_X1.computeCDF(xi[0]), dist_X2.computeCDF(xi[1])] zi = [ot.Normal().computeQuantile(ui[0])[0], ot.Normal().computeQuantile(ui[1])[0]] print(xi, "->", ui, "->", zi) .. rst-class:: sphx-glr-script-out .. code-block:: none [0.38973,-0.540729] -> [0.32276015315467066, 0.2943473026427408] -> [-0.4599943077266245, -0.5407286830219163] .. GENERATED FROM PYTHON SOURCE LINES 227-235 We also build the isoprobabilistic transform :math:`T_1` and its inverse :math:`T_1^{-1}` for the first marginal: .. math:: :label: detT1 T_1 = \Phi^{-1} \circ F_1 .. GENERATED FROM PYTHON SOURCE LINES 235-238 .. code-block:: Python transform_X1 = dist_X1.getIsoProbabilisticTransformation() inverse_transform_X1 = dist_X1.getInverseIsoProbabilisticTransformation() .. GENERATED FROM PYTHON SOURCE LINES 239-241 We can implement the transformation using :math:`T_1` on the first components directly using :math:`T` on both components `xi`: .. GENERATED FROM PYTHON SOURCE LINES 241-244 .. code-block:: Python zi1D = [transform_X1([xi[0]])[0], xi[1]] zi2D = transformation(xi) .. GENERATED FROM PYTHON SOURCE LINES 245-246 We can check the result of our experiment : we observe the results are the same. .. GENERATED FROM PYTHON SOURCE LINES 246-251 .. code-block:: Python print("zi = ", zi) print("zi1D = ", zi1D) print("zi2D = ", zi2D) .. rst-class:: sphx-glr-script-out .. code-block:: none zi = [-0.4599943077266245, -0.5407286830219163] zi1D = [-0.4599943077266245, -0.5407286830219165] zi2D = [-0.459994,-0.540729] .. GENERATED FROM PYTHON SOURCE LINES 252-259 The model in the standard space is defined by: .. math:: \widetilde{\model} = \model \circ T^{-1} We can define it using the capacities of the composition of functions implemented in the library. .. GENERATED FROM PYTHON SOURCE LINES 259-261 .. code-block:: Python g_tilde = ot.ComposedFunction(g, inverse_transformation) .. GENERATED FROM PYTHON SOURCE LINES 262-275 The failure domain in the standard space is defined by: .. math:: \set{\widetilde{D}} = \{ (z_1, z_2) \in [0,+\infty[ \times \Rset \, | \, \widetilde{\model}(z_1, z_2) \geq s \} and its boundary is defined by: .. math:: \partial \set{\widetilde{D}} = \{ (z_1, z_2) \in [0,+\infty[ \times \Rset \, | \, \widetilde{\model}(z_1, z_2) = s \} .. GENERATED FROM PYTHON SOURCE LINES 277-278 We draw the graph of :math:`\widetilde{g}` in the standard space. .. GENERATED FROM PYTHON SOURCE LINES 278-289 .. code-block:: Python graph_standard_space = g_tilde.draw([0.0, 0.0], [7.0, 7.0], [101] * 2) draw_boundary_stand_space = graph_standard_space.getDrawable(0) draw_boundary_stand_space.setLevels([s]) draw_boundary_stand_space.setLegend(r"Boundary $\partial \mathcal{\tilde{D}}$") graph_standard_space.setDrawables([draw_boundary_stand_space]) graph_standard_space.setXTitle(r"$z_1$") graph_standard_space.setYTitle(r"$z_2$") graph_standard_space.setTitle("Failure domain in the standard space") .. GENERATED FROM PYTHON SOURCE LINES 290-291 We add some annotations .. GENERATED FROM PYTHON SOURCE LINES 291-300 .. code-block:: Python texts = [r"$\mathcal{\tilde{D}} = \{(z_1, z_2)\, |\, \tilde{g}(z_1, z_2) \geq 10 \}$"] text_graph = ot.Text([[4.0, 3.0]], texts) text_graph.setTextSize(1) text_graph.setColor("black") graph_standard_space.add(text_graph) graph_standard_space.setLegendPosition("topright") view = otv.View(graph_standard_space, square_axes=True) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_005.png :alt: Failure domain 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 301-309 The design point ^^^^^^^^^^^^^^^^ Due to the spherical distribution in the standard space, the area that contributes the most to the integral defining the probability is the vicinity of the closest point of the failure domain to the origin of the standard space. 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 311-312 We configure the Cobyla solver that we use for the optimization : .. GENERATED FROM PYTHON SOURCE LINES 312-319 .. 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 320-322 We build the FORM algorithm with its basic constructor. The starting point for the optimization algorithm is the mean of the input distribution. .. GENERATED FROM PYTHON SOURCE LINES 322-324 .. code-block:: Python algo_FORM = ot.FORM(solver, event, dist_X.getMean()) .. GENERATED FROM PYTHON SOURCE LINES 325-326 We are ready to run the algorithm and store the result. .. GENERATED FROM PYTHON SOURCE LINES 326-329 .. code-block:: Python algo_FORM.run() result = algo_FORM.getResult() .. GENERATED FROM PYTHON SOURCE LINES 330-333 The design point can be retrieved in both physical and standard space with respectively the `getPhysicalSpaceDesignPoint` and `getStandardSpaceDesignPoint` methods. We denote them respectively :math:`\vect{x}^*` and :math:`\vect{z}^*`. .. GENERATED FROM PYTHON SOURCE LINES 333-338 .. code-block:: Python design_point_physical_space = result.getPhysicalSpaceDesignPoint() design_point_standard_space = result.getStandardSpaceDesignPoint() print("Design point in physical space : ", design_point_physical_space) print("Design point in standard space : ", design_point_standard_space) .. 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 339-341 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 341-344 .. code-block:: Python beta_HL = result.getHasoferReliabilityIndex() print("Hasofer index : ", beta_HL) .. rst-class:: sphx-glr-script-out .. code-block:: none Hasofer index : 3.176696193499833 .. GENERATED FROM PYTHON SOURCE LINES 345-346 We visualize the design point on the previous graph. .. GENERATED FROM PYTHON SOURCE LINES 346-364 .. code-block:: Python cloud = ot.Cloud([design_point_standard_space]) cloud.setColor("red") cloud.setPointStyle("fcircle") cloud.setLegend(r"design point $z^*$") graph_standard_space.add(cloud) graph_standard_space.setGrid(True) graph_standard_space.setLegendPosition("lower right") cc = ot.Curve( [0.0, design_point_standard_space[0]], [0.0, design_point_standard_space[1]], r"$\beta_{HL}$ distance", ) cc.setLineStyle("dashed") cc.setColor("black") graph_standard_space.add(cc) graph_standard_space.setLegendPosition("topright") view = otv.View(graph_standard_space, square_axes=True) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_006.png :alt: Failure domain 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 365-388 The FORM approximation ^^^^^^^^^^^^^^^^^^^^^^ The last step of the FORM algorithm is to replace the failure domain boundary by the hyperplane which is tangent to the failure domain at the design point in the standard space. To draw this hyperplane :math:`\mathcal{P}_{\vect{z}^*}`, we define the function from :math:`\Rset^2` to :math:`\Rset` defined by: .. math:: M \rightarrow \scalarproduct{\nabla \widetilde{\model}(\vect{z}^*)}{\vect{Z^*M}} where :math:`\nabla \vect{\widetilde{\model}(\vect{z}^*)}` is the gradient of the function :math:`\widetilde{\model}` at the design point :math:`Z^*(\vect{z}^*)`. Then, the tangent hyperplane is the isoline associated to the zero level of the previous function : .. math:: \mathcal{P}_{z^*} = \{ \vect{z} \in \Rset^2 \, | \, \scalarproduct{\nabla\widetilde{\model}(\vect{z}^*)}{\vect{Z^*M}} = 0\} We can use the :class:`~openturns.LinearFunction` class. .. GENERATED FROM PYTHON SOURCE LINES 388-406 .. code-block:: Python center = design_point_standard_space grad_design_point = g_tilde.gradient(design_point_standard_space) constant = [0.0] linear_mat = ot.Matrix(1, 2) linear_mat[0, 0] = grad_design_point[0, 0] linear_mat[0, 1] = grad_design_point[1, 0] linear_proj = ot.LinearFunction(center, constant, linear_mat) graph_tangent = linear_proj.getMarginal(0).draw([0.0, 0.0], [7.0, 7.0], [101] * 2) draw_tangent = graph_tangent.getDrawable(0) draw_tangent.setLevels([0]) draw_tangent.setLegend(r"$\mathcal{\Pi}_{z^*}$ (FORM)") draw_tangent.setColor("green") draw_tangent.setLineStyle("dashed") graph_standard_space.add(draw_tangent) graph_standard_space.setLegendPosition("topright") view = otv.View(graph_standard_space, square_axes=True) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_007.png :alt: Failure domain 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 407-421 Depending on whether the origin of the standard space :math:`\vect{0}` belongs to the failure domain, the FORM probability is defined by: .. math:: P_{FORM} = \begin{cases} E(+\beta_{HL}) & \text{ if } \vect{0} \in \set{\widetilde{D}} \\ E(-\beta_{HL}) & \text{ otherwise.} \end{cases} where :math:`E(.)` is the marginal cumulative distribution function along any direction of the spherical distribution in the standard space. In this example, this is the normal distribution. So we have: .. GENERATED FROM PYTHON SOURCE LINES 421-429 .. code-block:: Python isOriginFail = result.getIsStandardPointOriginInFailureSpace() normal = ot.Normal() if isOriginFail: pf_FORM = normal.computeCDF(beta_HL) else: pf_FORM = normal.computeCDF(-beta_HL) print("FORM : Pf_FORM = ", pf_FORM) .. rst-class:: sphx-glr-script-out .. code-block:: none FORM : Pf_FORM = 0.0007448149708283228 .. GENERATED FROM PYTHON SOURCE LINES 430-432 This failure probability is implemented but the FORM algorithm and can be obtained with the `getEventProbability` method. We check we have the same result. .. GENERATED FROM PYTHON SOURCE LINES 432-435 .. code-block:: Python pf = result.getEventProbability() print("Probability of failure (FORM) Pf_FORM = ", pf) .. rst-class:: sphx-glr-script-out .. code-block:: none Probability of failure (FORM) Pf_FORM = 0.0007448149708283228 .. GENERATED FROM PYTHON SOURCE LINES 436-460 The SORM approximation ^^^^^^^^^^^^^^^^^^^^^^ The SORM approximation uses the main curvatures :math:`\kappa_i^0` of the homothetic of the failure domain at distance 1 from the origin. These curvatures are calculated at the design point. They are linked to the curvatures :math:`\kappa_i` of the failure domain by: .. math:: \kappa_i^0 = \beta_{HL} \kappa_i The Breitung approximation is valid for :math:`\beta_{HL} \rightarrow +\infty` and is defined by : .. math:: P_{SORM, Breitung} = \begin{cases} E(+\beta_{HL}) \prod_{i=1}^{d-1} \dfrac{1}{\sqrt{1+\kappa_i^0}} & \text{ if } \vect{0} \in \set{\widetilde{D}} \\ E(-\beta_{HL}) \prod_{i=1}^{d-1} \dfrac{1}{\sqrt{1+\kappa_i^0}} & \text{ otherwise. } \end{cases} and approximates the boundary by the osculating paraboloid at the design point. Note that the term :math:`\kappa_i^0` does not depend on :math:`\beta_{HL}`. .. GENERATED FROM PYTHON SOURCE LINES 462-467 In this example, we can easily implement the boundary of the failure domain in the physical space, using the function :math:`h` defined in :eq:`defH`. In the standard space, the boundary is defined by the composed function :math:`z_1 \mapsto h \circ T_1^{-1}(z_1)`. .. GENERATED FROM PYTHON SOURCE LINES 467-471 .. code-block:: Python failure_boundary_physical_space = ot.SymbolicFunction(["x"], ["10.0 / x"]) failure_boundary_standard_space = ot.ComposedFunction( failure_boundary_physical_space, inverse_transform_X1 ) .. GENERATED FROM PYTHON SOURCE LINES 472-474 We need the value of the second derivative of the failure boundary function at the abscissa of the design point in the standard space: .. GENERATED FROM PYTHON SOURCE LINES 474-488 .. code-block:: Python z1_star = [design_point_standard_space[0]] dz1_star = failure_boundary_standard_space.getGradient().gradient(z1_star) d2z1_star = failure_boundary_standard_space.getHessian().hessian(z1_star) print("first component of the design point = ", z1_star[0]) print( "second component of the design point = ", failure_boundary_standard_space(z1_star)[0], ) print( "value of the hessian of the failure boundary at this abscissa= ", d2z1_star[0, 0, 0], ) .. rst-class:: sphx-glr-script-out .. code-block:: none first component of the design point = 2.413844253971037 second component of the design point = 2.065335164471684 value of the hessian of the failure boundary at this abscissa= 0.9401058369722994 .. GENERATED FROM PYTHON SOURCE LINES 489-497 In the standard space, the osculating parabola :math:`\mathcal{P}_{\vect{z}^*}` at :math:`\vect{z}^*` is the graph of the function defined by: .. math:: z_1 \mapsto h \circ T_1^{-1} (z_1^*) + \frac{d}{du_1} (h \circ T_1^{-1})(z_1^*) (z_1-z_1^*) + \frac{1}{2} \frac{d^2}{dz_1^2} (h \circ T_1^{-1})(z_1^*) (z_1-z_1^*)^2 .. GENERATED FROM PYTHON SOURCE LINES 497-511 .. code-block:: Python z = np.linspace(1.1, 4.0, 100) parabola = ( failure_boundary_standard_space(z1_star)[0] + dz1_star[0, 0] * (z - z1_star) + 0.5 * d2z1_star[0, 0, 0] * (z - z1_star) ** 2 ) curve_parabola = ot.Curve(z, parabola, r"$\mathcal{P}_{z^*}$ (SORM)") curve_parabola.setLineStyle("dashed") curve_parabola.setColor("orange") graph_standard_space.add(curve_parabola) graph_standard_space.setLegendPosition("topright") view = otv.View(graph_standard_space) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_form_explained_008.png :alt: Failure domain 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 512-523 The next step is to estimate the principal curvatures of the osculating paraboloid. For any regular function :math:`\ell: \Rset \rightarrow \Rset` the curvature :math:`\kappa(x)` at the point :math:`x` in cartesian coordinates reads as: .. math:: \kappa(x) = \frac{\ell''(x)}{(1+[\ell'(x)]^2)^{3/2}}. For the osculating parabola of concern we use the previously computed gradient and Hessian previously computed: .. GENERATED FROM PYTHON SOURCE LINES 523-527 .. code-block:: Python curvature = (d2z1_star[0, 0, 0]) / (1 + (dz1_star[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 528-529 We build the SORM algorithm and run it : .. GENERATED FROM PYTHON SOURCE LINES 529-532 .. code-block:: Python algo_SORM = ot.SORM(solver, event, dist_X.getMean()) algo_SORM.run() .. GENERATED FROM PYTHON SOURCE LINES 533-534 The SORM result is obtained with the `getResult` method: .. GENERATED FROM PYTHON SOURCE LINES 534-536 .. code-block:: Python result_SORM = algo_SORM.getResult() .. GENERATED FROM PYTHON SOURCE LINES 537-539 The principal curvatures of the osculating paraboloid at the design point is obtained by the `getSortedCurvatures` method: .. GENERATED FROM PYTHON SOURCE LINES 539-542 .. code-block:: Python print("Curvature (library) = ", result_SORM.getSortedCurvatures()[1]) .. rst-class:: sphx-glr-script-out .. code-block:: none Curvature (library) = 0.25761034541451805 .. GENERATED FROM PYTHON SOURCE LINES 543-546 Once the curvature is computed, there are several ways of approximating the failure probability :math:`P_f`. The library implements the Breitung, Hohenbichler and Tvedt estimates. We detail here the calculus of the Breitung approximation. .. GENERATED FROM PYTHON SOURCE LINES 546-553 .. code-block:: Python coeff = (1.0 + beta_HL * curvature) ** (-0.5) if isOriginFail: pf_SORM = (normal.computeCDF(beta_HL)) * coeff else: pf_SORM = (normal.computeCDF(-beta_HL)) * coeff print("SORM : Pf_SORM = ", pf_SORM) .. rst-class:: sphx-glr-script-out .. code-block:: none SORM : Pf_SORM = 0.0005523531956154584 .. GENERATED FROM PYTHON SOURCE LINES 554-555 We can compare with the different estimators: .. GENERATED FROM PYTHON SOURCE LINES 555-563 .. code-block:: Python pf_Breitung = result_SORM.getEventProbabilityBreitung() pf_Hohenbichler = result_SORM.getEventProbabilityHohenbichler() pf_Tvedt = result_SORM.getEventProbabilityTvedt() print("Probability of failure (SORM Breintung) Pf = ", pf_Breitung) print("Probability of failure (SORM Hohenbichler) Pf = ", pf_Hohenbichler) print("Probability of failure (SORM Tvedt) Pf = ", pf_Tvedt) .. 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 564-565 Display all figures .. GENERATED FROM PYTHON SOURCE LINES 565-567 .. code-block:: Python otv.View.ShowAll() .. GENERATED FROM PYTHON SOURCE LINES 568-569 Reset default settings .. GENERATED FROM PYTHON SOURCE LINES 569-570 .. 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 ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_form_explained.zip `