Structural reliability analysis

In this notebook we will explore the intricacies of performing structural reliability analysis.

Computational lab setup



In [42]:

    
% matplotlib inline

import sys

import numpy as np
import scipy as sp
import pandas as pd
import matplotlib as mpl

from scipy.stats import norm, lognorm
from scipy.integrate import quad
import matplotlib.pyplot as plt

print(sys.version)
for module in (np, sp, pd, mpl):
    print('{:.<15}{}'.format(module.__name__, module.__version__))









    



3.5.2 |Anaconda 4.3.0 (32-bit)| (default, Jul  2 2016, 17:49:02) 
[GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]
numpy..........1.11.3
scipy..........0.18.1
pandas.........0.19.2
matplotlib.....2.0.0

Basic reliability problem

Describe problem here!

Failure probability computation

$$PF = \int_{-\infty}^{+\infty}\phi_{demand}(x)\cdot\Phi_{capacity}(x)\;dx = \int_{-\infty}^{+\infty}\left[1-\Phi_{demand}(x)\right]\cdot\phi_{capacity}(x)\;dx$$

First example: normally distributed random variables

Consider a very simple case where both capacity and demand (measured in terms of action-effects) are represented by independent random variables with a normal distribution.



In [43]:

    
capacity = norm(loc=1.5, scale=0.1) # capacity
demand = norm(loc=1.0, scale=0.1) # demand
x = np.linspace(demand.ppf(0.001), capacity.ppf(0.999), 100)

fig, ax = plt.subplots()
ax.plot(x, capacity.pdf(x), label='capacity', color='r')
ax.plot(x, demand.pdf(x), label='demand', color='g')
ax.axvline(capacity.stats(moments='m'), linestyle='--', color='r')
ax.axvline(demand.stats(moments='m'), linestyle='--', color='g')
ax.fill_between(x, 0., np.minimum(capacity.pdf(x),demand.pdf(x)))
ax.set_xlabel('x')
ax.set_ylabel('$\phi(x)$')
ax.legend()









    Out[43]:





<matplotlib.legend.Legend at 0xa9d4d78c>

First numerical solution:

$$PF = \int_{-\infty}^{+\infty}\phi_{demand}(x)\cdot\Phi_{capacity}(x)\;dx$$



In [44]:

    
def f1(x):
    """Objective function for numerical quadrature."""
    return demand.pdf(x)*capacity.cdf(x)


val, err = quad(f1, demand.ppf(0.001), capacity.ppf(0.999))
print('Value = {:1.6e}, error = {:1.6e}'.format(val, err))









    



Value = 2.034760e-04, error = 1.759145e-09



In [45]:

    
fig, ax = plt.subplots()
ax1 = ax.twinx() # secondary axis for capacity
ax.plot(x, demand.pdf(x), label='$\phi_{demand}(x)$', color='g')
ax1.plot(x, capacity.cdf(x), label='$\Phi_{capacity}(x)$', color='r')
ax.axvline(demand.stats(moments='m'), linestyle='--', color='g')
ax1.axvline(capacity.stats(moments='m'), linestyle='--', color='r')
ax1.plot((capacity.stats(moments='m'),x[-1]),(0.5, 0.5), linestyle='--', color='r')
ax.set_xlabel('x')
ax.set_ylabel('$\phi_{demand}(x)$', color='g')
ax.tick_params('y', colors='g')
ax1.set_ylabel('$\Phi_{capacity}(x)$', color='r')
ax1.tick_params('y', colors='r')

Second numerical solution:

$$PF = \int_{-\infty}^{+\infty}\left[1-\Phi_{demand}(x)\right]\cdot\phi_{capacity}(x)\;dx$$



In [46]:

    
def f2(x):
    """Objective function for numerical quadrature."""
#    return (1.-demand.cdf(x))*capacity.pdf(x)
    return demand.sf(x)*capacity.pdf(x) # sf is the survival function


val, err = quad(f2, demand.ppf(0.001), capacity.ppf(0.999))
print('Value = {:1.6e}, error = {:1.6e}'.format(val, err))









    



Value = 2.034760e-04, error = 1.759145e-09



In [47]:

    
fig, ax = plt.subplots()
ax1 = ax.twinx() # secondary axis for capacity
ax.plot(x, 1.-demand.cdf(x), label='$1-\Phi_{demand}(x)$', color='g')
ax1.plot(x, capacity.pdf(x), label='$\phi_{capacity}(x)$', color='r')
ax.axvline(demand.stats(moments='m'), linestyle='--', color='g')
ax1.axvline(capacity.stats(moments='m'), linestyle='--', color='r')
ax.plot((x[0], demand.stats(moments='m')),(0.5, 0.5), linestyle='--', color='g')
ax.set_xlabel('x')
ax.set_ylabel('$1-\Phi_{demand}(x)$', color='g')
ax.tick_params('y', colors='g')
ax1.set_ylabel('$\phi_{capacity}(x)$',color='r')
ax1.tick_params('y', colors='r')

For this first case there is a closed form solution for the failure probability ($PF$):

$$PF = \Phi(\frac{0-\mu}{\sigma}) = \Phi(-\beta)$$

where $\beta=\frac{\mu}{\sigma}$ is called the safety index.

See pages 47-48 of: http://www.km.fgg.uni-lj.si/coste24/data/coimbradocuments/coimbra-faber.pdf



In [48]:

    
miu = capacity.mean()-demand.mean()
sigma = np.sqrt(capacity.std()**2+demand.std()**2)
print('miu = {:1.6e}, sigma = {:1.6e}'.format(miu, sigma))
beta = miu/sigma
print('beta = {:1.6}'.format(beta))
PF = norm.cdf(-beta) # probability of failure
print('Value = {:1.6e}'.format(PF))









    



miu = 5.000000e-01, sigma = 1.414214e-01
beta = 3.53553
Value = 2.034760e-04

Second example: lognormally distributed random variables

Now consider a similar problem but with lorgnormally distributed random variables:



In [49]:

    
capacity2 = lognorm(loc=1.5, s=0.1) # capacity
demand2 = lognorm(loc=1.0, s=0.1) # demand
x = np.linspace(demand2.ppf(0.001), capacity2.ppf(0.999), 100)

fig, ax = plt.subplots()
ax.plot(x, capacity2.pdf(x), label='capacity', color='r')
ax.plot(x, demand2.pdf(x), label='demand', color='g')
ax.axvline(capacity2.stats(moments='m'), linestyle='--', color='r')
ax.axvline(demand2.stats(moments='m'), linestyle='--', color='g')
ax.fill_between(x, 0., np.minimum(capacity2.pdf(x),demand2.pdf(x)))
ax.set_xlabel('x')
ax.set_ylabel('$\phi(x)$')
ax.legend()









    Out[49]:





<matplotlib.legend.Legend at 0xa9b336ac>

First numerical solution:

$$PF = \int_{-\infty}^{+\infty}\phi_{demand}(x)\cdot\Phi_{capacity}(x)\;dx$$



In [50]:

    
def f1(x):
    """Objective function for numerical quadrature."""
    return demand2.pdf(x)*capacity2.cdf(x)


val, err = quad(f1, demand2.ppf(0.001), capacity2.ppf(0.999))
print('Value = {:1.6e}, error = {:1.6e}'.format(val, err))









    



Value = 3.214780e-04, error = 2.015274e-09



In [51]:

    
fig, ax = plt.subplots()
ax1 = ax.twinx() # secondary axis for capacity
ax.plot(x, demand2.pdf(x), label='$\phi_{demand}(x)$', color='g')
ax1.plot(x, capacity2.cdf(x), label='$\Phi_{capacity}(x)$', color='r')
ax.axvline(demand2.stats(moments='m'), linestyle='--', color='g')
ax1.axvline(capacity2.stats(moments='m'), linestyle='--', color='r')
ax1.plot((capacity2.stats(moments='m'),x[-1]),(0.5, 0.5), linestyle='--', color='r')
ax.set_xlabel('x')
ax.set_ylabel('$\phi_{demand}(x)$', color='g')
ax.tick_params('y', colors='g')
ax1.set_ylabel('$\Phi_{capacity}(x)$', color='r')
ax1.tick_params('y', colors='r')

Second numerical solution:

$$PF = \int_{-\infty}^{+\infty}\left[1-\Phi_{demand}(x)\right]\cdot\phi_{capacity}(x)\;dx$$



In [52]:

    
def f2(x):
    """Objective function for numerical quadrature."""
#    return (1.-demand2.cdf(x))*capacity2.pdf(x)
    return demand2.sf(x)*capacity2.pdf(x) # sf is the survival function


val, err = quad(f2, demand2.ppf(0.001), capacity2.ppf(0.999))
print('Value = {:1.6e}, error = {:1.6e}'.format(val, err))









    



Value = 3.214783e-04, error = 1.204997e-09

For this second case there is also a closed form solution for the failure probability ($PF$):

$$PF = \Phi(\frac{0-log\left(\frac{\mu_C}{\mu_D}\sqrt{\frac{1+COV_D^2}{1+COV_C^2}}\right)}{\sqrt{log\left[\left(1+COV_D^2\right)\left(1+COV_C^2\right)\right]}}) = \Phi(-\beta)$$

where $COV=\frac{\sigma}{\mu}$ is the covariance and $\beta=\frac{\mu}{\sigma}$ is called the safety index.

See pages 47-48 of: http://www.km.fgg.uni-lj.si/coste24/data/coimbradocuments/coimbra-faber.pdf



In [53]:

    
cov_C = capacity2.std()/capacity2.mean()
cov_D = demand2.std()/demand2.mean()
miu = np.log(capacity2.mean()/demand2.mean()*np.sqrt((1+cov_D**2)/(1+cov_C**2)))
sigma = np.sqrt(np.log((1+cov_D**2)*(1+cov_C**2)))
print('miu = {:1.6e}, sigma = {:1.6e}'.format(miu, sigma))
beta = miu/sigma
print('beta = {:1.6}'.format(beta))
PF = norm.cdf(-beta) # probability of failure
print('Value = {:1.6e}'.format(PF))









    



miu = 2.230962e-01, sigma = 6.432981e-02
beta = 3.46801
Value = 2.621673e-04

Reliability index

To be done!

Example from:

http://www.springer.com/cda/content/document/cda_downloaddocument/9783319209456-c2.pdf%3FSGWID%3D0-0-45-1544683-p177543534&usd=2&usg=AFQjCNFO1huDK6zcCJeCqWsgSdD6EqkSMg

Definition of COV taken from: https://en.wikipedia.org/wiki/Log-normal_distribution



In [101]:

    
mu_C = 15.
COV_C = 0.2
sigma_C = np.log1p(COV_C)
C = lognorm(sigma_C, loc=0, scale=mu_C) # capacity
print('C_median={}, C_mean={}'.format(C.ppf(0.5), C.stats(moments='m')))

mu_D = 10.
COV_D = 0.1
sigma_D = np.log1p(COV_D)
D = lognorm(sigma_D, loc=0, scale=mu_D) # demand
print('D_median={}, D_mean={}'.format(D.ppf(0.5), D.stats(moments='m')))

x = np.linspace(D.ppf(1e-12), C.ppf(0.9999), 100)
#x = np.linspace(0.01, 30, 1000)

fig, ax = plt.subplots()
ax.plot(x, C.pdf(x), label='$\phi_{C}(x)$', color='r')
ax.plot(x, D.cdf(x), label='$\Phi_{D}(x)$', color='g')
ax.axvline(mu_C, linestyle='--', color='r')
ax.axvline(mu_D, linestyle='--', color='g')
ax.grid(b=True)
ax.legend()
ax.set_xlabel('x')
ax.set_ylabel('$\Phi_D(x)/\phi_C(x)$')


def f1(x):
    """Objective function for numerical quadrature."""
    return D.cdf(x)*C.pdf(x)


def f2(x):
    """Objective function for numerical quadrature."""
    return D.pdf(x)*C.sf(x)


val, err = quad(f1, D.ppf(1e-12), C.ppf(0.9999))
print('Value = {:.3f}, error = {:1.6e}'.format(val, err))
val, err = quad(f2, D.ppf(0.001), C.ppf(0.999))
print('Value = {:.3f}, error = {:1.6e}'.format(val, err))
num = np.log((mu_C/np.sqrt(1.+COV_C**2))/(mu_D/np.sqrt(1.+COV_D**2)))
den = np.sqrt(np.log((1.+COV_D**2)*(1.+COV_C**2)))
R = 1.-norm.cdf(-num/den)
print(R)









    



C_median=15.0, C_mean=15.251391978065396
D_median=10.0, D_mean=10.045523457727839
Value = 0.976, error = 1.090195e-12
Value = 0.975, error = 7.385874e-12
0.961008928246