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