1#! /usr/bin/env python
2
3from __future__ import print_function
4from openturns import *
5from math import *
6
7TESTPREAMBLE()
8
9
10def cleanScalar(inScalar):
11    if (fabs(inScalar) < 1.e-10):
12        inScalar = 0.0
13    return inScalar
14
15#
16# initialize the random generator
17#
18
19
20RandomGenerator.SetSeed(0)
21
22try:
23
24    #
25    # Physical model
26    #
27
28    inputFunction = Description(2)
29    inputFunction[0] = "u1"
30    inputFunction[1] = "u2"
31
32    formulas = Description(1)
33    formulas[
34        0] = "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))"
35
36    limitState = SymbolicFunction(inputFunction, formulas)
37    limitState.setGradient(CenteredFiniteDifferenceGradient(ResourceMap.GetAsScalar(
38        "CenteredFiniteDifferenceGradient-DefaultEpsilon"), limitState.getEvaluation()))
39    limitState.setHessian(CenteredFiniteDifferenceHessian(ResourceMap.GetAsScalar(
40        "CenteredFiniteDifferenceHessian-DefaultEpsilon"), limitState.getEvaluation()))
41    dim = limitState.getInputDimension()
42
43    #
44    # Probabilistic model
45    #
46
47    mean = Point(dim, 0.0)
48
49    sigma = Point(dim, 1.0)
50
51    R = IdentityMatrix(dim)
52    myDistribution = Normal(mean, sigma, R)
53
54    start = myDistribution.getMean()
55
56    #
57    # Limit state
58    #
59
60    vect = RandomVector(myDistribution)
61
62    output = CompositeRandomVector(limitState, vect)
63
64    myEvent = ThresholdEvent(output, Less(), 0.0)
65
66    #
67    # Computation
68    #
69
70    #
71    # FORM/SORM Cobyla
72    myCobyla = Cobyla()
73    myCobyla.setMaximumEvaluationNumber(100 * dim)
74    myCobyla.setMaximumAbsoluteError(1.0e-10)
75    myCobyla.setMaximumRelativeError(1.0e-10)
76    myCobyla.setMaximumResidualError(1.0e-10)
77    myCobyla.setMaximumConstraintError(1.0e-10)
78
79    myAlgoC = FORM(myCobyla, myEvent, start)
80    myAlgoC2 = SORM(myCobyla, myEvent, start)
81
82    myAlgoC.run()
83    myAlgoC2.run()
84
85    resultC = FORMResult(myAlgoC.getResult())
86    resultC2 = SORMResult(myAlgoC2.getResult())
87
88    #
89    # FORM/SORM Abdo Rackwitz
90    myAbdoRackwitz = AbdoRackwitz()
91    myAbdoRackwitz.setMaximumIterationNumber(100)
92    myAbdoRackwitz.setMaximumAbsoluteError(1.0e-10)
93    myAbdoRackwitz.setMaximumRelativeError(1.0e-10)
94    myAbdoRackwitz.setMaximumResidualError(1.0e-10)
95    myAbdoRackwitz.setMaximumConstraintError(1.0e-10)
96
97    myAlgoAR = FORM(myAbdoRackwitz, myEvent, start + Point(2, 1))
98    myAlgoAR2 = SORM(myAbdoRackwitz, myEvent, start + Point(2, 1))
99
100    myAlgoAR.run()
101    myAlgoAR2.run()
102
103    resultAR = FORMResult(myAlgoAR.getResult())
104    resultAR2 = SORMResult(myAlgoAR2.getResult())
105
106    #
107    # Monte Carlo
108    CoV_MC = 0.1
109    myMC = ProbabilitySimulationAlgorithm(myEvent, MonteCarloExperiment())
110    myMC.setMaximumOuterSampling(1000000)
111    myMC.setBlockSize(1)
112    myMC.setMaximumCoefficientOfVariation(CoV_MC)
113    myMC.run()
114
115    #
116    # LHS
117    CoV_LHS = 0.1
118    myLHS = LHS(myEvent)
119    myLHS.setMaximumOuterSampling(1000000)
120    myLHS.setBlockSize(1)
121    myLHS.setMaximumCoefficientOfVariation(CoV_LHS)
122    myLHS.run()
123
124    #
125    # Results
126    #
127    #
128    # FORM/SORM Cobyla
129    PfC = resultC.getEventProbability()
130    Beta_generalizedC = resultC.getGeneralisedReliabilityIndex()
131    u_starC = resultC.getStandardSpaceDesignPoint()
132    x_starC = resultC.getPhysicalSpaceDesignPoint()
133    PtC = resultC.getIsStandardPointOriginInFailureSpace(
134    ) and "true" or "false"
135    gammaC = resultC.getImportanceFactors()
136    beta_hasoferC = resultC.getHasoferReliabilityIndex()
137    SensitivityC = resultC.getEventProbabilitySensitivity()
138
139    PFBreitC2 = resultC2.getEventProbabilityBreitung()
140    BetaBreitC2 = resultC2.getGeneralisedReliabilityIndexBreitung()
141    PFHBC2 = resultC2.getEventProbabilityHohenbichler()
142    BetaHBC2 = resultC2.getGeneralisedReliabilityIndexHohenbichler()
143    PFTvedtC2 = resultC2.getEventProbabilityTvedt()
144    BetaTvedtC2 = resultC2.getGeneralisedReliabilityIndexTvedt()
145    CurvC2 = resultC2.getSortedCurvatures()
146    u_starC2 = resultC2.getStandardSpaceDesignPoint()
147    x_starC2 = resultC2.getPhysicalSpaceDesignPoint()
148    PtC2 = resultC2.getIsStandardPointOriginInFailureSpace(
149    ) and "true" or "false"
150    gammaC2 = resultC2.getImportanceFactors()
151    beta_hasoferC2 = resultC2.getHasoferReliabilityIndex()
152
153    #
154    # FORM/SORM Abdo Rackwitz
155    PfAR = resultAR.getEventProbability()
156    Beta_generalizedAR = resultAR.getGeneralisedReliabilityIndex()
157    u_starAR = resultAR.getStandardSpaceDesignPoint()
158    x_starAR = resultAR.getPhysicalSpaceDesignPoint()
159    PtAR = resultAR.getIsStandardPointOriginInFailureSpace(
160    ) and "true" or "false"
161    gammaAR = resultAR.getImportanceFactors()
162    beta_hasoferAR = resultAR.getHasoferReliabilityIndex()
163    SensitivityAR = resultAR.getEventProbabilitySensitivity()
164
165    PFBreitAR2 = resultAR2.getEventProbabilityBreitung()
166    BetaBreitAR2 = resultAR2.getGeneralisedReliabilityIndexBreitung()
167    PFHBAR2 = resultAR2.getEventProbabilityHohenbichler()
168    BetaHBAR2 = resultAR2.getGeneralisedReliabilityIndexHohenbichler()
169    PFTvedtAR2 = resultAR2.getEventProbabilityTvedt()
170    BetaTvedtAR2 = resultAR2.getGeneralisedReliabilityIndexTvedt()
171    CurvAR2 = resultAR2.getSortedCurvatures()
172    u_starAR2 = resultAR2.getStandardSpaceDesignPoint()
173    x_starAR2 = resultAR2.getPhysicalSpaceDesignPoint()
174    PtAR2 = resultAR2.getIsStandardPointOriginInFailureSpace(
175    ) and "true" or "false"
176    gammaAR2 = resultAR2.getImportanceFactors()
177    beta_hasoferAR2 = resultAR2.getHasoferReliabilityIndex()
178
179    #
180    # Monte Carlo
181    ResultMC = myMC.getResult()
182    PFMC = ResultMC.getProbabilityEstimate()
183    CVMC = ResultMC.getCoefficientOfVariation()
184    Variance_PF_MC = ResultMC.getVarianceEstimate()
185    length90MC = ResultMC.getConfidenceLength(0.90)
186
187    #
188    # LHS
189    ResultLHS = myLHS.getResult()
190    PFLHS = ResultLHS.getProbabilityEstimate()
191    CVLHS = ResultLHS.getCoefficientOfVariation()
192    Variance_PF_LHS = ResultLHS.getVarianceEstimate()
193    length90LHS = ResultLHS.getConfidenceLength(0.90)
194
195    #
196
197    #
198    # Outputs
199    #
200    print("")
201    print("")
202    print(
203        "************************************************************************************************")
204    print(
205        "***************************************** FORM  COBYLA *****************************************")
206    print(
207        "************************************************************************************************")
208    print("event probability = %.5e" % PfC)
209    print("generalized reliability index = %.5f" % Beta_generalizedC)
210    print(
211        "************************************************************************************************")
212    for i in range(u_starC.getDimension()):
213        print("standard space design point = %.5f" % u_starC[i])
214    print(
215        "************************************************************************************************")
216    for i in range(x_starC.getDimension()):
217        print("physical space design point = %.5f" % x_starC[i])
218    print(
219        "************************************************************************************************")
220    print("is standard point origin in failure space? ", PtC)
221    print(
222        "************************************************************************************************")
223    for i in range(gammaC.getDimension()):
224        print("importance factors = %.5f" % gammaC[i])
225    print(
226        "************************************************************************************************")
227    print("Hasofer reliability index = %.5f" % beta_hasoferC)
228    print(
229        "************************************************************************************************")
230    for i in range(SensitivityC.getSize()):
231        for j in range(SensitivityC[i].getDimension()):
232            print("Pf sensitivity = %.5f" % SensitivityC[i][j])
233    print(
234        "************************************************************************************************")
235    print("")
236    print(
237        "************************************************************************************************")
238    print(
239        "************************************** FORM ABDO RACKWITZ **************************************")
240    print(
241        "************************************************************************************************")
242    print("event probability = %.5e" % PfAR)
243    print("generalized reliability index = %.5f" % Beta_generalizedAR)
244    print(
245        "************************************************************************************************")
246    for i in range(u_starAR.getDimension()):
247        print("standard space design point = %.5f" % u_starAR[i])
248    print(
249        "************************************************************************************************")
250    for i in range(x_starAR.getDimension()):
251        print("physical space design point = %.5f" % x_starAR[i])
252    print(
253        "************************************************************************************************")
254    print("is standard point origin in failure space? ", PtAR)
255    print(
256        "************************************************************************************************")
257    for i in range(gammaAR.getDimension()):
258        print("importance factors = %.5f" % gammaAR[i])
259    print(
260        "************************************************************************************************")
261    print("Hasofer reliability index = %.5f" % beta_hasoferAR)
262    print(
263        "************************************************************************************************")
264    for i in range(SensitivityAR.getSize()):
265        for j in range(SensitivityAR[i].getDimension()):
266            print("Pf sensitivity = %.5f" % SensitivityAR[i][j])
267    print(
268        "************************************************************************************************")
269    print("")
270    print(
271        "************************************************************************************************")
272    print(
273        "***************************************** SORM  COBYLA *****************************************")
274    print(
275        "************************************************************************************************")
276    print("Breitung event probability = %.5e" % PFBreitC2)
277    print("Breitung generalized reliability index = %.5f" % BetaBreitC2)
278    print("Hohenbichler event probability = %.5e" % PFHBC2)
279    print("Hohenbichler generalized reliability index = %.5f" % BetaHBC2)
280    print("Tvedt event probability = %.5e" % PFTvedtC2)
281    print("Tvedt generalized reliability index = %.5f" % BetaTvedtC2)
282    print(
283        "************************************************************************************************")
284    for i in range(CurvC2.getDimension()):
285        print("sorted curvatures = %.5f" % cleanScalar(CurvC2[i]))
286    print(
287        "************************************************************************************************")
288    for i in range(u_starC2.getDimension()):
289        print("standard space design point = %.5f" % u_starC2[i])
290    print(
291        "************************************************************************************************")
292    for i in range(x_starC2.getDimension()):
293        print("physical space design point = %.5f" % x_starC2[i])
294    print(
295        "************************************************************************************************")
296    print(
297        "************************************************************************************************")
298    print("is standard point origin in failure space? ", PtC2)
299    print(
300        "************************************************************************************************")
301    for i in range(gammaC2.getDimension()):
302        print("importance factors = %.5f" % gammaC2[i])
303    print(
304        "************************************************************************************************")
305    print("Hasofer reliability index = %.5f" % beta_hasoferC2)
306    print(
307        "************************************************************************************************")
308    print("")
309    print(
310        "************************************************************************************************")
311    print(
312        "************************************** SORM ABDO RACKWITZ **************************************")
313    print(
314        "************************************************************************************************")
315    print("Breitung event probability = %.5e" % PFBreitAR2)
316    print("Breitung generalized reliability index = %.5f" % BetaBreitAR2)
317    print("Hohenbichler event probability = %.5e" % PFHBAR2)
318    print("Hohenbichler generalized reliability index = %.5f" % BetaHBAR2)
319    print("Tvedt event probability = %.5e" % PFTvedtAR2)
320    print("Tvedt generalized reliability index = %.5f" % BetaTvedtAR2)
321    print(
322        "************************************************************************************************")
323    for i in range(CurvAR2.getDimension()):
324        print("sorted curvatures = %.5f" % cleanScalar(CurvAR2[i]))
325    print(
326        "************************************************************************************************")
327    for i in range(u_starAR2.getDimension()):
328        print("standard space design point = %.5f" % u_starAR2[i])
329    print(
330        "************************************************************************************************")
331    for i in range(x_starAR2.getDimension()):
332        print("physical space design point = %.5f" % x_starAR2[i])
333    print(
334        "************************************************************************************************")
335    print(
336        "************************************************************************************************")
337    print("is standard point origin in failure space? ", PtAR2)
338    print(
339        "************************************************************************************************")
340    for i in range(gammaAR2.getDimension()):
341        print("importance factors = %.5f" % gammaAR2[i])
342    print(
343        "************************************************************************************************")
344    print("Hasofer reliability index = %.5f" % beta_hasoferAR2)
345    print(
346        "************************************************************************************************")
347    print("")
348    print(
349        "************************************************************************************************")
350    print(
351        "**************************************** MONTE CARLO *******************************************")
352    print(
353        "************************************************************************************************")
354    print("Pf estimation = %.5e" % PFMC)
355    print("Pf Variance estimation = %.5e" % Variance_PF_MC)
356    print("CoV = %.5f" % CVMC)
357    print("90% Confidence Interval =", "%.5e" % length90MC)
358    print("CI at 90% =[", "%.5e" % (PFMC - 0.5 * length90MC),
359          "; %.5e" % (PFMC + 0.5 * length90MC), "]")
360    print(
361        "************************************************************************************************")
362    print("")
363    print(
364        "************************************************************************************************")
365    print(
366        "******************************************* L H S **********************************************")
367    print(
368        "************************************************************************************************")
369    print("Pf estimation = %.5e" % PFLHS)
370    print("Pf Variance estimation = %.5e" % Variance_PF_LHS)
371    print("CoV = %.5f" % CVLHS)
372    print("90% Confidence Interval =", "%.5e" % length90LHS)
373    print("CI at 90% =[", "%.5e" % (PFLHS - 0.5 * length90LHS),
374          "; %.5e" % (PFLHS + 0.5 * length90LHS), "]")
375    print(
376        "************************************************************************************************")
377    print("")
378
379except:
380    import sys
381    print("t_Waarts_system_series.py", sys.exc_info()[0], sys.exc_info()[1])
382