Monte Carlo Simulation

The uncertainty associated with a finite element analysis is as important, if not more important, than the results of the analysis itself. Thanks to Terje Haukaas, OpenSees has several modules for finite element reliability analysis: FORM, FOSM, SORM, and several other methods to quantify uncertainty. Unfortunately, those methods have not yet made their way into OpenSeesPy, although you can check out this branch for progress.

However, a few reliability commands have made it into OpenSeesPy, including randomVariable and probabilityTransformation, which are all you need in order to perform Monte Carlo simulation.

In this post, I’ll show Monte Carlo simulation for a simple bar model–a “Little MCS” if you will. You can find the Tcl version of this analysis on GitHub and also on YouTube. Seweryn Kokot helped me get the Python version of the analysis up and running, so thanks, Seweryn!

The bar has capacity and is subjected to a load (both normal distributions). The performance, or limit state, function g=R-S, defines the fail (g<=0) and safe (g>0) states for the bar.

The Monte Carlo simulation will be carried out with two approaches: 1) with basic random variables and 2) with random variables mapped to a finite element model of the bar.

Using only basic random variables in the first approach, you don’t really need OpenSees–you could do it all in Python. However, one advantage with OpenSees is its implementation of the Nataf probability transformation for correlated random variables.

Within the Monte Carlo loop, realizations of standard normal random variables (zero mean, unit standard deviation) are generated, then transformed to the original random variable space using the transformUtoX command. Then, the performance function is evaluated and failure states are counted.

import openseespy.opensees as ops
from scipy.stats import norm

meanR = 2000; sigR = 135  # Resistance
meanS = 1700; sigS = 200  # Load

ops.randomVariable(1, 'normal', '-mean', meanR, '-stdv', sigR)
ops.randomVariable(2, 'normal', '-mean', meanS, '-stdv', sigS)

# pf ~ 0.04 or ~ 0.1 if commented (uncorrelated RVs)
ops.correlate(1, 2, 0.5)
ops.probabilityTransformation('Nataf')

nrv = len(ops.getRVTags())
nTrials = 50000
nFail = 0
for i in range(nTrials):
    U = list(norm.rvs(size=nrv)) # RVs on N(0,1) 
    X = ops.transformUtoX(*U)
    r = X[0]
    s = X[1]
    
    g = r-s
    if (g <= 0.0):
        nFail += 1

MCpf = nFail / nTrials
print(f'Monte Carlo simulation, pf_MC = {nFail} / {nTrials} = {MCpf}')

The probability of failure should be about 0.04. If you make the random variables uncorrelated, the probability of failure will be about 0.1. Try this analysis out for yourself.

For the second approach, a finite element model of the bar is built with length and cross-section area both equal to one. Because OpenSees solves equilibrium, evaluating the performance function, g=R-S, is not very meaningful. Instead, the force-deformation response of the bar is defined as bilinear with yield strength R and the performance function becomes g=R/k-U where k is the bar stiffness and U is the bar displacement. Because the bar is statically determinate, the magnitude of the bar stiffness doesn’t really matter.

After defining the model, two parameters are defined–one for the yield strength of the bar and the other for the applied load. Then, just like in the first approach, realizations of the random variables are generated with standard normals and the Nataf transformation. These realizations are fed into the OpenSees model with the updateParameter command, then the performance function is evaluated and fail/safe states are determined. Each time through the Monte Carlo loop, the OpenSees model is brought back to its initial state using the reset command.

import openseespy.opensees as ops
from scipy.stats import norm

meanR = 2000; sigR = 135  # Resistance
meanS = 1700; sigS = 200  # Load

ops.randomVariable(1, 'normal', '-mean', meanR, '-stdv', sigR)
ops.randomVariable(2, 'normal', '-mean', meanS, '-stdv', sigS)

# pf ~ 0.04 or ~ 0.1 if commented (uncorrelated RVs)
ops.correlate(1, 2, 0.5)
ops.probabilityTransformation('Nataf')

ops.wipe()
ops.model('basic','-ndm',1)

ops.node(1,0); ops.fix(1,1)
ops.node(2,1)

k = 1000 # Bar stiffness (doesn't really matter)
ops.uniaxialMaterial('Hardening',1,k,meanR,0,0.01*k)
ops.element('truss',1,1,2,1.0,1)

ops.timeSeries('Constant',1)
ops.pattern('Plain',1,1)
ops.load(2,meanS)

ops.parameter(1,'element',1,'Fy')
ops.parameter(2,'loadPattern',1,'loadAtNode',2,1)

ops.analysis('Static')

nrv = len(ops.getRVTags())
nTrials = 50000
nFail = 0
for i in range(nTrials):
    ops.reset()
    
    U = list(norm.rvs(size=nrv))  
    X = ops.transformUtoX(*U)
    r = X[0]
    s = X[1]
    ops.updateParameter(1,r)
    ops.updateParameter(2,s)
    
    ops.analyze(1)
    
    g = r/k - ops.nodeDisp(2,1)
    if (g <= 0.0):
        nFail += 1

MCpf = nFail / nTrials
print(f'Monte Carlo simulation, pf_MC = {nFail} / {nTrials} = {MCpf}')

You can verify that you get the same probability of failure as in the first approach.

For larger finite element models with more random variables, there are better ways to handle the parameter updating, but this simple bar analysis shows the important aspects of Monte Carlo simulation with OpenSees.