1#! /usr/bin/env python
2
3from __future__ import print_function
4from openturns import *
5# from math import *
6
7TESTPREAMBLE()
8RandomGenerator.SetSeed(0)
9
10try:
11
12    #
13    # Physical model
14    #
15
16    inputFunction = Description(2)
17    inputFunction[0] = "X1"
18    inputFunction[1] = "X2"
19
20    outputFunction = Description(1)
21    outputFunction[0] = "G"
22
23    formulas = Description(outputFunction.getSize())
24    formulas[0] = "2.5 + 0.1*(X1-X2)^2.0 - (X1+X2)/sqrt(2.0)"
25
26    EtatLimite = Function(inputFunction, outputFunction, formulas)
27
28    dim = EtatLimite.getInputDimension()
29    print(dim)
30
31    #
32    # Probabilistic model
33    #
34
35    mean = Point(dim, 0.0)
36    mean[0] = 0.0
37    mean[1] = 0.0
38
39    sigma = Point(dim, 0.0)
40    sigma[0] = 1.0
41    sigma[1] = 1.0
42
43    R = IdentityMatrix(dim)
44
45    myDistribution = Normal(mean, sigma, R)
46
47    start = myDistribution.getMean()
48    Covariance = myDistribution.getCovariance()
49
50    #
51    # Limit state
52    #
53
54    vect = RandomVector(myDistribution)
55
56    output = RandomVector(EtatLimite, vect)
57
58    myEvent = ThresholdEvent(output, Less(), 0.0)
59
60    #
61    # Calculs
62    #
63    #
64    # FORM/SORM Cobyla
65    myCobyla = Cobyla()
66    myCobyla.setMaximumEvaluationNumber(100 * dim)
67    myCobyla.setMaximumAbsoluteError(1.0e-10)
68    myCobyla.setMaximumRelativeError(1.0e-10)
69    myCobyla.setMaximumResidualError(1.0e-10)
70    myCobyla.setMaximumConstraintError(1.0e-10)
71
72    myAlgoC = FORM(myCobyla, myEvent, start)
73    myAlgoC2 = SORM(myCobyla, myEvent, start)
74
75    myAlgoC.run()
76    myAlgoC2.run()
77
78    resultC = FORMResult(myAlgoC.getResult())
79    resultC2 = SORMResult(myAlgoC2.getResult())
80
81    #
82    # FORM/SORM Abdo Rackwitz
83    myAbdoRackwitz = AbdoRackwitz()
84    myAbdoRackwitz.setMaximumIterationNumber(100 * dim)
85    myAbdoRackwitz.setMaximumAbsoluteError(1.0e-10)
86    myAbdoRackwitz.setMaximumRelativeError(1.0e-10)
87    myAbdoRackwitz.setMaximumResidualError(1.0e-10)
88    myAbdoRackwitz.setMaximumConstraintError(1.0e-10)
89
90    myAlgoAR = FORM(myAbdoRackwitz, myEvent, start)
91    myAlgoAR2 = SORM(myAbdoRackwitz, myEvent, start)
92
93    myAlgoAR.run()
94    myAlgoAR2.run()
95
96    resultAR = FORMResult(myAlgoAR.getResult())
97    resultAR2 = SORMResult(myAlgoAR2.getResult())
98
99    #
100    # Monte Carlo
101    CoV_MC = 0.5
102    myMC = MonteCarlo(myEvent)
103    myMC.setMaximumOuterSampling(100000)
104    myMC.setBlockSize(1)
105    myMC.setMaximumCoefficientOfVariation(CoV_MC)
106    myMC.run()
107
108    #
109    # LHS
110    CoV_LHS = 0.1
111    myLHS = LHS(myEvent)
112    myLHS.setMaximumOuterSampling(100000)
113    myLHS.setBlockSize(1)
114    myLHS.setMaximumCoefficientOfVariation(CoV_LHS)
115    myLHS.run()
116
117    #
118
119    #
120    # Outputs
121    #
122    print("")
123    print("")
124    print(
125        "************************************************************************************************")
126    print(
127        "***************************************** FORM  COBYLA *****************************************")
128    print(
129        "************************************************************************************************")
130    print("event probability = %.5e" % PfC)
131    print("generalized reliability index = %.5f" % Beta_generalizedC)
132    print(
133        "************************************************************************************************")
134    for i in range(u_starC.getDimension()):
135        print("standard space design point = %.5f" % u_starC[i])
136    print(
137        "************************************************************************************************")
138    for i in range(x_starC.getDimension()):
139        print("physical space design point = %.5f" % x_starC[i])
140    print(
141        "************************************************************************************************")
142    print("is standard point origin in failure space? ", PtC)
143    print(
144        "************************************************************************************************")
145    for i in range(gammaC.getDimension()):
146        print("importance factors = %.5f" % gammaC[i])
147    print(
148        "************************************************************************************************")
149    print("Hasofer reliability index = %.5f" % beta_hasoferC)
150    print(
151        "************************************************************************************************")
152    for i in range(SensitivityC.getSize()):
153        for j in range(SensitivityC[i].getDimension()):
154            print("Pf sensitivity = %.5f" % SensitivityC[i][j])
155    print(
156        "************************************************************************************************")
157    print("")
158    print(
159        "************************************************************************************************")
160    print(
161        "************************************** FORM ABDO RACKWITZ **************************************")
162    print(
163        "************************************************************************************************")
164    print("event probability = %.5e" % PfAR)
165    print("generalized reliability index = %.5f" % Beta_generalizedAR)
166    print(
167        "************************************************************************************************")
168    for i in range(u_starAR.getDimension()):
169        print("standard space design point = %.5f" % u_starAR[i])
170    print(
171        "************************************************************************************************")
172    for i in range(x_starAR.getDimension()):
173        print("physical space design point = %.5f" % x_starAR[i])
174    print(
175        "************************************************************************************************")
176    print("is standard point origin in failure space? ", PtAR)
177    print(
178        "************************************************************************************************")
179    for i in range(gammaAR.getDimension()):
180        print("importance factors = %.5f" % gammaAR[i])
181    print(
182        "************************************************************************************************")
183    print("Hasofer reliability index = %.5f" % beta_hasoferAR)
184    print(
185        "************************************************************************************************")
186    for i in range(SensitivityAR.getSize()):
187        for j in range(SensitivityAR[i].getDimension()):
188            print("Pf sensitivity = %.5f" % SensitivityAR[i][j])
189    print(
190        "************************************************************************************************")
191    print("")
192    print(
193        "************************************************************************************************")
194    print(
195        "***************************************** SORM  COBYLA *****************************************")
196    print(
197        "************************************************************************************************")
198    print("Breitung event probability = %.5e" % PFBreitC2)
199    print("Breitung generalized reliability index = %.5f" % BetaBreitC2)
200    print("Hohenbichler event probability = %.5e" % PFHBC2)
201    print("Hohenbichler generalized reliability index = %.5f" % BetaHBC2)
202    print("Tvedt event probability = %.5e" % PFTvedtC2)
203    print("Tvedt generalized reliability index = %.5f" % BetaTvedtC2)
204    print(
205        "************************************************************************************************")
206    for i in range(CurvC2.getDimension()):
207        print("sorted curvatures = %.5f" % cleanScalar(CurvC2[i]))
208    print(
209        "************************************************************************************************")
210    for i in range(u_starC2.getDimension()):
211        print("standard space design point = %.5f" % u_starC2[i])
212    print(
213        "************************************************************************************************")
214    for i in range(x_starC2.getDimension()):
215        print("physical space design point = %.5f" % x_starC2[i])
216    print(
217        "************************************************************************************************")
218    print(
219        "************************************************************************************************")
220    print("is standard point origin in failure space? ", PtC2)
221    print(
222        "************************************************************************************************")
223    for i in range(gammaC2.getDimension()):
224        print("importance factors = %.5f" % gammaC2[i])
225    print(
226        "************************************************************************************************")
227    print("Hasofer reliability index = %.5f" % beta_hasoferC2)
228    print(
229        "************************************************************************************************")
230    print("")
231    print(
232        "************************************************************************************************")
233    print(
234        "************************************** SORM ABDO RACKWITZ **************************************")
235    print(
236        "************************************************************************************************")
237    print("Breitung event probability = %.5e" % PFBreitAR2)
238    print("Breitung generalized reliability index = %.5f" % BetaBreitAR2)
239    print("Hohenbichler event probability = %.5e" % PFHBAR2)
240    print("Hohenbichler generalized reliability index = %.5f" % BetaHBAR2)
241    print("Tvedt event probability = %.5e" % PFTvedtAR2)
242    print("Tvedt generalized reliability index = %.5f" % BetaTvedtAR2)
243    print(
244        "************************************************************************************************")
245    for i in range(CurvAR2.getDimension()):
246        print("sorted curvatures = %.5f" % cleanScalar(CurvAR2[i]))
247    print(
248        "************************************************************************************************")
249    for i in range(u_starAR2.getDimension()):
250        print("standard space design point = %.5f" % u_starAR2[i])
251    print(
252        "************************************************************************************************")
253    for i in range(x_starAR2.getDimension()):
254        print("physical space design point = %.5f" % x_starAR2[i])
255    print(
256        "************************************************************************************************")
257    print(
258        "************************************************************************************************")
259    print("is standard point origin in failure space? ", PtAR2)
260    print(
261        "************************************************************************************************")
262    for i in range(gammaAR2.getDimension()):
263        print("importance factors = %.5f" % gammaAR2[i])
264    print(
265        "************************************************************************************************")
266    print("Hasofer reliability index = %.5f" % beta_hasoferAR2)
267    print(
268        "************************************************************************************************")
269    print("")
270    print(
271        "************************************************************************************************")
272    print(
273        "**************************************** MONTE CARLO *******************************************")
274    print(
275        "************************************************************************************************")
276    print("Pf estimation = %.5e" % PFMC)
277    print("Pf Variance estimation = %.5e" % Variance_PF_MC)
278    print("CoV = %.5f" % CVMC)
279    print("90% Confidence Interval =", "%.5e" % length90MC)
280    print("CI at 90% =[", "%.5e" % (PFMC - 0.5 * length90MC),
281          "; %.5e" % (PFMC + 0.5 * length90MC), "]")
282    print(
283        "************************************************************************************************")
284    print("")
285    print(
286        "************************************************************************************************")
287    print(
288        "******************************************* L H S **********************************************")
289    print(
290        "************************************************************************************************")
291    print("Pf estimation = %.5e" % PFLHS)
292    print("Pf Variance estimation = %.5e" % Variance_PF_LHS)
293    print("CoV = %.5f" % CVLHS)
294    print("90% Confidence Interval =", "%.5e" % length90LHS)
295    print("CI at 90% =[", "%.5e" % (PFLHS - 0.5 * length90LHS),
296          "; %.5e" % (PFLHS + 0.5 * length90LHS), "]")
297    print(
298        "************************************************************************************************")
299    print("")
300
301except:
302    import sys
303    print("t_Waarts_convex.py", sys.exc_info()[0], sys.exc_info()[1])
304