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