<jupyter_start><jupyter_text>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<jupyter_code>q_lsf = lambda h: w_max*h**3/c0     # LSF, eq. (3.36) (rewritten), X-space
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()<jupyter_output><empty_output><jupyter_text>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. ![image.png](attachment:image.png) 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<jupyter_code>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
# 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))
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()<jupyter_output><empty_output><jupyter_text>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<jupyter_code>gu = lambda u: w_max*(u[1]*sigma_H + mu_H)**3 - c0*(u[0]*sigma_Q + mu_Q) # eq. (3.41)<jupyter_output><empty_output><jupyter_text>Next, we define the function that gives us the next, improved $\alpha$-vector<jupyter_code>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])<jupyter_output><empty_output><jupyter_text>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$<jupyter_code>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.<jupyter_output><empty_output><jupyter_text>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`<jupyter_code>from scipy.optimize import fsolve<jupyter_output><empty_output><jupyter_text>Since `fsolve` is a numerical solver it needs a starting guess for $\beta$. The logical guess would be the previous $\beta$ we calculated.<jupyter_code>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<jupyter_output><empty_output><jupyter_text>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).<jupyter_code>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
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))<jupyter_output>Reliability index after 1 iteration: 1.08012
Reliability index after 5 iterations: 1.07938<jupyter_text>Lets compare the results in a graph:<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>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
alpha, beta, its = FORM_tol(alpha_start=alpha0,beta_start=beta0, itreturn=True)
print('Reliability index after %i iteration(s): %.5f' % (its, beta))<jupyter_output>Reliability index after 2 iteration(s): 1.07938<jupyter_text>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)$:<jupyter_code>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<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>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<jupyter_output><empty_output><jupyter_text>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) $$<jupyter_code>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<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>Q = {}
Q['type'] = "Gumbel"
Q['moments'] = [mu_Q,sigma_Q]
addparams(Q)

H = {}
H['type'] = "Lognormal"
H['moments'] = [mu_H,sigma_H]
addparams(H)
gX = lambda x: w_max*x[1]**3 - c0*x[0]
gU = lambda u: gX([x2u(Q)(u[0]), x2u(H)(u[1])])
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()<jupyter_output><empty_output><jupyter_text>Limit state function in U-space, isolated for $u_Q$<jupyter_code>uQ_lsf = lambda uH: fsolve(lambda uQ: gU([uQ,uH]), x0=0)
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()<jupyter_output><empty_output><jupyter_text>We now have to update or FORM iteration scheme<jupyter_code>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])
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
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
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)<jupyter_output>Reliability index = 1.15159
Probability of failure = 1.247e-01
{'fraction_non_alphanumeric': 0.15121496416292757, 'fraction_numerical': 0.028611386282850648, 'mean_word_length': 4.235509456985967, 'pattern_counts': {'":': 3, '<': 62, '<?xml version=': 0, '>': 64, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '52431', 'license': 'MIT', 'n_tokens_mistral': 6692, 'n_tokens_neox': 6492, 'n_words': 1769, 'path': 'E3_1.ipynb', 'repo_name': 'jochenkohler/TKT4196_material', 'size': '17161'}