ExampleAxialStressedBeam: AxialStressedBeam-v1.2.py

File AxialStressedBeam-v1.2.py, 4.1 KB (added by schueller, 4 years ago)
Line 
1#! /usr/bin/env python
2
3# import OpenTURNS module
4from openturns import *
5
6# import optional OpenTURNS viewer module
7from openturns.viewer import View
8
9# dimension
10dim = 2
11
12# Random generator initialization:
13RandomGenerator.SetSeed(0)
14
15# Analytical model definition:
16limitState = NumericalMathFunction(['R', 'F'], ['G'], ['R-F/(_pi*100.0)'])
17
18# Test of the limit state function:
19x = [300., 75000.]
20print 'x=', x
21print '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
27R_dist = LogNormal(300., 30., 0., LogNormal.MUSIGMA)
28R_dist.setName('Yield strength')
29R_dist.setDescription('R')
30# Graphical output of the PDF
31R_pdf = R_dist.drawPDF()
32View(R_pdf).save('pdf_R.png')
33
34# Create a second marginal : Normal distribution 1D
35F_dist = Normal(75000., 5000.)
36F_dist.setName('Traction_load')
37F_dist.setDescription('F')
38# Graphical output of the PDF
39F_pdf = F_dist.drawPDF()
40View(F_pdf).save('pdf_F.png')
41
42# Create a copula : IndependentCopula (no correlation)
43aCopula = IndependentCopula(dim)
44aCopula.setName('Independent copula')
45
46# Instanciate one distribution object
47myDistribution = ComposedDistribution([R_dist, F_dist], aCopula)
48myDistribution.setName('myDist')
49
50# We create a 'usual' RandomVector from the Distribution
51vect = RandomVector(myDistribution)
52
53# We create a composite random vector
54G = RandomVector(limitState, vect)
55
56# We create an Event from this RandomVector
57myEvent = Event(G, Less(), 0.0, 'Event 1')
58
59# Using Monte Carlo simulations
60cv = 0.05
61NbSim = 100000
62
63algoMC = MonteCarlo(myEvent)
64algoMC.setMaximumOuterSampling(NbSim)
65algoMC.setBlockSize(1)
66algoMC.setMaximumCoefficientOfVariation(cv)
67# For statistics about the algorithm
68initialNumberOfCall = limitState.getEvaluationCallsNumber()
69
70# Perform the analysis:
71algoMC.run()
72
73# Results:
74result = algoMC.getResult()
75probability = result.getProbabilityEstimate()
76print 'MonteCarlo result=', result
77print 'Number of executed iterations =', result.getOuterSampling()
78print 'Number of calls to the limit state =', limitState.getEvaluationCallsNumber() - initialNumberOfCall
79print 'Pf = ', probability
80print 'CV =', result.getCoefficientOfVariation()
81MCcvgraph = algoMC.drawProbabilityConvergence()
82View(MCcvgraph).save('montecarlo_convergence.png')
83
84# Using FORM analysis
85
86# We create a NearestPoint algorithm
87myCobyla = Cobyla()
88# Resolution options:
89eps = 1e-3
90myCobyla.setMaximumIterationsNumber(100)
91myCobyla.setMaximumAbsoluteError(eps)
92myCobyla.setMaximumRelativeError(eps)
93myCobyla.setMaximumResidualError(eps)
94myCobyla.setMaximumConstraintError(eps)
95
96# For statistics about the algorithm
97initialNumberOfCall = 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
104algoFORM = FORM(myCobyla, myEvent, myDistribution.getMean())
105
106# Perform the analysis:
107algoFORM.run()
108
109# Results:
110result = algoFORM.getResult()
111print 'Number of calls to the limit state =', limitState.getEvaluationCallsNumber() - initialNumberOfCall
112print 'Pf =', result.getEventProbability()
113
114# Graphical result output
115importanceFactorsGraph = result.drawImportanceFactors()
116View(importanceFactorsGraph).save('importance_factors.png')
117
118# Using Directional sampling
119
120# Resolution options:
121cv = 0.05
122NbSim = 100000
123
124algoDS = DirectionalSampling(myEvent)
125algoDS.setMaximumOuterSampling(NbSim)
126algoDS.setBlockSize(1)
127algoDS.setMaximumCoefficientOfVariation(cv)
128# For statistics about the algorithm
129initialNumberOfCall = limitState.getEvaluationCallsNumber()
130
131# Perform the analysis:
132algoDS.run()
133
134# Results:
135result = algoDS.getResult()
136probability = result.getProbabilityEstimate()
137print 'Number of executed iterations =', result.getOuterSampling()
138print 'Number of calls to the limit state =', limitState.getEvaluationCallsNumber() - initialNumberOfCall
139print 'Pf = ', probability
140print 'CV =', result.getCoefficientOfVariation()
141DScvgraph = algoDS.drawProbabilityConvergence()
142View(DScvgraph).save('directionalsampling_convergence.png')