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