```python
import numpy as np
import scipy as sp
from scipy.stats import multivariate_normal as mvn
import matplotlib.pyplot as plt
#from matplotlib.figure import figaspect
import time
import sympy as sym
fontsizes=18
plt.rcParams.update({'font.size': fontsizes})
plt.rcParams.update({"font.family": "serif"})
plt.rcParams.update({"mathtext.fontset" : "cm"})
plt.rcParams.update({'font.serif': 'Times New Roman'})
plt.close('all')
```


```python
# =============================================================================
# INPUT

b = 50                          # [mm] Cross-section width
E = 205000                      # [MPa] Elastic modulus
l = 3000                        # [mm] Beam span
mu_H = 100                      # [mm] Cross-section height mean
sigma_H = 5                     # [mm] Cross-section height std. dev.
mu_Q = 5                        # [N/mm] Mean load
sigma_Q = 1                     # [N/mm] Std. load
w_max = 8                       # [mm] 5/384*ql^4/EI

c0 = 5*l**4/(384*E*b/12)
```


```python
# Plot limit state function
# plot x-space

MU = np.array([mu_Q,mu_H])
CX = np.array([[sigma_Q**2, 0],           #Covariation Matrix
              [0, sigma_H**2]])

fx = lambda X: mvn(MU,CX).pdf(X)          #Probability density function

dum_Q = np.linspace(0,mu_Q+5*sigma_Q,100) #Dummy variables (only used for plotting)
dum_H = np.linspace(0,mu_H+10*sigma_H,100)
X1, X2 = np.meshgrid(dum_Q,dum_H, indexing='ij')
X = np.dstack((X1, X2))
```


```python
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X1,X2,fx(X),cmap=plt.cm.coolwarm)
ax.set_xlabel(r'$q$', labelpad=10)
ax.set_ylabel('$h$', labelpad=10)
ax.set_zlabel(r'$f_{\mathbf{X}}(q,h)$', labelpad=20, rotation=90)
ax.set_title('Probability density function (PDF)')
fig.tight_layout()
plt.show()
```

Now we want to add the limit state function to the plot. Since the plot can look a little crowded in 3D we move over to a contour plot. To be able to plot the limit state function as a curve we must rewrite it so that one of the random variables is a function of the other. (Think $y=f(x)$). We get the nicest expression when we write $Q$ as a function of $H$, but the other way around is also viable. Be sure to be consistent with what goes along the x- and y-axis when plotting if opting for the latter.

$$w_{max} h^3 - c_0 q = 0$$
$$\Rightarrow q(h) = \frac{w_{max}h^3}{c_0}$$

Limit state function in the X-space


```python
q_lsf = lambda h: w_max*h**3/c0     # LSF, eq. (3.36) (rewritten), X-space
```


```python
fig = plt.figure(figsize=(8,8))
#Plot H in 2D
axH = fig.add_subplot(224)
axH.set_box_aspect(2)
fH = lambda h: sp.stats.norm.pdf(h,loc=mu_H,scale=sigma_H)
axH.plot(fH(dum_H),dum_H, 'k')
axH.set_ylim(0,dum_H[-1])
axH.axis('off')

#plot Q in 2D
axQ = fig.add_subplot(221)
axQ.set_box_aspect(0.5)
fQ = lambda q: sp.stats.norm.pdf(q,loc=mu_Q,scale=sigma_Q)
axQ.plot(dum_Q,fQ(dum_Q),'k')
axQ.set_xlim(0,dum_Q[-1])
axQ.axis('off')

#Plot Q, H and LSF together in contour plot
ax1 = fig.add_subplot(223)
ax1.set_box_aspect(1)
ax1.contour(X1, X2, fx(X))
ax1.plot(q_lsf(dum_H),dum_H);
ax1.set_xlim(0,dum_Q[-1])
ax1.set_ylim(0,dum_H[-1])
ax1.set_title('X-space')
ax1.set_xlabel('$q$')
ax1.set_ylabel('$h$')
plt.tight_layout()
fig.subplots_adjust(wspace=-0.2,hspace=-0.2)
plt.show()
```

To simplify calculations greatly, and take use of the standard normal distribution's useful properties we should first of all transform our coordinate system, i.e. map our space to the so-called u-space. The procedure is illustrated below. One can observe that there is a 1-1 correspondance between $x$ and $u$ when going through the cumulative distribution functions. Also, note that the y-axis in the two plots are equal, both dimensionless and spanning from 0 to 1 as they represent probabilities.



The procedure can be represented by three steps:
- Find the y-value in the second graph, depending on the u-value: $y = \Phi(u_i)$
- Transfer the y-value to the first graph
- Find the x-value, depending on y-value: $x_i = F_{x_i}^{-1}(y)$

In total we get:

$$x_i = F_{X_i}^{-1}(\Phi(u_i))$$

For a normal distributed variable we have that 

$$x_i = F_{x_i}^{-1}(\Phi(u_i)) = u_i \sigma_i + \mu_i$$

Therefore:

$$ h = u_H \sigma_H + \mu_H $$
$$ q = u_Q \sigma_Q + \mu_Q $$

Substituing this into the limit state function yields

$$w_{max}(u_H \sigma_H + \mu_H)^3 - c0(u_Q \sigma_Q + \mu_Q) = 0$$

Once again, for plotting, we want to write one variable as a function of the other. Isolating $u_Q$ gives

$$ u_Q(u_H) = \frac{1}{\sigma_Q}\left(\frac{w_{max}(u_H \sigma_H + \mu_H)^3}{c_0} - \mu_Q\right)$$

Limit state function in U-space


```python
uQ_lsf = lambda uH: 1/(sigma_Q)*(w_max*(uH*sigma_H + mu_H)**3/c0 - mu_Q) #LSF, eq. (3.41) (rewritten), U-space
```


```python
# Plot limit state function
# plot u-space

MU_U = np.array([0,0])
CU = np.array([[1, 0],
              [0, 1]])

fu = lambda U: mvn(MU_U,CU).pdf(U)

dum_U = np.linspace(-5,5,100)
U1, U2 = np.meshgrid(dum_U,dum_U, indexing='ij')
U = np.dstack((U1, U2))
```


```python
fig = plt.figure(figsize=(8,8))

#Plot UH in 2D
axUH = fig.add_subplot(224)
axUH.set_box_aspect(2)
axUH.plot(sp.stats.norm.pdf(dum_U),dum_U, 'k')
axUH.axis('off')

#Plot UQ in 2D
axUQ = fig.add_subplot(221)
axUQ.set_box_aspect(0.5)
axUQ.plot(dum_U,sp.stats.norm.pdf(dum_U), 'k')
axUQ.axis('off')

ax1 = fig.add_subplot(223)
ax1.contour(U1, U2, fu(U))
ax1.plot(uQ_lsf(dum_U),dum_U);
ax1.set_xlim(dum_U[0],dum_U[-1])
ax1.set_ylim(dum_U[0],dum_U[-1])
ax1.set_title('U-space')
ax1.set_xlabel('$u_Q$')
ax1.set_ylabel('$u_H$')
ax1.set_aspect('equal')
plt.tight_layout()
fig.subplots_adjust(wspace=-0.2,hspace=-0.2)
plt.show()
```

Note that the contours now are symmetric, moreover symmetric around (0,0). This is of course because the standard normal distribution is symmetric along the y-axis (x=0).

Now we want to find the reliability index graphically by finding the shortest path from the origin to the LSF. To do this we adapt the iteration scheme from the compendium, eqs. (3.10) and (3.11). First we define the limit state function. This can be done by defining $g_u = g_u(u_H, u_Q)$, however we do it more compact by introducing $\mathbf{u}=[u_Q, u_H]$, therefore

$$g_u = g_u(\mathbf{u})$$

where

$$\mathbf{u}[0] = u_Q \: , \:\mathbf{u}[1] = u_H$$

We are now ready to define our limit state function, in the U-space


```python
gu = lambda u: w_max*(u[1]*sigma_H + mu_H)**3 - c0*(u[0]*sigma_Q + mu_Q) # eq. (3.41)
```

Next, we define the function that gives us the next, improved $\alpha$-vector


```python
def alpha_next(u):
    dgdu0 = -c0*sigma_Q                                 #dgu/du[0] = dgdu/duQ
    dgdu1 = 3*w_max*(u[1]*sigma_H + mu_H)**2*sigma_H    #dgu/du[1] = dgu/duH
    k = np.sqrt(dgdu0**2+dgdu1**2)                      #Normalization factor, makes ||alpha|| = 1
    return np.array([-dgdu0/k,-dgdu1/k])
```

Note that this does not correspond exactly to eq. (3.10) as this equation requires us to insert the previous values of $\beta$ and $\mathrm{\alpha}$, in other words, a new  $\mathrm{\alpha}$ is not found until we do this. This means that we need some starting values, $\beta_0$ and $\mathrm{\alpha}_0$


```python
beta0 = 2
alpha0 = np.array([1/np.sqrt(2),-1/np.sqrt(2)]) #we need ||alpha|| = 1, thereby the sqrt()'s. Starting alpha-value should be positive for load variables and negative for resistance variables.
```

Now we need a way to solve $g_u(\beta \cdot \mathbf{\alpha}) = 0$ with respect to $\beta$. This can be done in Python by utilizing `fsolve`


```python
from scipy.optimize import fsolve
```

Since `fsolve` is a numerical solver it needs a starting guess for $\beta$. The logical guess would be the previous $\beta$ we calculated.


```python
def beta_next(alpha, beta_prev):
    equation = lambda beta: gu(beta*alpha)
    beta_new = fsolve(equation, x0=beta_prev) #Finds beta that gives equtaion = 0
    return beta_new
```

Now, lets put them togheter to make the iteration scheme complete. First we want a function that takes in starting guesses, as well as number of iterations we want. We want it to return the final calculated values (or a list of all calculated values).


```python
def FORM_it(alpha_start, beta_start, it):
    alphas = [alpha_start]
    betas = [beta_start]
    for i in range(0,it):
        alpha_new = alpha_next(betas[-1]*alphas[-1]) #calculates new, improved alpha, using the last calculated values (by using [-1])
        alphas.append(alpha_new)                     #adds the new, improved alpha to the list of alphas
        beta_new = beta_next(alphas[-1],betas[-1])   #calculates new, improved beta
        betas.append(beta_new)                       #adds the new, improved beta to the list of betas
    #return alphas, betas                             #returns the whole list, letting us see the iteration progress better
    return alphas[-1], betas[-1]                     #returns only the final values
```


```python
alpha1, beta1 = FORM_it(alpha0, beta0, it=1)
alpha5, beta5 = FORM_it(alpha0, beta0, it=5)
print('Reliability index after 1 iteration: %.5f' % (beta1))
print('Reliability index after 5 iterations: %.5f' % (beta5))
```

    Reliability index after 1 iteration: 1.08012
    Reliability index after 5 iterations: 1.07938


Lets compare the results in a graph:


```python
uQ0, uH0 = beta0*alpha0
uQ1, uH1 = beta1*alpha1
uQ5, uH5 = beta5*alpha5

fig = plt.figure(figsize=(8,8))

ax1 = fig.add_subplot(111)
ax1.contour(U1, U2, fu(U))
ax1.plot(uQ_lsf(dum_U),dum_U,'k')
ax1.plot(uQ0,uH0,'mo')
ax1.plot([0,uQ0],[0,uH0],'m', label='0 iterations')
ax1.plot(uQ1,uH1,'bo')
ax1.plot([0,uQ1],[0,uH1],'b', label='1 iteration')
ax1.plot(uQ5,uH5,'ro')
ax1.plot([0,uQ5],[0,uH5],'r', label='5 iterations')
ax1.set_xlim(-2,2)
ax1.set_ylim(-2,2)
ax1.set_title('U-space')
ax1.set_xlabel('$u_Q$')
ax1.set_ylabel('$u_H$')
ax1.set_aspect('equal')
ax1.legend(loc = 'upper left', bbox_to_anchor=(1, 1))
plt.tight_layout()
plt.show()
```

The natural question to now ask yourself is "how many iterations is enough?". We do not know this beforehand, and is dependent on the initial guess and of course the accuracy we want. We should modify our iteration scheme as to take in a error tolerance rather than number of iterations.

This way, the algorithm iterates until the error between the last and second to last calculation is small enough.


```python
def FORM_tol(alpha_start, beta_start, tol=1e-5, itreturn=False):
    alphas = [alpha_start]
    beta1 = beta_next(alpha_start, beta0)
    betas = [beta_start, beta1]                      #Need at least two values of beta for the next line to work, first step is therefore done manually
    while abs((betas[-1] - betas[-2])/betas[-1]) >= tol:
        alpha_new = alpha_next(betas[-1]*alphas[-1]) #calculates new, improved alpha, using the last calculated values (by using [-1])
        alphas.append(alpha_new)                     #adds the new, improved alpha to the list of alphas
        beta_new = beta_next(alphas[-1],betas[-1])   #calculates new, improved beta
        betas.append(beta_new)                       #adds the new, improved beta to the list of betas
    if itreturn:
        return alphas[-1], betas[-1], len(alphas)-1  #returns the final values, and number of iterations
    else:
        return alphas[-1], betas[-1]                 #returns only the final values
```


```python
alpha, beta, its = FORM_tol(alpha_start=alpha0,beta_start=beta0, itreturn=True)
print('Reliability index after %i iteration(s): %.5f' % (its, beta))
```

    Reliability index after 2 iteration(s): 1.07938


We see that 5 iterations was abundant, only 2 was needed to achieve sufficient accuracy.

# General procedure

Earlier, you saw that we could define the limit state function as a function of $u$ directly. However, when working with random variables that are not normal distributed, this expression can become quite messy. Instead we can utilize mapping, and the chain rule.

Recall that we turn our limit state function from a function of $x$, to a function of $u$

$$g=g(x)=g(x(u)) \Rightarrow g=g(u)$$

To do this mapping, we need the to know the function $x(u)$

We recall that

$$x = F_{X}^{-1}(\Phi(u))$$

Note that that $F_{x_i}$ depends on the distribution type of variable $X$, whilst $\Phi$ is always the same.


The compendium covers the distribution types which we represent variables with in this course, we program their $x(u)$:


```python
def x2u(X):
    #X is here a dictionary, which describes the type and stores its parameters
    if X['type'] == 'Normal':
        # If X is normal distributed, use Eq. (3.8) (isolated x) 
        x = lambda u: u*X['theta'][1] + X['theta'][0]
        # X['theta'][1] = sigma
        # X['theta'][0] = mu
    if X['type'] == 'Lognormal':
        # If X is lognormal distributed, use Eq. (3.26)
        x = lambda u: np.exp(X['theta'][1]*u+X['theta'][0])
        # X['theta'][1] = sigma_L
        # X['theta'][0] = mu_L
    if X['type'] == 'Gumbel':
        # If X is Gumbel distributed, use Eq. (3.27)
        x = lambda u: X['theta'][1]-1/X['theta'][0]*np.log(-np.log(sp.stats.norm.cdf(u)))
        # X['theta'][1] = b
        # X['theta'][0] = a
    return x
```

We are often given moments (mean and standard deviation), rather than the parameters. However, we can calculate the parameters given the moments. We program functions for these calculations.


```python
def addparams(X):
    moments = X['moments']
    if X['type'] == 'Normal':
        theta = X['moments']
    if X['type'] == 'Lognormal':
        sigma_ln = np.sqrt(np.log(moments[1]**2/moments[0]**2+1))
        mu_ln = np.log(moments[0])-0.5*sigma_ln**2
        theta = [mu_ln,sigma_ln]
    if X['type'] == 'Gumbel':
        a = np.pi/np.sqrt(6)/moments[1]
        b = moments[0] - 0.5772/a
        theta = [a,b]
    X['theta']=theta
        
```

We want to differentiate the limit state with respects to $u$ in our iteration. We can do this by utilizing the chain rule:

$$\frac{\partial g}{\partial u} = \frac{\partial g}{ \partial x}\frac{\partial x}{\partial u}$$

Note that $X$'s distribution only affects $\frac{\partial x}{\partial u}$, not $\frac{\partial g}{ \partial x}$. We will now make a function that calculates $\frac{\partial x}{\partial u}$, given $X$'s distribution. Recall that

$$x = x(u) = F_{x}^{-1}(\Phi(u))$$

therefore,

$$ \frac{\partial x}{\partial u} = \frac{\partial}{\partial u} \left( F_{x}^{-1}(\Phi(u) \right) $$


```python
def dxdufun(X):
    if X['type'] == "Normal":
        dxdu = lambda u: X['theta'][1]
    if X['type'] == "Lognormal":
        dxdu = lambda u: X['theta'][1] * np.exp(X['theta'][0]+X['theta'][1] * u)
    if X['type'] == "Gumbel":
        dxdu = lambda u: - sp.stats.norm.pdf(u)/(X['theta'][0] * sp.stats.norm.cdf(u)*np.log(sp.stats.norm.cdf(u)))
    return dxdu
```

# Example cont.

We can now perform FORM with non-normal distributed variables. We try it out on our example by now assuming $Q$ and $H$ to be gumbel and lognormal distributed, respectively.


```python
Q = {}
Q['type'] = "Gumbel"
Q['moments'] = [mu_Q,sigma_Q]
addparams(Q)

H = {}
H['type'] = "Lognormal"
H['moments'] = [mu_H,sigma_H]
addparams(H)
```


```python
gX = lambda x: w_max*x[1]**3 - c0*x[0]
gU = lambda u: gX([x2u(Q)(u[0]), x2u(H)(u[1])])
```


```python
fH = lambda h: sp.stats.lognorm.pdf(h, s = H['theta'][1], scale = np.exp(H['theta'][0]))
fQ = lambda q: sp.stats.gumbel_r.pdf(q, loc = Q['theta'][1], scale = 1/Q['theta'][0])

fX = lambda X1,X2: fQ(X1)*fH(X2)

fig = plt.figure(figsize=(8,8))

#Plot H in 2D
axH = fig.add_subplot(224)
axH.set_box_aspect(2)
axH.plot(fH(dum_H),dum_H, 'k')
axH.set_ylim(0,dum_H[-1])
axH.axis('off')

#plot Q in 2D
axQ = fig.add_subplot(221)
axQ.set_box_aspect(0.5)
axQ.plot(dum_Q,fQ(dum_Q),'k')
axQ.set_xlim(0,dum_Q[-1])
axQ.axis('off')


ax1 = fig.add_subplot(223)
ax1.plot(q_lsf(dum_H),dum_H)
ax1.contour(X1, X2, fX(X1,X2))
ax1.set_xlim(0,dum_Q[-1])
ax1.set_ylim(0,dum_H[-1])
ax1.set_title('X-space')
ax1.set_xlabel('$q$')
ax1.set_ylabel('$h$')
plt.tight_layout()
fig.subplots_adjust(wspace=-0.2,hspace=-0.2)
plt.show()
```

Limit state function in U-space, isolated for $u_Q$


```python
uQ_lsf = lambda uH: fsolve(lambda uQ: gU([uQ,uH]), x0=0)
```


```python
fig = plt.figure(figsize=(8,8))

#Plot UH in 2D
axUH = fig.add_subplot(224)
axUH.set_box_aspect(2)
axUH.plot(sp.stats.norm.pdf(dum_U),dum_U, 'k')
axUH.axis('off')

#Plot UQ in 2D
axUQ = fig.add_subplot(221)
axUQ.set_box_aspect(0.5)
axUQ.plot(dum_U,sp.stats.norm.pdf(dum_U), 'k')
axUQ.axis('off')

ax1 = fig.add_subplot(223)
ax1.contour(U1, U2, fu(U))
ax1.plot(np.vectorize(uQ_lsf)(dum_U), dum_U)
ax1.set_xlim(dum_U[0],dum_U[-1])
ax1.set_ylim(dum_U[0],dum_U[-1])
ax1.set_title('U-space')
ax1.set_xlabel('$u_Q$')
ax1.set_ylabel('$u_H$')
ax1.set_aspect('equal')
plt.tight_layout()
fig.subplots_adjust(wspace=-0.2,hspace=-0.2)
plt.show()
```

We now have to update or FORM iteration scheme


```python
def alpha_next(u):
    dgdu0 = -c0 * dxdufun(Q)(u[0])
    dgdu1 = w_max * 3 * (x2u(H)(u[1]))**2 * dxdufun(H)(u[1])
    k = np.sqrt(dgdu0**2+dgdu1**2)
    return np.array([-dgdu0/k,-dgdu1/k])
```


```python
def beta_next(alpha, beta_prev):
    equation = lambda beta: gU(beta*alpha)
    beta_new = fsolve(equation, x0=beta_prev) #Finds beta that gives equtaion = 0
    return beta_new
```


```python
def FORM_tol(alpha_start, beta_start, tol=1e-5, itreturn=False):
    alphas = [alpha_start]
    beta1 = beta_next(alpha_start, beta0)
    betas = [beta_start, beta1]                      #Need at least two values of beta for the next line to work, first beta is therefore calculated manually
    while abs((betas[-1] - betas[-2])/betas[-1]) >= tol:
        alpha_new = alpha_next(betas[-1]*alphas[-1]) #calculates new, improved alpha, using the last calculated values (by using [-1])
        alphas.append(alpha_new)                     #adds the new, improved alpha to the list of alphas
        beta_new = beta_next(alphas[-1],betas[-1])   #calculates new, improved beta
        betas.append(beta_new)                       #adds the new, improved beta to the list of betas
    if itreturn:
        return alphas[-1], betas[-1], len(alphas)-1  #returns the final values, and number of iterations
    else:
        return alphas[-1], betas[-1]                 #returns only the final values
```


```python
alpha, beta = FORM_tol(alpha_start=alpha0, beta_start=beta0)
pf = sp.stats.norm.cdf(-beta)
print('Reliability index = %.5f' % beta)
print('Probability of failure = %.3e' % pf)
```

    Reliability index = 1.15159
    Probability of failure = 1.247e-01


Test the mapping yourselves, what happens if you change $Q$'s and $H$'s type to "Normal", run that cell + the last cell? Is the result what you expect?


```python

```

5096009c54dbd28be1af8847a9c3432711fe9c91
ipynb
Jupyter Notebook
E3_1.ipynb
jochenkohler/TKT4196_material
2e3be73beaafd22387ea27cd8d49f04a776e3bce
MIT
E3_1.ipynb
jochenkohler/TKT4196_material
2e3be73beaafd22387ea27cd8d49f04a776e3bce
MIT
E3_1.ipynb
jochenkohler/TKT4196_material
2e3be73beaafd22387ea27cd8d49f04a776e3bce
MIT
Qwen/Qwen-72B
1. YES
2. YES