Monte Carlo simulation [elastic portal frame]

[mgrubisic]

I'm trying to perform a simple Monte Carlo simulation on an existing example "portalframe.tcl" available in OpenSees repository. The model is simple and clean, composed of elastic elements and also runs FORM analysis. The values of the modulus of elasticity (E) and a horizontal force (P) are defined as random variables.

Monte Carlo simulation works if only one random variable is active (E or P in this case), and thus MCS (for a reasonable number of trials) gives equal probability of failure as a FORM analysis. After activation of both random variables (E and P) - the value of horizontal displacement of node #2 (on the basis of which the performance function is defined) increased by almost a million times in MCS, while FORM analysis works just fine. I do not understand why this is happening?

The problem may be trivial, but please help me in any way! The model is given below.

----------------------------------------------------------------------------------------------------------------------------
wipe
wipeReliability

# +------------------------+
# | MODEL BUILD |
# +------------------------+

model basic -ndm 2 -ndf 3

node 1 0 0
node 2 0 144
node 3 240 144
node 4 240 0

fix 1 1 1 1
fix 4 1 1 1

set E 30000
set Ag 25
set Ig 1500
set Ac 29
set Ic 2000

section Elastic 1 $E $Ag $Ig
section Elastic 2 $E $Ac $Ic

geomTransf Linear 1

element forceBeamColumn 1 1 2 1 Lobatto 2 3
element forceBeamColumn 2 2 3 1 Lobatto 1 3
element forceBeamColumn 3 3 4 1 Lobatto 2 3

set P 25.0
set w 1.0e-1
pattern Plain 1 Constant {
load 2 $P 0 0
eleLoad -ele 2 -type beamUniform -$w
}

# ======== Static Analysis Setup ========
system BandGeneral
constraints Transformation
numberer RCM

set Tol 1.e-8;
set maxNumIter 6;
set printFlag 0;
set TestType EnergyIncr;
test $TestType $Tol $maxNumIter $printFlag;

algorithm Newton
integrator LoadControl 1.0
analysis Static


# +----------------------------------+
# | RELIABILITY ANALYSIS |
# +----------------------------------+

reliability

# randomVariable 62 lognormal -mean $E -stdv [expr 0.3*$E]
randomVariable 32 normal -mean $P -stdv [expr 0.3*$P]

# parameter 62 randomVariable 62 element 1 E
# addToParameter 62 element 2 E
# addToParameter 62 element 3 E
parameter 32 randomVariable 32 loadPattern 1 loadAtNode 2 1

parameter 23 node 2 disp 1

performanceFunction 76 "0.15-\$par(23)"


# ======== Reliability Analysis Setup ========
sensitivityIntegrator -static
sensitivityAlgorithm -computeAtEachStep

randomNumberGenerator CStdLib
probabilityTransformation Nataf -print 3
reliabilityConvergenceCheck Standard -e1 1.0e-2 -e2 1.0e-2 -print 1
functionEvaluator Tcl -file "analyze 1"
gradientEvaluator Implicit
searchDirection iHLRF
meritFunctionCheck AdkZhang -multi 2.0 -add 50 -factor 0.5
# stepSizeRule Armijo -maxNum 10 -base 0.5 -initial 0.3 5
stepSizeRule Fixed -stepSize 1.0
startPoint Mean
findDesignPoint StepSearch -maxNumIter 30; # -printDesignPointX designPointX.out

# ======== Run FORM Analysis ========
runFORMAnalysis portalframe.out

foreach perf [getLSFTags] {
puts "Performance Function $perf"
puts "FORM beta = [format %.7f $betaFORM($perf)]"
foreach rv [getRVTags] {
puts "\t x*($rv) = [format %7.4f $designPointXFORM($perf,$rv)], alpha($rv) = [format %.7f $alphaFORM($perf,$rv)], gamma($rv) = [format %.7f $gammaFORM($perf,$rv)]"
}
}

puts ""
set beta $betaFORM(76); # betaFORM (lsfTag) gives the beta value for LSF $lsfTag
puts "FORM beta = [format %.7f $betaFORM($perf)]"
puts "FORM Pf = [getStdNormalCDF -$beta]"; # pf = Phi(-beta)


# +----------------------------------------+
# | MONTE CARLO SIMULATION |
# +----------------------------------------+

set Ntrials 10000
set Nfail 0

puts ""

for {set i 1} {$i <= $Ntrials} {incr i} {
reset
set U ""
foreach rv [getRVTags] {
set val [expr rand()]
lappend U [getStdNormalInverseCDF $val]
}

set X [transformUtoX $U]
set irv 0

foreach rv [getRVTags] {
updateParameter $rv [lindex $X $irv]
incr irv
}

analyze 1
print -node -flag 1 2
set D [nodeDisp 2 1]
if {[expr 0.15 - $D] <= 0.0} {incr Nfail}
}

set MCSPf [expr double($Nfail)/($Ntrials)];
puts ""
puts "Number of Trials = $Ntrials"
puts "Monte Carlo Pf = [format %.10f $MCSPf]"
puts ""
puts "FORM beta = [format %.7f $betaFORM($perf)]"
puts "FORM Pf = [getStdNormalCDF -$beta]"; # pf = Phi(-beta)

puts ""
puts "========== End of Reliability Analysis =========="
puts ""

[mgrubisic]

Reply by Michael H. Scott via e-mail, quote:

"It seems you have come across a bug in the transformUtoX (and transformXtoU) functions. Basically, it is switching the order of your random variables internally, resulting in E taking on the values of P and vice versa. I’ll work on a fix. In the meantime, you can fix your problem easily by defining random variable 32 before 62."

Problem solved!