1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
|
"""
Estimate a buckling probability
===============================
"""
# %%
#
# In this example, we estimate the probability that the output of a function
# exceeds a given threshold with the FORM method, the SORM method and an advanced
# sampling method.
#
# We consider the :ref:`stiffened panel model <use-case-stiffened-panel>`.
# %%
# Define the model
# ----------------
# %%
from openturns.usecases import stiffened_panel
import openturns as ot
import openturns.viewer as otv
# %%
# We load the stiffened panel model from the usecases module :
panel = stiffened_panel.StiffenedPanel()
distribution = panel.distribution
model = panel.model
# %%
# See the input distribution
distribution
# %%
# See the model
model.getOutputDescription()
# %%
# Draw the distribution of a sample of the output.
sampleSize = 1000
inputSample = distribution.getSample(sampleSize)
outputSample = model(inputSample)
graph = ot.HistogramFactory().build(outputSample).drawPDF()
_ = otv.View(graph)
# %%
# Define the event
# ----------------
# %%
# Then we create the event whose probability we want to estimate.
# %%
vect = ot.RandomVector(distribution)
criticalLoad = ot.CompositeRandomVector(model, vect)
minimumCriticalLoad = 165.0
event = ot.ThresholdEvent(criticalLoad, ot.Less(), minimumCriticalLoad)
event.setName("buckling")
# %%
# Estimate the probability with FORM
# ----------------------------------
# %%
# Define a solver.
# %%
optimAlgo = ot.Cobyla()
optimAlgo.setMaximumCallsNumber(1000)
optimAlgo.setMaximumAbsoluteError(1.0e-10)
optimAlgo.setMaximumRelativeError(1.0e-10)
optimAlgo.setMaximumResidualError(1.0e-10)
optimAlgo.setMaximumConstraintError(1.0e-10)
# %%
# Run FORM.
# %%
optimAlgo.setStartingPoint(distribution.getMean())
algo = ot.FORM(optimAlgo, event)
n0 = model.getCallsNumber()
algo.run()
n1 = model.getCallsNumber()
result = algo.getResult()
standardSpaceDesignPoint = result.getStandardSpaceDesignPoint()
# %%
# Retrieve results.
# %%
result = algo.getResult()
probability = result.getEventProbability()
print("Pf (FORM)=%.3e" % probability, "nb evals=", n1 - n0)
# %%
# Importance factors.
# %%
graph = result.drawImportanceFactors()
view = otv.View(graph)
# %%
# Estimate the probability with SORM
# ----------------------------------
# %%
# Run SORM.
# %%
algo = ot.SORM(optimAlgo, event)
n0 = model.getCallsNumber()
algo.run()
n1 = model.getCallsNumber()
# %%
# Retrieve results.
# %%
result = algo.getResult()
probability = result.getEventProbabilityBreitung()
print("Pf (SORM)=%.3e" % probability, "nb evals=", n1 - n0)
# %%
# We see that the FORM and SORM approximations give significantly different
# results. Use a simulation algorithm to get a confidence interval.
# %%
# Estimate the probability with PostAnalyticalControlledImportanceSampling
# ------------------------------------------------------------------------
# %%
algo = ot.PostAnalyticalControlledImportanceSampling(result)
algo.setBlockSize(100)
algo.setMaximumOuterSampling(100)
algo.setMaximumCoefficientOfVariation(0.1)
n0 = model.getCallsNumber()
algo.run()
n1 = model.getCallsNumber()
result = algo.getResult()
Pf = result.getProbabilityEstimate()
print("Pf (sim) = %.3e" % Pf, "nb evals=", n1 - n0)
width = result.getConfidenceLength(0.95)
print("C.I (95%)=[" + "%.3e" % (Pf - 0.5 * width), ",%.3e" % (Pf + 0.5 * width), "]")
# %%
# Display all figures
otv.View.ShowAll()
|
