Workshop 5: Deterministic Analysis (notebook 2 of 2)

Workshop 5: Deterministic Analysis (notebook 2 of 2)#

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15, 6

import openturns as ot
ot.Log.Show(ot.Log.NONE)

import Probabilistic

Task 1: Set up deterministic calculation#

Everything needed to describe stress in the tunnels is included in the *.py files that are imported below. See notebook Analysis_deterministic.ipynb for additional information.

Task 1: confirm that the Python modules can be imported successfully.

import FEM
import Morrisson_forces

Task 2: Determine statistical properties of Random Variables#

Task 2.1: define the marginal distributions using OpenTURNs (the tuple containing the random variables and their descriptions is already provided).

# YOUR_CODE_HERE

# Solution:
X1 = ot.Exponential(2.0309)
X2 = ot.Exponential(2.0663)
X3 = ot.Normal(2, 0.2)

X = (X1, X2, X3)
descriptions = ["windsea","swell", "u_current"]

Task 2.2: visualize each of the random variables to confirm that they are implemented correctly, and that you get a sense for their magnitude and distribution.

X1_pdf = X1.drawPDF()
graph = X1.drawPDF()
graph.setTitle("PDF of X1, significant wave height, swell waves")
graph.setXTitle("x")
graph.setYTitle("Density")
view = ot.viewer.View(graph)
plt.show()

../../../_images/9cdbc21d55fe9c6fea5e2b402c20afc158b96023b68da2b692b02462c7543f56.png

Task 3: Determine the Limit State Function#

To keep things easy to follow, we will divide this into two parts:

write a function to calculate the total axial: total_calculation_axial_stresses
write a function to represent the vectorized LSF, as in WS04: myLSF

Task 3.1: complete the function as done in notebook Analysis_deterministic.ipynb to evaluate the axial stress.

# def total_calculation_axial_stresses(significant_wave_height_swell,
#                                      significant_wave_height_windsea,
#                                      U_current_velocity):
#     '''Find axial stresses due to wind and swell waves.
    
#     Inputs: three load random variables
#     Returns: stress
#     '''
    
#     YOUR_CODE_FROM_DETERMINISTIC_NOTEBOOK_HERE
#     # Note: you do not need to include the `time` tracking here.

# SOLUTION:
def total_calculation_axial_stresses(significant_wave_height_swell,
                                     significant_wave_height_windsea,
                                     U_current_velocity):
    '''Find axial stresses due to wind and swell waves.
    
    Inputs: three load random variables
    Returns: stress
    '''
    sigma_wind = FEM.calculate_stress_FEM(
        Morrisson_forces.F_morison_wind(
            significant_wave_height_windsea))
    
    sigma_swell_and_current = FEM.calculate_stress_FEM(
        Morrisson_forces.F_morison_swell_current(
            significant_wave_height_swell,
            U_current_velocity))
    
    # Total stress
    sigma = sigma_wind + sigma_swell_and_current 
    
    return sigma 

Task 3.2: based on the definition of failure provided above, and your understanding of the stress calculation algorithm (notebook Analysis_deterministic.ipynb), definte the limit-state function. Remember to vectorize the function.

# def myLSF(x):
#     ''' 
#     Vectorized limit-state function.

#     Inputs: x (array of random variables)
#     Output: limit-state function value
#     '''
#     YOUR_CODE_HERE

# Solution:
def myLSF(x):
    ''' 
    Vectorized limit-state function.

    Inputs: x (array of random variables)
    Output: limit-state function value
    '''
    sigma_max = 75 # MPa, maximum allowable compression stress in special concrete
    g = [sigma_max - np.max(total_calculation_axial_stresses(x[0], x[1], x[2]))]

    return g

Task 3.3: we need to define failure for OpenTURNs. Use the cell below to define the failure criteria for use in OpenTURN's.

Hint: look at the inputs for function input_OpenTurns, which is defined inside the file Probabilistic.py. Executing the code cell below will print out the docstring.

help(Probabilistic.input_OpenTurns)
# YOUR_CODE_HERE

# Solution:
failure_threshold = 0

Help on function input_OpenTurns in module Probabilistic:

input_OpenTurns(X, descriptions, myLSF, failure_threshold)

Task 4: Run Monte Carlo Simulation and FORM analysis#

Most of the work has been done for you inside the file Probabilistic.py. It reuses the same code for using OpenTURNs that you saw previously in WS03 and WS04.

Task 4.1: read through the functions in Probabilistic.py until you understand what each one does and how you can use it to solve the component reliability problem. Test your understanding by taking turns explaining the contents to your fellow group members.

You now have everything you need to complete the component reliability problem. All you have to do is execute the right code, interpret the results, then write about it in Report.md.

Note that you will have to decide yourself how to use each method (e.g., FORM versus MCS). Note that the FEM model takes some time to run, so using a lot of simulations with MCS is probably not a good way to start.

Before you start running code, take a few minutes to review the following overview pages in the workbook:

Task 4.2:
Use the functions in Probabilistic.py to complete the MCS and FORM analysis. Because our target probability is 1%, we should run MCS for at least 1000 simulations. Think carefully about how long that analysis might take before you start it (for example, how long did each analysis take in the deterministic notebook?). Perhaps it is better to start with a different (faster) method?

Note that the task of setting up the reliability analysis will be much easier if you actually spent some time on Task 4.1 😜

You are free to modify the code in Probabilistic.py, but it is not required in order to complete the assignment.

# YOUR_CODE_HERE

# Solution:
Probabilistic.input_OpenTurns(X, descriptions, myLSF, failure_threshold)
result, x_star, u_star, pf_FORM, beta = Probabilistic.run_FORM_analysis()
alpha_ot, alpha, sens = Probabilistic.importance_factors(result)

# for illustration purposes. don't do this because it will take 20-30 minutes!
# pf_mc = Probabilistic.run_MonteCarloSimulation(mc_size=1000)

The FORM analysis took 138.861 seconds
FORM result, pf = 0.1845
FORM result, beta = 0.898

The design point in the u space:  [0.594652,0.672709,0.0296273]
The design point in the x space:  [0.633816,0.669812,2.00593]
--- FORM Importance Factors (alpha) ---


Importance factors, from OpenTURNs:
   0.438
   0.561
   0.001

Importance factors, based on normal vector in U-space = 
   0.662
   0.749
   0.033
Note: this will be different from      result.getImportanceFactors()      if there are resistance variables.

Sensitivity of Reliability Index to Multivariate Distribution

Distribution item number: 0
  Item name: windsea

Distribution item number: 1
  Item name: swell

Distribution item number: 2
  Item name: u_current

Distribution item number: 3
  Item name: NormalCopula
    +0.000e+00 for parameter R_2_1_copula
    +0.000e+00 for parameter R_3_1_copula
    +0.000e+00 for parameter R_3_2_copula

Task 5: Design Case#

The consortium of Colomes & Agarwal Inc. and Frank & Co is still trying to decide what material to use for the tunnel: concrete or steel. They think the current concrete design is insufficient, so they have asked you to look into the matter for them. There are two options:

Improve the quality of the concrete production process, which reduces \(\sigma\) from 5 to 1
Change the material to steel (distribution defined in README.md)

It is not immediately clear which will have a bigger impact, as the standard deviation is much higher for steel, despite the higher mean value.

Task 5: add a fourth variable to the reliability analysis and run a few more reliability calculations. Note that you will need to modify several Python objects:the create a new random varable X4, redefine the vector of random variables X, redefine the descriptions and create a new limit-state function.

# Solution:
# X4 = ot.Normal(75, 5)
# X4 = ot.Normal(75, 1)
# X4 = ot.Normal(250, 20)

X = (X1, X2, X3, X4)
descriptions = ["windsea","swell", "u_current", "sigma_max"]

def myLSF2(x):
    ''' 
    Vectorized limit-state function.

    Inputs: x (array of random variables)
    Output: limit-state function value
    '''
    sigma_max = x[3] # MPa, maximum allowable compression stress in special concrete
    g = [sigma_max - np.max(total_calculation_axial_stresses(x[0], x[1], x[2]))]

    return g

Probabilistic.input_OpenTurns(X, descriptions, myLSF2, failure_threshold)
result, x_star, u_star, pf_FORM, beta = Probabilistic.run_FORM_analysis()
alpha_ot, alpha, sens = Probabilistic.importance_factors(result)

# for illustration purposes. don't do this because it will take 20-30 minutes!
# pf_mc = Probabilistic.run_MonteCarloSimulation(mc_size=1000)

The FORM analysis took 150.177 seconds
FORM result, pf = 0.1846
FORM result, beta = 0.898

The design point in the u space:  [0.594536,0.672218,0.0294498,-0.0201069]
The design point in the x space:  [0.633747,0.66951,2.00589,74.9799]
--- FORM Importance Factors (alpha) ---


Importance factors, from OpenTURNs:
   0.438
   0.560
   0.001
   0.001

Importance factors, based on normal vector in U-space = 
   0.662
   0.749
   0.033
  -0.022
Note: this will be different from      result.getImportanceFactors()      if there are resistance variables.

Sensitivity of Reliability Index to Multivariate Distribution

Distribution item number: 0
  Item name: windsea

Distribution item number: 1
  Item name: swell

Distribution item number: 2
  Item name: u_current

Distribution item number: 3
  Item name: sigma_max

Distribution item number: 4
  Item name: NormalCopula
    +0.000e+00 for parameter R_2_1_copula
    +0.000e+00 for parameter R_3_1_copula
    +0.000e+00 for parameter R_3_2_copula
    +0.000e+00 for parameter R_4_1_copula
    +0.000e+00 for parameter R_4_2_copula
    +0.000e+00 for parameter R_4_3_copula