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