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