Permalink
Branch: master
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
154 lines (142 sloc) 5.65 KB
#!/usr/bin/env python
# -*- encoding: utf-8 -*-
"""Load and Resistance probability density and Safety Margins per ISO 19902:2007
2018 ckunte
Usage: pdf.py ( --gom | --nns | --cns | --aus ) [--exp=L]
pdf.py --help
pdf.py --version
Options:
-h, --help Show help screen
--gom Plot for Gulf of Mexico
--nns Plot for Northern North Sea
--cns Plot for Central North Sea (also applicable to SNS)
--aus Plot for Australian Northwest shelf
--exp=L Exposure level: 1 for L1, 2 for L2 [default: 1]
--version Show version
"""
from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt
from docopt import docopt
args = docopt(__doc__, version='Load and Resistance probability density and Safety Margins per ISO 19902:2007, ver:0.3')
cat = int(args['--exp'])
"""
Legend:
Em -- Mean load for a reference return period
Ve -- Load COV
Rm -- Mean resistance
Vr -- Resistance COV
COV = [0.212, 0.265, 0.32, 0.33] # For CNS & SNS, NNS, GOM, AUS
Ve = map(sqrt(x^2 + 0.08^2), COV) => [0.2266, 0.2768, 0.3298, 0.3396]
Rm = RSR_mean = [1.73, 1.92, 1.85, 2.18] for target Pf = 3E-5/annum (used in these plots)
Rm = RSR_mean = [1.40, 1.49, 1.60, 1.60] for target Pf = 5E-4/annum
"""
def ER(Em,Ve,Rm,Vr):
mu_e = np.log(Em / np.sqrt(1 + Ve**2)) # mean
mu_r = np.log(Rm / np.sqrt(1 + Vr**2)) # mean
sigma_e = np.sqrt(np.log(1 + Ve**2))
sigma_r = np.sqrt(np.log(1 + Vr**2))
#
s_e = np.random.lognormal(mu_e, sigma_e, 1000)
s_r = np.random.lognormal(mu_r, sigma_r, 1000)
count_e, bins_e, ignored_e = plt.hist(s_e, 100, normed=True, align='mid', alpha=0.09, color='r')
count_r, bins_r, ignored_r = plt.hist(s_r, 100, normed=True, align='mid', alpha=0.09, color='g')
x_e = np.linspace(min(bins_e), max(bins_e), 10000)
x_r = np.linspace(min(bins_r), max(bins_r), 10000)
#
pdf_E = (np.exp(-(np.log(x_e) - mu_e)**2 / (2 * sigma_e**2)) / (x_e * sigma_e * np.sqrt(2 * np.pi)))
pdf_R = (np.exp(-(np.log(x_r) - mu_r)**2 / (2 * sigma_r**2)) / (x_r * sigma_r * np.sqrt(2 * np.pi)))
#
plt.plot(x_e, pdf_E, linewidth=1, color='r')
plt.plot(x_r, pdf_R, linewidth=1, color='g')
plt.ylabel('Probability density')
misc()
pass
# Data: Em, Ve, Rm, Vr
d = [
0.79,0.3298,1.85,0.05, # L1: GOM
0.79,0.3298,1.60,0.05, # L2: GOM
0.81,0.2768,1.92,0.05, # L1: NNS
0.81,0.2768,1.49,0.05, # L2: NNS
0.84,0.2266,1.73,0.05, # L1: CNS and SNS
0.84,0.2266,1.40,0.05, # L2: CNS and SNS
0.78,0.3396,2.18,0.05, # L1: AUS
0.78,0.3396,1.60,0.05 # L2: AUS
]
def gom():
plt.title('Gulf of Mexico')
if cat == 1:
ER(d[0],d[1],d[2],d[3])
plt.xlabel('Load or resistance as times the nominal load, $x, E_{mean}(%.2f); x, R_{mean}(%.2f)$' %(d[0],d[2]))
plt.axvline(x=d[0],color='r',linestyle=':')
plt.axvline(x=d[2],color='g',linestyle=':')
plt.savefig('GoM_L1.png')
if cat == 2:
ER(d[4],d[5],d[6],d[7])
plt.xlabel('Load or resistance as times the nominal load, $x, E_{mean}(%.2f); x, R_{mean}(%.2f)$' %(d[4],d[6]))
plt.axvline(x=d[4],color='r',linestyle=':')
plt.axvline(x=d[6],color='g',linestyle=':')
plt.savefig('GoM_L2.png')
pass
def nns():
plt.title('Northern North Sea')
if cat == 1:
ER(d[8],d[9],d[10],d[11])
plt.xlabel('Load or resistance as times the nominal load, $x, E_{mean}(%.2f); x, R_{mean}(%.2f)$' %(d[8],d[10]))
plt.axvline(x=d[8],color='r',linestyle=':')
plt.axvline(x=d[10],color='g',linestyle=':')
plt.savefig('NNS_L1.png')
if cat == 2:
ER(d[12],d[13],d[14],d[15])
plt.xlabel('Load or resistance as times the nominal load, $x, E_{mean}(%.2f); x, R_{mean}(%.2f)$' %(d[12],d[14]))
plt.axvline(x=d[12],color='r',linestyle=':')
plt.axvline(x=d[14],color='g',linestyle=':')
plt.savefig('NNS_L2.png')
pass
def cns_sns():
plt.title('Central North Sea (also applicable to Southern North Sea)')
if cat == 1:
ER(d[16],d[17],d[18],d[19])
plt.xlabel('Load or resistance as times the nominal load, $x, E_{mean}(%.2f); x, R_{mean}(%.2f)$' %(d[16],d[18]))
plt.axvline(x=d[16],color='r',linestyle=':')
plt.axvline(x=d[18],color='g',linestyle=':')
plt.savefig('CNS-SNS_L1.png')
if cat == 2:
ER(d[20],d[21],d[22],d[23])
plt.xlabel('Load or resistance as times the nominal load, $x, E_{mean}(%.2f); x, R_{mean}(%.2f)$' %(d[20],d[22]))
plt.axvline(x=d[20],color='r',linestyle=':')
plt.axvline(x=d[22],color='g',linestyle=':')
plt.savefig('CNS-SNS_L2.png')
pass
def aus():
plt.title('Australian Northwest Shelf')
if cat == 1:
ER(d[24],d[25],d[26],d[27])
plt.xlabel('Load or resistance as times the nominal load, $x, E_{mean}(%.2f); x, R_{mean}(%.2f)$' %(d[24],d[26]))
plt.axvline(x=d[24],color='r',linestyle=':')
plt.axvline(x=d[26],color='g',linestyle=':')
plt.savefig('AUS_L1.png')
if cat == 2:
ER(d[28],d[29],d[30],d[31])
plt.xlabel('Load or resistance as times the nominal load, $x, E_{mean}(%.2f); x, R_{mean}(%.2f)$' %(d[28],d[30]))
plt.axvline(x=d[28],color='r',linestyle=':')
plt.axvline(x=d[30],color='g',linestyle=':')
plt.savefig('AUS_L2.png')
pass
def misc():
plt.rcParams['grid.linestyle'] = ':'
plt.rcParams['grid.linewidth'] = 0.5
plt.grid(True)
plt.xlim(0.25,2.75)
plt.ylim(0,5)
if __name__ == '__main__':
if args['--gom']:
gom()
elif args['--nns']:
nns()
elif args['--cns']:
cns_sns()
elif args['--aus']:
aus()
else:
print "No option was selected. For help, try: python pdf.py -h"