1""" 2Use the FORM - SORM algorithms 3============================== 4""" 5# %% 6# In this example we estimate a failure probability with the `FORM` algorithm on the :ref:`cantilever beam <use-case-cantilever-beam>` example. More precisely, we show how to use the associated results: 7# 8# - the design point in both physical and standard space, 9# - the probability estimation according to the FORM approximation, and the following SORM ones: Tvedt, Hohenbichler and Breitung, 10# - the Hasofer reliability index and the generalized ones evaluated from the Breitung, Tvedt and Hohenbichler approximations, 11# - the importance factors defined as the normalized director factors of the design point in the :math:`U`-space 12# - the sensitivity factors of the Hasofer reliability index and the FORM probability. 13# - the coordinates of the mean point in the standard event space. 14# 15# The coordinates of the mean point in the standard event space is: 16# 17# .. math:: 18# \frac{1}{E_1(-\beta)}\int_{\beta}^{\infty} u_1 p_1(u_1)du_1 19# 20# 21# where :math:`E_1` is the spheric univariate distribution of the standard space and :math:`\beta` is the reliability index. 22 23# %% 24# Model definition 25# ---------------- 26 27# %% 28from __future__ import print_function 29from openturns.usecases import cantilever_beam as cantilever_beam 30import openturns as ot 31import openturns.viewer as viewer 32from matplotlib import pylab as plt 33ot.Log.Show(ot.Log.NONE) 34 35# %% 36# We load the model from the usecases module : 37cb = cantilever_beam.CantileverBeam() 38 39# %% 40# We use the input parameters distribution from the data class : 41distribution = cb.distribution 42distribution.setDescription(['E', 'F', 'L', 'I']) 43 44# %% 45# We define the model 46model = cb.model 47 48# %% 49# Create the event whose probability we want to estimate. 50 51# %% 52vect = ot.RandomVector(distribution) 53G = ot.CompositeRandomVector(model, vect) 54event = ot.ThresholdEvent(G, ot.Greater(), 0.3) 55event.setName("deviation") 56 57# %% 58# FORM Analysis 59# ------------- 60 61# %% 62# Define a solver 63optimAlgo = ot.Cobyla() 64optimAlgo.setMaximumEvaluationNumber(1000) 65optimAlgo.setMaximumAbsoluteError(1.0e-10) 66optimAlgo.setMaximumRelativeError(1.0e-10) 67optimAlgo.setMaximumResidualError(1.0e-10) 68optimAlgo.setMaximumConstraintError(1.0e-10) 69 70# %% 71# Run FORM 72algo = ot.FORM(optimAlgo, event, distribution.getMean()) 73algo.run() 74result = algo.getResult() 75 76# %% 77# Analysis of the results 78# ----------------------- 79 80# %% 81# Probability 82result.getEventProbability() 83 84# %% 85# Hasofer reliability index 86result.getHasoferReliabilityIndex() 87 88# %% 89# Design point in the standard U* space. 90 91# %% 92print(result.getStandardSpaceDesignPoint()) 93 94# %% 95# Design point in the physical X space. 96 97# %% 98print(result.getPhysicalSpaceDesignPoint()) 99 100# %% 101# Importance factors 102graph = result.drawImportanceFactors() 103view = viewer.View(graph) 104 105# %% 106marginalSensitivity, otherSensitivity = result.drawHasoferReliabilityIndexSensitivity() 107marginalSensitivity.setLegends(["E", "F", "L", "I"]) 108marginalSensitivity.setLegendPosition('bottom') 109view = viewer.View(marginalSensitivity) 110 111# %% 112marginalSensitivity, otherSensitivity = result.drawEventProbabilitySensitivity() 113marginalSensitivity.setLegends(["E", "F", "L", "I"]) 114marginalSensitivity.setLegendPosition('bottom') 115view = viewer.View(marginalSensitivity) 116 117# %% 118# Error history 119optimResult = result.getOptimizationResult() 120graphErrors = optimResult.drawErrorHistory() 121graphErrors.setLegendPosition('bottom') 122graphErrors.setYMargin(0.0) 123view = viewer.View(graphErrors) 124 125# %% 126# Get additional results with SORM 127algo = ot.SORM(optimAlgo, event, distribution.getMean()) 128algo.run() 129sorm_result = algo.getResult() 130 131# %% 132# Reliability index with Breitung approximation 133sorm_result.getGeneralisedReliabilityIndexBreitung() 134 135# %% 136# ... with Hohenbichler approximation 137sorm_result.getGeneralisedReliabilityIndexHohenbichler() 138 139# %% 140# .. with Tvedt approximation 141sorm_result.getGeneralisedReliabilityIndexTvedt() 142 143# %% 144# SORM probability of the event with Breitung approximation 145sorm_result.getEventProbabilityBreitung() 146 147# %% 148# ... with Hohenbichler approximation 149sorm_result.getEventProbabilityHohenbichler() 150 151# %% 152# ... with Tvedt approximation 153sorm_result.getEventProbabilityTvedt() 154 155plt.show() 156 157# %% 158# FORM analysis with finite difference gradient 159# --------------------------------------------- 160 161# %% 162# When the considered function has no analytical expression, the gradient may not be known. 163# In this case, a constant step finite difference gradient definition may be used. 164 165# %% 166 167 168def cantilever_beam_python(X): 169 E, F, L, I = X 170 return [F*L**3/(3*E*I)] 171 172 173cbPythonFunction = ot.PythonFunction(4, 1, func=cantilever_beam_python) 174epsilon = [1e-8]*4 # Here, a constant step of 1e-8 is used for every dimension 175gradStep = ot.ConstantStep(epsilon) 176cbPythonFunction.setGradient(ot.CenteredFiniteDifferenceGradient(gradStep, 177 cbPythonFunction.getEvaluation())) 178G = ot.CompositeRandomVector(cbPythonFunction, vect) 179event = ot.ThresholdEvent(G, ot.Greater(), 0.3) 180event.setName("deviation") 181 182# %% 183# However, given the different nature of the model variables, a blended (variable) 184# finite difference step may be preferable: 185# The step depends on the location in the input space 186gradStep = ot.BlendedStep(epsilon) 187cbPythonFunction.setGradient(ot.CenteredFiniteDifferenceGradient(gradStep, 188 cbPythonFunction.getEvaluation())) 189G = ot.CompositeRandomVector(cbPythonFunction, vect) 190event = ot.ThresholdEvent(G, ot.Greater(), 0.3) 191event.setName("deviation") 192 193# %% 194# We can then run the FORM analysis in the same way as before: 195algo = ot.FORM(optimAlgo, event, distribution.getMean()) 196algo.run() 197result = algo.getResult() 198