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
|
#!/usr/bin/env python
import openturns as ot
#
# Physical model
#
formulas = [
"min(0.1 * (u1 - u2)^2.0 - (u1 + u2) / sqrt(2.0) + 3.0, 0.1 * (u1 - u2)^2.0 + (u1 + u2) / sqrt(2.0) + 3.0, u1 - u2 + 3.5 * sqrt(2.0), -u1 + u2 + 3.5 * sqrt(2.0))"
]
limitState = ot.SymbolicFunction(["u1", "u2"], formulas)
dim = limitState.getInputDimension()
#
# Probabilistic model
#
mean = ot.Point(dim, 0.0)
sigma = ot.Point(dim, 1.0)
R = ot.IdentityMatrix(dim)
myDistribution = ot.Normal(mean, sigma, R)
#
# Limit state
#
vect = ot.RandomVector(myDistribution)
output = ot.CompositeRandomVector(limitState, vect)
myEvent = ot.ThresholdEvent(output, ot.Less(), 0.0)
#
# Computation
#
bs = 1
# Monte Carlo
experiment = ot.MonteCarloExperiment()
myMC = ot.ProbabilitySimulationAlgorithm(myEvent, experiment)
myMC.setMaximumOuterSampling(int(1e6) // bs)
myMC.setBlockSize(bs)
myMC.setMaximumCoefficientOfVariation(-1.0)
myMC.run()
#
# SubsetSampling
mySS = ot.SubsetSampling(myEvent)
mySS.setMaximumOuterSampling(10000 // bs)
mySS.setBlockSize(bs)
mySS.run()
#
# Results
#
# Monte Carlo
resultMC = myMC.getResult()
PFMC = resultMC.getProbabilityEstimate()
CVMC = resultMC.getCoefficientOfVariation()
variance_PF_MC = resultMC.getVarianceEstimate()
length90MC = resultMC.getConfidenceLength(0.90)
N_MC = resultMC.getOuterSampling() * resultMC.getBlockSize()
#
# SubsetSampling
resultSS = mySS.getResult()
PFSS = resultSS.getProbabilityEstimate()
CVSS = resultSS.getCoefficientOfVariation()
variance_PF_SS = resultSS.getVarianceEstimate()
length90SS = resultSS.getConfidenceLength(0.90)
N_SS = resultSS.getOuterSampling() * resultSS.getBlockSize()
#
print("")
print(
"************************************************************************************************"
)
print(
"**************************************** MONTE CARLO *******************************************"
)
print(
"************************************************************************************************"
)
print("Pf estimation = %.5e" % PFMC)
print("Pf Variance estimation = %.5e" % variance_PF_MC)
print("CoV = %.5f" % CVMC)
print("90% Confidence Interval =", "%.5e" % length90MC)
print(
"CI at 90% =[",
"%.5e" % (PFMC - 0.5 * length90MC),
"; %.5e" % (PFMC + 0.5 * length90MC),
"]",
)
print("Limit state calls =", N_MC)
print(
"************************************************************************************************"
)
print("")
print(
"************************************************************************************************"
)
print(
"******************************************* SUBSET SAMPLING **********************************************"
)
print(
"************************************************************************************************"
)
print("Pf estimation = %.5e" % PFSS)
print("Pf Variance estimation = %.5e" % variance_PF_SS)
print("CoV = %.5f" % CVSS)
print("90% Confidence Interval =", "%.5e" % length90SS)
print(
"CI at 90% =[",
"%.5e" % (PFSS - 0.5 * length90SS),
"; %.5e" % (PFSS + 0.5 * length90SS),
"]",
)
print("Limit state calls =", N_SS)
print(
"************************************************************************************************"
)
print("")
|
