# ExampleAxialStressedBeam: AxialStressedBeam-v1.4.py

File AxialStressedBeam-v1.4.py, 4.0 KB (added by schueller, 4 years ago) |
---|

Line | |
---|---|

1 | #! /usr/bin/env python |

2 | |

3 | # import OpenTURNS module |

4 | from openturns import * |

5 | |

6 | # import optional OpenTURNS viewer module |

7 | from openturns.viewer import View |

8 | |

9 | # dimension |

10 | dim = 2 |

11 | |

12 | # Random generator initialization: |

13 | RandomGenerator.SetSeed(0) |

14 | |

15 | # Analytical model definition: |

16 | limitState = NumericalMathFunction(['R', 'F'], ['G'], ['R-F/(_pi*100.0)']) |

17 | |

18 | # Test of the limit state function: |

19 | x = [300., 75000.] |

20 | print 'x=', x |

21 | print 'G(x)=', limitState(x) |

22 | |

23 | # Stochastic model definition: |

24 | |

25 | # Create a first marginal : LogNormal distribution 1D, parameterized by |

26 | # its mean and standard deviation |

27 | R_dist = LogNormal(300., 30., 0., LogNormal.MUSIGMA) |

28 | R_dist.setName('Yield strength') |

29 | R_dist.setDescription('R') |

30 | # Graphical output of the PDF |

31 | R_pdf = R_dist.drawPDF() |

32 | View(R_pdf).save('pdf_R.png') |

33 | |

34 | # Create a second marginal : Normal distribution 1D |

35 | F_dist = Normal(75000., 5000.) |

36 | F_dist.setName('Traction_load') |

37 | F_dist.setDescription('F') |

38 | # Graphical output of the PDF |

39 | F_pdf = F_dist.drawPDF() |

40 | View(F_pdf).save('pdf_F.png') |

41 | |

42 | # Create a copula : IndependentCopula (no correlation) |

43 | aCopula = IndependentCopula(dim) |

44 | aCopula.setName('Independent copula') |

45 | |

46 | # Instanciate one distribution object |

47 | myDistribution = ComposedDistribution([R_dist, F_dist], aCopula) |

48 | myDistribution.setName('myDist') |

49 | |

50 | # We create a 'usual' RandomVector from the Distribution |

51 | vect = RandomVector(myDistribution) |

52 | |

53 | # We create a composite random vector |

54 | G = RandomVector(limitState, vect) |

55 | |

56 | # We create an Event from this RandomVector |

57 | myEvent = Event(G, Less(), 0.0) |

58 | |

59 | # Using Monte Carlo simulations |

60 | cv = 0.05 |

61 | NbSim = 100000 |

62 | |

63 | algoMC = MonteCarlo(myEvent) |

64 | algoMC.setMaximumOuterSampling(NbSim) |

65 | algoMC.setBlockSize(1) |

66 | algoMC.setMaximumCoefficientOfVariation(cv) |

67 | # For statistics about the algorithm |

68 | initialNumberOfCall = limitState.getEvaluationCallsNumber() |

69 | |

70 | # Perform the analysis: |

71 | algoMC.run() |

72 | |

73 | # Results: |

74 | result = algoMC.getResult() |

75 | probability = result.getProbabilityEstimate() |

76 | print 'MonteCarlo result=', result |

77 | print 'Number of executed iterations =', result.getOuterSampling() |

78 | print 'Number of calls to the limit state =', limitState.getEvaluationCallsNumber() - initialNumberOfCall |

79 | print 'Pf = ', probability |

80 | print 'CV =', result.getCoefficientOfVariation() |

81 | MCcvgraph = algoMC.drawProbabilityConvergence() |

82 | View(MCcvgraph).save('montecarlo_convergence.png') |

83 | |

84 | # Using FORM analysis |

85 | |

86 | # We create a NearestPoint algorithm |

87 | myCobyla = Cobyla() |

88 | # Resolution options: |

89 | eps = 1e-3 |

90 | myCobyla.setMaximumIterationsNumber(100) |

91 | myCobyla.setMaximumAbsoluteError(eps) |

92 | myCobyla.setMaximumRelativeError(eps) |

93 | myCobyla.setMaximumResidualError(eps) |

94 | myCobyla.setMaximumConstraintError(eps) |

95 | |

96 | # For statistics about the algorithm |

97 | initialNumberOfCall = limitState.getEvaluationCallsNumber() |

98 | |

99 | # We create a FORM algorithm |

100 | # The first parameter is a NearestPointAlgorithm |

101 | # The second parameter is an event |

102 | # The third parameter is a starting point for the design point research |

103 | |

104 | algoFORM = FORM(myCobyla, myEvent, myDistribution.getMean()) |

105 | |

106 | # Perform the analysis: |

107 | algoFORM.run() |

108 | |

109 | # Results: |

110 | result = algoFORM.getResult() |

111 | print 'Number of calls to the limit state =', limitState.getEvaluationCallsNumber() - initialNumberOfCall |

112 | print 'Pf =', result.getEventProbability() |

113 | |

114 | # Graphical result output |

115 | importanceFactorsGraph = result.drawImportanceFactors() |

116 | View(importanceFactorsGraph).save('importance_factors.png') |

117 | |

118 | # Using Directional sampling |

119 | |

120 | # Resolution options: |

121 | cv = 0.05 |

122 | NbSim = int(1e5) |

123 | |

124 | algoDS = DirectionalSampling(myEvent) |

125 | algoDS.setMaximumOuterSampling(NbSim) |

126 | algoDS.setBlockSize(1) |

127 | algoDS.setMaximumCoefficientOfVariation(cv) |

128 | # For statistics about the algorithm |

129 | initialNumberOfCall = limitState.getEvaluationCallsNumber() |

130 | |

131 | # Perform the analysis: |

132 | algoDS.run() |

133 | |

134 | # Results: |

135 | result = algoDS.getResult() |

136 | probability = result.getProbabilityEstimate() |

137 | print 'Number of executed iterations =', result.getOuterSampling() |

138 | print 'Number of calls to the limit state =', limitState.getEvaluationCallsNumber() - initialNumberOfCall |

139 | print 'Pf = ', probability |

140 | print 'CV =', result.getCoefficientOfVariation() |

141 | DScvgraph = algoDS.drawProbabilityConvergence() |

142 | View(DScvgraph).save('directionalsampling_convergence.png') |