Perspective Chapter: Probabilistic Modeling of Failure – Nonlinear Approximation

Abstract

The failure of systems to meet the specified requirements may have adverse effects on their integrity and reliability. The systems could be mechanical, electrical, structural, telecommunications, or electronic that are designed and built to satisfy certain technical specifications and operational requirements. Failure does not necessarily mean the occurrence of a disaster or damage to the system, but also the degraded performance of such systems is considered a failure. One of the essential indicators of the performance and reliability of a system is the probability of failure which is computed by probabilistic methods. One of these methods is the first-order reliability method (FORM). Using FORM to estimate the probability of failure of systems having a nonlinear or a higher-order performance function may provide inaccurate results that may lead to misleading conclusions. To resolve this issue, the second-order reliability method (SORM) is recommended to estimate the probability of failure. This chapter presents commonly used probabilistic approximation methods to estimate the probability of failure for nonlinear performance functions. Illustrative examples to demonstrate the application of these methods are provided at the end of the chapter.

1. Introduction

Probability of failure can be calculated using analytical as well as simulation methods. Simulation methods such as the Monte Carlo technique can be used to calculate the probability of failure. It is an essential step in the calculation process to validate the analytical results. One of the analytical methods to compute the probability of failure is the first-order reliability method (FORM) which is based on the first-order expansion of the Taylor series. This method may provide inaccurate results when approximating nonlinear or higher-order performance functions. It should be noted that the performance function can be referred to as the limit state function (LSF).

The second-order reliability method, abbreviated as SORM, is used to resolve the nonlinearity issues of the performance function or the LSF. It uses the second-order expansion of the Taylor series to include the curvature of the LSF in the calculation to achieve better accuracy of the results.

Figure 1 shows the failure regions for linear and nonlinear LSFs. It also shows the most probable point of failure (MPPF) as the tangent point on the limit-state surface, $Z\left(U^{\ast }\right)=0$, and the reliability index, $\beta$, is the shortest distance from the origin to the limit-state surface.

Figure 1.

Representation of the linear and nonlinear LSF.

Breitung [1], Hohenbichler and Rackwitz [2], Tvedt [3, 4], and many other researchers came up with methods to calculate the probability of failure for nonlinear second-order LSFs. These methods use $\beta$ obtained by FORM and the principal curvature, $k_i$, of the limit state surface in the calculation. Essentially, $k_i$ estimates how much the curve diverges from the straight line, and it is an essential part of the calculation process. It is added to the calculation to provide a more accurate estimate of the probability of failure obtained by FORM. This chapter presents commonly used techniques that provide solutions to nonlinear second-order LSFs. The focus of this chapter is on SORM with the assumption that the random variables of the LSF are uncorrelated. A brief discussion about SORM and the commonly used methods for calculating the probability of failure are provided in the next section. Numerical examples are provided in Section 3 to demonstrate the calculation steps and provide in-depth understanding of these methods. Appendix A provides an overview of FORM and the steps involved in its calculation.

2. Second Order Reliability Method (SORM)

The LSF may exhibit nonlinearity when random variables have a non-normal distribution or a nonlinear relationship. Moreover, transforming uncorrelated variables to correlated variables may present nonlinearity in the LSF. SORM is used to deal with the nonlinearity issues to provide more accurate results. It uses the second-order Taylor series expansion of the LSF at the most probable point of failure (MPPF). The vectorized form of the second-order Taylor Series expansion is used to express the LSF around the MPPF as:

$$Z\left(U\right)\approx Z\left(U^{\ast }\right)+\nabla Z{\left(U^{\ast }\right)}^t\left(U-U^{\ast }\right)+0.5{\left(U-U^{\ast }\right)}^t{\nabla }^{2}Z\left(U^{\ast }\right)\left(U-U^{\ast }\right)$$E1

where:

• $Z\left(U\right)$: LSF.

• $U$: represents the vector of the input variables.

• $U^{\ast }$: is the MPPF.

• $Z\left(U^{\ast }\right):$ constant term.

• $\nabla Z\left(U^{\ast }\right)$: is the first-order partial derivative of the LSF evaluated at MPPF.

• $\nabla Z{\left(U^{\ast }\right)}^t\left(U-U^{\ast }\right):$ is the linear term, that is, the first-order derivative.

• 0.5${\left(U-U^{\ast }\right)}^t{\nabla }^{2}Z\left(U^{\ast }\right)\left(U-U^{\ast }\right)$: is the quadratic term.

• ${\nabla }^{2}Z\left(U^{\ast }\right):$ is a square matrix of the second-order partial derivative of the LSF evaluated at MPPF. It is called the Hessian matrix.

Since $U^{\ast }$ is located on the limit-state surface, $Z\left(U^{\ast }\right)$ equals zero.

$Z\left(U\right)$ is normalized by dividing Eq. (1) by $\left|\nabla Z\left(U^{\ast }\right)\right|\text{which yields}:$

$$Z_{\mathit{norm}}\left(U\right)\approx \frac{\nabla Z{\left(U^{\ast }\right)}^t\left(U-U^{\ast }\right)}{\left|\nabla Z\left(U^{\ast }\right)\right|}+\frac{0.5{\left(U-U^{\ast }\right)}^t{\nabla }^{2}Z\left(U^{\ast }\right)\left(U-U^{\ast }\right)}{\left|\nabla Z\left(U^{\ast }\right)\right|}$$E2

where $\left|\nabla Z\left(U^{\ast }\right)\right|$ is the length of the gradient vector of the LSF at the MPPF and is equal to:

$$\left|\nabla Z\left(U^{\ast }\right)\right|=\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(U^{\ast }\right)}{\partial u_i}\right)}^2}$$E3

The approximation of LSF at the MPPF is expressed as:

$$Z_{\mathit{norm}}\left(U\right)\approx \frac{\nabla Z{\left(U^{\ast }\right)}^t\left(U-U^{\ast }\right)}{\left|\nabla Z\left(U^{\ast }\right)\right|}+\frac{0.5{\left(U-U^{\ast }\right)}^t{\nabla }^{2}Z\left(U^{\ast }\right)\left(U-U^{\ast }\right)}{\left|\nabla Z\left(U^{\ast }\right)\right|}=0$$E4

From the first-order reliability approximation (FORM approximation) we have:

$$\frac{\nabla Z{\left(U^{\ast }\right)}^t{U}^{\ast }}{\left|\nabla Z\left(U^{\ast }\right)\right|}=-\beta$$E5

$$\frac{\nabla Z\left(U^{\ast }\right)}{\left|\nabla Z\left(U^{\ast }\right)\right|}=\mathit{\alpha }$$E6

$\alpha$: the directional cosines, ${\alpha }_1,\dots ,{\alpha }_n$ the components of the unit gradient vector $\alpha .$ These are the cosines of the angles between the vector $\beta$ and the axes.

Denoting $\frac{{\nabla }^{2}Z\left(U^{\ast }\right)}{\left|\nabla Z\left(U^{\ast }\right)\right|}$ as H:

$$\mathit{H}\mathbf{=}\frac{{\nabla }^{2}Z\left(U^{\ast }\right)}{\left|\nabla Z\left(U^{\ast }\right)\right|}$$E7

Using Eqs. (5) and (7), Eq. (4) yields:

$$Z_{\mathit{norm}}\left(U\right)\approx \frac{\nabla Z{\left(U^{\ast }\right)}^{t}U}{\left|\nabla Z\left(U^{\ast }\right)\right|}+\beta +0.5{\left(U-U^{\ast }\right)}^t\mathit{H}\left(U-U^{\ast }\right)=0$$E8

To go further with the solution of SORM, the $U$ variables are transformed into a new set of standard normal random variables denoted as $U_i^{\prime }$ where the axis of the last variable coincides with the $\beta$ vector as indicated in Figure 2 [5, 6, 7].

Figure 2.

Second-order failure surface shown in the rotated coordinates.

This is an orthogonal transformation where:

$${\mathit{U}}^{\prime }=\mathit{RU}$$E9

$\mathit{R}$: orthogonal n by n rotation matrix, where n is the number of variables. Since $\mathit{R}$ is an orthogonal rotation matrix, ${\mathit{R}}^{\mathit{-}\mathbf{1}}\mathit{=}{\mathit{R}}^{\mathit{t}}$

$$\mathit{U}={\mathit{R}}^{\mathit{t}}{\mathit{U}}^{\prime }$$E10

The matrix $\mathit{R}$ is constructed in two steps. First, define an initial matrix, ${\mathit{R}}_{\mathbf{01}}$ that consists of rows representing the unit vectors of the axes of the input variables. Since, we want to make $\beta$ vector coincide with the axis of the last variable, we can simply substitute the last row of matrix ${\mathit{R}}_{\mathbf{01}}$ with the directional cosines of $\beta$. The following row vectors are the components of the ${\mathit{R}}_{\mathbf{01}}$ matrix:

$$\begin{array}{l}{r}_{01}=\left[\begin{array}{ccccccccccccc}1& 0& 0& .& .& .& .& .& .& .& .& .& 0\end{array}\right]\\ r_{02}=\left[\begin{array}{ccccccccccccc}0& 1& 0& .& .& .& .& .& .& .& .& .& 0\end{array}\right]\\ \begin{array}{c}\begin{array}{c}.\\ .\\ .\end{array}\\ r_{0n}=\left[\begin{array}{cccc}-\frac{\raisebox{1ex}{${\displaystyle \partial Z\left(U^{\ast }\right)}$}\!\left/ \!\raisebox{-1ex}{${\partial U}_1$}\right.}{\left|\nabla Z\left(U^{\ast }\right)\right|}& .& .& .-\frac{\raisebox{1ex}{${\displaystyle \partial Z\left(U^{\ast }\right)}$}\!\left/ \!\raisebox{-1ex}{${\partial U}_n$}\right.}{\left|\nabla Z\left(U^{\ast }\right)\right|}\end{array}\right]\end{array}\end{array}$$E11

The Gram-Schmidt process is applied to orthogonalize the matrix [8]. The next step is to orthonormalize the matrix by dividing each element of the matrix by the length of its corresponding row vector.

Perform the Gram-Schmidt process for the matrix ${\mathbf{R}}_{\mathbf{01}}$ using Eqs. (12) and (13), in reverse order:

$$r_n=r_{0n}$$E12

$$r_i=r_{0i}-\sum _{j=i+1}^n\frac{r_j{r}_{0i}^t}{r_j{r}_j^t}{r}_{j}i=1,2,\dots .,n-1$$E13

Then orthonormalization is performed for each vector to produce the matrix $\mathit{R}$. Using Eq. (10), Eq. (8) yields:

$$Z_{\mathit{norm}}\left(U^{\prime }\right)\approx -U_n^{\prime }+\beta +0.5{\left(U^{\prime }-{U^{\prime }}^{\ast }\right)}^t\mathit{RH}{\mathit{R}}^{\mathit{t}}\left(U^{\prime }-{U^{\prime }}^{\ast }\right)=0$$E14

The $\mathit{H}$ matrix is the second-order derivative of the LSF at the design point $U^{\ast }$ divided by the length of the gradient vector, $\left|\nabla Z\left(U^{\ast }\right)\right|$, and is formulated as:

$$\mathit{H}\mathbf{=}\frac{1}{\left|\nabla Z\left({U^{\prime }}^{\ast }\right)\right|}\left[\begin{array}{cccccc}\frac{{\partial }^{2}Z\left({U^{\prime }}^{\ast }\right)}{{\partial u}_1^2}& \frac{{\partial }^{2}Z\left({U^{\prime }}^{\ast }\right)}{\partial u_1\partial u_2}& .& .& .& \frac{{\partial }^{2}Z\left({U^{\prime }}^{\ast }\right)}{\partial u_1\partial u_n}\\ \frac{{\partial }^{2}Z\left({U^{\prime }}^{\ast }\right)}{\partial u_2\partial u_1}& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ \frac{{\partial }^{2}Z\left({U^{\prime }}^{\ast }\right)}{\partial u_n\partial u_1}& .& .& .& .& \frac{{\partial }^{2}Z\left({U^{\prime }}^{\ast }\right)}{{\partial u}_n^2}\end{array}\right]$$E15

Denoting $\mathit{RH}{\mathit{R}}^{\mathit{t}}$ as $\mathit{B}$, Eq. (14) becomes:

$$Z_{\mathit{norm}}\left(U^{\prime }\right)\mathbf{\approx }\mathbf{-}{U}_n^{\prime }\mathbf{+}\beta \mathbf{+}0.5\left({{\mathit{U}}^{\prime }}^{\mathit{t}}\mathit{B}{\mathit{U}}^{\prime }\right)\mathbf{=}0$$E16

Then the solution of the second-order approximation is given as:

$$U_n^{\prime }=\beta +0.5\left({{\mathit{U}}^{\prime }}^t\mathit{B}{\mathit{U}}^{\prime }\right)$$E17

$\mathit{B}$ is an (n − 1) by (n − 1) matrix and ${\mathit{U}}^{\prime }$ is a 1 by (n − 1) vector. The eigenvalues of matrix $\mathit{B}$ are computed to obtain the main curvatures, ${k_i}^{\prime }s$ , of the LSF at the MPPF. These ${k_i}^{\prime }s$ are used in the calculation of the probability of failure for nonlinear LSF as indicated in the next subsection [5, 6, 7]. The $0.5\left({{\mathit{U}}^{\prime }}^t\mathit{B}{\mathit{U}}^{\prime }\right)$ term is simplified to $0.5\sum _{i=1}^{n-1}{k}_i{U}_i^{\prime 2}$. Thus, Eq. (17) is expressed in terms of the curvatures, ${k_i}^{\prime }s$, as:

$$U_n^{\prime }=\beta +0.5\sum _{i=1}^{n-1}{k}_i{U}_i^{\prime 2}$$E18

If $k_i>0$, the LSF will have a convex failure region, and if $k_i<0$, it will have a concave failure region as shown in Figure 3. The figure shows the failure surfaces for an arbitrary LSF with three random variables.

Figure 3.

Convex and concave failure regions for an arbitrary LSF with three random variables.

The calculation of $\nabla Z\left(U^{\ast }\right)$ and ${\nabla }^{2}Z\left(U^{\ast }\right)$ can be performed in the X domain using the following equations:

$$\nabla Z\left({U^{\prime }}^{\ast }\right)=\frac{\partial Z\left({U^{\prime }}^{\ast }\right)}{\partial u_i}=\frac{\partial Z\left(X^{\ast }\right)}{\partial x_i}{{\sigma }_x}_i$$E19

$${\nabla }^{2}Z\left({U^{\prime }}^{\ast }\right)=\frac{{\partial }^{2}Z\left({U^{\prime }}^{\ast }\right)}{{\partial u}_i^2}=\frac{{\partial }^{2}Z\left(X^{\ast }\right)}{{\partial x}_i^2}{{\sigma }_x^2}_i$$E20

where ${{\sigma }_x}_i$ is the standard deviation of the random variable $x_i$.

To summarize, the main objective of using SORM is to improve the estimate of the probability of failure that was obtained by FORM. The calculation includes three key steps. These are, the calculation of the reliability index and its associated design points either in the U domain or X domain using FORM, determining the principal curvatures, and determining the probability of failure using one of the commonly used methods for nonlinear LSFs.

2.1 Probability of failure

Breitung, Hohenbichler, and Tvedt methods use a correction factor which is expressed in terms of $k_i$ and $\beta$ to adjust the probability of failure obtained by FORM. The probability of failure for nonlinear LSF is then formulated as:

$$P_f=\Phi \left(-\beta \right)\left(\mathit{CF}\right)$$E21

The first term is the cumulative distribution function (cdf) of beta, $\mathit{CF}$ is the SORM correction factor and $\beta$ is the reliability index obtained by FROM.

2.1.1 Breitung formulation

Breitung formulation is expressed as:

$${P_f}_{\mathit{\text{Breitung}}}=\Phi \left(-\beta \right)\prod _{i=1}^{n-1}\frac{1}{\sqrt{\left(1+\beta k_i\right)}}$$E22

2.1.2 Hohenbichler formulation

Hohenbichler formulation is expressed in terms of the cdf of beta (the first term), probability density function of beta (pdf) which is the upper term in the denominator and the cdf of beta which is the lower term in the denominator as:

$${P_f}_{\mathit{\text{Hohenbichler}}}=\Phi \left(-\beta \right)\prod _{i=1}^{n-1}\frac{1}{\sqrt{\left(1+\frac{\varphi \left(-\beta \right)}{\Phi \left(-\beta \right)}{k}_i\right)}}$$E23

2.1.3 Tvedt formulation

Tvedt formulation has three parts as indicated below:

$${P_f}_{\mathit{\text{Tvedt}}}=\Phi \left(-\beta \right)f1+f2+f3$$E24

$$f1=\prod _{i=1}^{n-1}\frac{1}{\sqrt{\left(1+\beta k_i\right)}}$$

$$f2=\left(\beta \mathrm{\Phi }\left(-\beta \right)-\varphi \left(-\beta \right)\right)\left(f1-\left(\prod _{i=1}^{n-1}\frac{1}{\sqrt{\left(1+\left(1+\beta \right)k_i\right)}}\right)\right)$$

$$f3=\left(1+\beta \right)\left(\beta \mathrm{\Phi }\left(-\beta \right)-\varphi \left(-\beta \right)\right)\left(f1-\mathit{\text{Real}}\left(\prod _{i=1}^{n-1}\frac{1}{\sqrt{\left(1+\left(j+\beta \right)k_i\right)}}\right)\right)$$

2.2 Calculation steps

• Calculate the reliability index $\beta$ using FORM.

• Calculate the design points $U^{\ast }$ associated with the reliability index $\beta$ following the steps in Appendix A.

• Determine the length of the gradient vector at the MPPF: $\left|\nabla Z\left({U^{\prime }}^{\ast }\right)\right|=\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(U^{\ast }\right)}{\partial u_i}\right)}^2}$ .

• Construct the matrix ${\mathbf{R}}_{\mathbf{01}}$ using Eq. (11).

• Perform Gram-Schmidt orthogonalization for the matrix ${\mathbf{R}}_{\mathbf{01}}$ using Eqs. (12) and (13) and perform orthonormalization of each row vector to come up with the matrix $\mathit{R}$.

• Compute the second-order derivative of the LSF at the design point, $U^{\ast }$ , using Eq. (15) to obtain the $\mathit{H}$ matrix.

• Compute the $\mathit{B}$ matrix, $\mathit{B}\mathit{=}\mathit{RH}{\mathit{R}}^{\mathit{t}}$.

• Compute the eigenvalues of the matrix $\mathit{B}$ to obtain the principal curvatures ($k_i$’s).

• Calculate the probability of failure $P_f$ using Eqs. (22)–(24).

3. Examples

3.1 Example 1

The limit state function for a system has been formulated as:

$$Z\left(X\right)=2.5x_1^2-x_2^2-0.167\text{where}{x}_1\sim N\left(2,0.2\right),x_2\sim N\left(2,0.32\right).$$

1. Determine the reliability index ($\beta$) and the Most Probable Point of Failure (MPPF) in the $X$ space and the $U$ space. Use a tolerance of 0.0001 for convergence.

2. Determine the probability of failure using SORM by applying Breitung, Hohenbichler, and Tvedt methods.

3. Compare the results of parts A and B with results obtained by Monte Carlo simulation. Monte Carlo results are shown in Table 1.

Number of simulation cycles	${\mathrm{P}}_{\mathrm{f}}$
2e5	0.00596500
1e6	0.00601600

Table 1.

MCS results – Example 1.0, Part C.

Solution

Part A

Determine $\beta$ and the MPPF in the $X$ space and the $U$ space. FORM is used to solve part A. The computation continues until the solution converges with a tolerance error of 0.0001 $\left(\frac{\left|{\beta }_i-{\beta }_{i-1}\right|}{{\beta }_{i-1}}<0.0001\right).$

Assume the initial MPPF as the given mean of each variable.

Transform the LSF $Z\left(\mathit{X}\right)$ to $Z\left(\mathit{U}\right)$:

$$Z\left(\mathit{U}\right)=2.5{\left(0.2u_1+2\right)}^2-{\left(0.32u_2+2\right)}^2-0.167$$

Compute the initial values of the design point in $U$ space:

$$x_1={\mu }_{x_1}=2$$

$$x_2={\mu }_{x_2}=2$$

$$u_1=\frac{x_1-{\mu }_{x_1}}{{{\sigma }_x}_1}=\frac{2-2}{0.2}=0$$

$$u_2=\frac{x_2-{\mu }_{x_2}}{{{\sigma }_x}_2}=\frac{2-2}{0.32}=0$$

Compute the initial estimate of $Z\left(U\right)$:

$$Z\left(0,0\right)=5.83300$$

Compute the partial derivatives of the LSF:

$$\frac{\partial Z\left(U^{\ast }\right)}{\partial u_1}=0.2u_1+2=2$$

$$\frac{\partial Z\left(U^{\ast }\right)}{\partial u_2}=-0.2048u_2-1.28=-1.28$$

Compute the standard deviation of the LSF:

$${\sigma }_z=\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(U^{\ast }\right)}{{\partial u}_i}\right)}^2}=\sqrt{{\left(2\right)}^2+{\left(-1.28\right)}^2}=2.3745$$

Compute the initial reliability index β:

$$\beta =\frac{{\mu }_z}{{\sigma }_z}=\frac{5.8330}{2.3745}=2.4565$$

Compute the directional cosines ${\alpha }_i$:

$${\alpha }_1=-\frac{\left(\frac{\partial Z\left(U^{\ast }\right)}{\partial u_1}\right)}{\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(U^{\ast }\right)}{\partial u_i}\right)}^2}}=-\frac{2}{2.3745}=-0.8423$$

$${\alpha }_2=-\frac{\left(\frac{\partial Z\left(U^{\ast }\right)}{\partial u_2}\right)}{\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(U^{\ast }\right)}{{\partial u}_i}\right)}^2}}=-\frac{-1.28}{2.3745}=0.5391$$

Iteration 1

Determine the new MPPF/design point in $U$ space and X space:

$$u_1=\beta {\alpha }_1=2.4565\ast \left(-0.8423\right)=-2.0690$$

$$u_2=\beta {\alpha }_2=2.4565\ast \left(0.5391\right)=1.3242$$

$$x_1=\beta {\alpha }_1{\sigma }_{x_1}+{\mu }_{x_1}=u_1{\sigma }_{x_1}+{\mu }_{x_1}=-2.0690\ast 0.2+2=1.5862$$

$$x_2=\beta {\alpha }_2{\sigma }_{x_2}+{\mu }_{x_2}=u_2{\sigma }_{x_2}+{\mu }_{x_2}=1.3242\ast 0.32+2=2.4237$$

Compute the LSF in terms of the new design points:

$$Z\left(U^{\ast }\right)=Z\left(-2.0690,1.3242\right)=0.2485$$

Compute the partial derivatives and the standard deviation of the LSF at the new design points:

The partial derivatives:

$$\frac{\partial Z\left(U^{\ast }\right)}{\partial u_1}=1.5862$$

$$\frac{\partial Z\left(U^{\ast }\right)}{\partial u_2}=-1.5512$$

The standard deviation of Z:

$${\sigma }_z=\sqrt{{\left(1.5862\right)}^2+{\left(-1.5512\right)}^2}=2.2186$$

Compute the new β in terms of the new design points:

$$\beta =\frac{Z\left({\mathit{U}}^{\ast }\right)-{\sum }_{i=1}^n\frac{\partial Z\left({\mathit{U}}^{\mathit{\ast }}\right)}{\partial u_i}\left(u_i^{\ast }\right)}{\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left({\mathit{U}}^{\ast }\right)}{{\partial u}_i}\right)}^2}}$$

$$=\frac{0.2485-\left(1.5862\ast \left(-2.0690\right)+\left(-1.5512\ast 1.3242\right)\right)}{2.2186}$$

$$\beta =2.5171$$

We examine the tolerance:

$$\frac{\left|{\beta }_i-{\beta }_{i-1}\right|}{{\beta }_{i-1}}=\frac{\left|2.5171-2.45650\right|}{2.45650}=0.02466$$

$$0.02466>0.0001$$

Since the error is greater than the established tolerance (0.0001), the computation should continue until convergence is obtained.

The next step is to compute the directional cosines:

$${\alpha }_1=-\frac{1.5862}{2.2186}=-0.7150$$

$${\alpha }_2=-\frac{-1.5512}{2.2186}=0.6992$$

The computation continues until β converges following the steps mentioned in Appendix A. Once convergence is obtained the probability of failure is calculated.

Table 2 shows the results for the next iterations. As the table indicates, convergence occurred at the third iteration with an error less than the tolerance.

Iteration	0	1	2	3
u1	0.0000	−2.0690	−1.7996	−1.7758
u2	0.0000	1.3242	1.7599	1.7762
x1	2.0000	1.5862	1.6401	1.6448
x2	2.0000	2.4237	2.5632	2.5684
σ1	0.2000	0.2000	0.2000	0.2000
σ2	0.3200	0.3200	0.3200	0.3200
Z(U^*)	5.8330	0.2485	−0.0122	0.0000
∂Z/∂u1	2.0000	1.5862	1.6401	1.6448
∂Z/∂u2	−1.2800	−1.5512	−1.6404	−1.6438
α1	−0.8423	−0.7150	−0.7070	−0.7073
α2	0.5391	0.6992	0.7072	0.7069
σ	2.3745	2.2186	2.3197	2.3254
β	2.45648	2.51711	2.51170	2.51171
Pf	0.007015	0.005916	0.006008	0.006007

Table 2.

Summary of results – Example 1.0, Part A.

$$\frac{\left|{\beta }_i-{\beta }_{i-1}\right|}{{\beta }_{i-1}}=\frac{\left|2.51171-2.51170\right|}{2.51170}=0.000003981$$

$$0.000003981\ll 0.0001$$

$\beta$ is calculated to be 2.51171 and the MPPF is located at $\mathit{U}$(-1.7758, 1.7762) and $\mathit{X}$(1.6448, 2.5684).

The probability of failure is computed in terms of the final value of $\beta :$

$$P_f=\mathrm{\Phi }\left(-\beta \right)=1-\mathrm{\Phi }\left(\beta \right)=1-\mathrm{\Phi }\left(2.51171\right)=0.006007$$

Part B

Determine the probability of failure using SORM by applying Breitung, Hohenbichler, and Tvedt methods.

Determine the length of the gradient vector at the MPPF:

The partial derivatives at the third iteration are:

$$\frac{\partial Z\left(U^{\ast }\right)}{\partial u_1}=1.6448$$

$$\frac{\partial Z\left(U^{\ast }\right)}{\partial u_2}=-1.6438$$

$$\left|\nabla Z\left(U^{\ast }\right)\right|=\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(U^{\ast }\right)}{\partial u_i}\right)}^2}=\sqrt{{\left(1.6448\right)}^2+{\left(-1.6438\right)}^2}=2.3254$$

Construct the matrix ${\mathbf{R}}_{\mathbf{01}}$ using Eq. (11):

$${\mathbf{R}}_{\mathbf{01}}=\left[\begin{array}{cc}1& 0\\ -\frac{\raisebox{1ex}{$\partial Z\left(U^{\ast }\right)$}\!\left/ \!\raisebox{-1ex}{${\partial u}_1$}\right.}{\left|\nabla Z\left(U^{\ast }\right)\right|}& -\frac{\raisebox{1ex}{$\partial Z\left(U^{\ast }\right)$}\!\left/ \!\raisebox{-1ex}{${\partial u}_2$}\right.}{\left|\nabla Z\left(U^{\ast }\right)\right|}\\ & \end{array}\right]=\left[\begin{array}{cc}1& 0\\ -0.7073& 0.7069\end{array}\right]$$

The second row is the directional cosines of the reliability index, $\beta$, at the MPPF.

Perform Gram–Schmidt orthogonalization for the matrix ${\mathit{R}}_{\mathbf{01}}$ using Eqs. (12) and (13), and perform orthonormalization of each row vector to come up with the matrix $\mathit{R}$:

$$r_n=r_{02}$$

The last element of the matrix is calculated as:

$$r_i=r_{0i}-\sum _{j=i+1}^n\frac{r_j{r}_{0i}^t}{r_j{r}_j^t}{r}_j$$

i = 2

$$r_2=r_{02}=\left[-0.70730.7069\right]$$

i = 1

$$r_1=r_{01}-\sum _{j=i+1}^n\frac{r_j{r}_{02}^t}{r_j{r}_j^t}{r}_j=\left[10\right]-\frac{r_1{r}_{02}^t}{r_1{r}_1^t}{r}_1$$

$$=\left[10\right]-\frac{\left[-0.70730.7069\right]\left[\begin{array}{c}1\\ 0\end{array}\right]}{\left[-0.70730.7069\right]\left[\begin{array}{c}-0.7073\\ 0.7069\end{array}\right]}\left[-0.70730.7069\right]$$

$$r_1=\left[0.49970.5000\right]$$

Normalizing the elements of the row vector, $r_2$ becomes:

$$r_2=\frac{\left[0.49970.5000\right]}{\sqrt{{0.4997}^2+{0.5000}^2}}=\left[0.70730.7069\right]$$

The matrix $\mathit{R}$ becomes:

$$\mathit{R}\mathit{=}\left[\begin{array}{cc}0.7069& 0.7073\\ -0.7073& 0.7069\end{array}\right]$$

Compute the second-order derivative of the LSF at the design point, $U^{\ast }$, using Equation (20) to obtain the $\mathit{H}$ matrix:

$$\mathit{H}\mathbf{=}\frac{1}{\left|\nabla Z\left(U^{\ast }\right)\right|}\left[\begin{array}{cc}\frac{{\partial }^{2}Z\left(U^{\ast }\right)}{{\partial u}_1^2}& \frac{{\partial }^{2}Z\left(U^{\ast }\right)}{\partial u_1\partial u_2}\\ \frac{{\partial }^{2}Z\left(U^{\ast }\right)}{\partial u_2\partial u_1}& \frac{{\partial }^{2}Z\left(U^{\ast }\right)}{{\partial u}_2^2}\end{array}\right]=\frac{1}{2.3254}\left[\begin{array}{cc}0.2& 0\\ 0& -0.2048\end{array}\right]=\left[\begin{array}{cc}0.0860& 0\\ 0& -0.0881\end{array}\right]$$

Compute $\mathit{RH}{\mathit{R}}^{\mathit{t}}$ matrix:

$\mathit{RH}{\mathit{R}}^{\mathit{t}}$ = $\left[\begin{array}{cc}-0.0011& -0.0870\\ -0.0870& -0.0010\end{array}\right]$

Compute the eigenvalues of the matrix $\mathit{B}$ to obtain principal curvatures ($k_i$’s):

The principal curvatures, ${k_i}^{\prime }s$, are computed by solving the eigenvalues of $\mathit{RH}{\mathit{R}}^{\mathit{t}}$. The last column and last row are dropped from the above matrix, then the eigenvalues are obtained to determine the $k_i$ values. From the above matrix, there is only one $k$at the MPPF, with a value of (−0.0011).

Calculate the probability of failure $P_f$ using Eqs. (22)–(24):

Breitung method

$${P_f}_{\mathit{\text{Breitung}}}=\mathrm{\Phi }\left(-\beta \right)\prod _{i=1}^{n-1}\frac{1}{\sqrt{\left(1+\beta k_i\right)}}$$

$$=\mathrm{\Phi }\left(-2.5117\right)\frac{1}{\sqrt{\left(1+2.5117(-0.0011\right))}}$$

$${P_f}_{\mathit{\text{Breitung}}}=0.00601565$$

Hohenbichler method

$${P_f}_{\mathit{\text{Hohenbichler}}}=\mathrm{\Phi }\left(-\beta \right)\prod _{i=1}^{n-1}\frac{1}{\sqrt{\left(1+\frac{\varphi \left(-\beta \right)}{\mathrm{\Phi }\left(-\beta \right)}{k}_i\right)}}$$

$$=\mathrm{\Phi }\left(-2.5117\right)\frac{1}{\sqrt{\left(1+\frac{\varphi \left(-2.5117\right)}{\mathrm{\Phi }\left(-2.5117\right)}\left(-0.0011\right)\right)}}$$

$${P_f}_{\mathit{\text{Hohenbichler}}}=0.00601671$$

Tvedt method

$$f1=\prod _{i=1}^{n-1}\frac{1}{\sqrt{\left(1+\beta k_i\right)}}=1.00137$$

$$f2=\left(\beta \mathrm{\Phi }\left(-\beta \right)-\varphi \left(-\beta \right)\right)\left(f1-\left(\prod _{i=1}^{n-1}\frac{1}{\sqrt{\left(1+\left(1+\beta \right)k_i\right)}}\right)\right)$$

$$=\left(-0.00193\right)\left(1.00137-\left(1.00192\right)\right)=1.05670E-06$$

$$f2=1.05670E-06$$

$$f3=\left(1+\beta \right)\left(\beta \mathrm{\Phi }\left(-\beta \right)-\varphi \left(-\beta \right)\right)\left(f1-\mathit{\text{Real}}\left(\prod _{i=1}^{n-1}\frac{1}{\sqrt{\left(1+\left(j+\beta \right)k_i\right)}}\right)\right)$$

$$=\left(1+2.5117\right)\left(2.5117\left(\mathrm{\Phi }\left(-2.5117\right)\right)-\varphi \left(-2.5117\right)\right)\left(1.00137-\mathit{\text{Real}}\left(\frac{1}{\sqrt{\left(1+(j+2.5117\right)\left(-0.0011\right))}}\right)\right)$$

$$f3=-3.03430E-9$$

$${P_f}_{\mathit{\text{Tvedt}}}=\mathrm{\Phi }\left(-\beta \right)f1+f2+f3$$

$$=\left(6.00743\mathrm{E}-3\right)\left(1.00137\right)+\left(1.05670E-06\right)+\left(-3.03430E-9\right)$$

$${P_f}_{\mathit{\text{Tvedt}}}=6.01671E-03$$

See Table 3.

Method	Reliability Index - β	${\mathrm{P}}_{\mathrm{f}}$
FORM	2.5117077	0.00600743
SORM – Breitung		0.00601565
SORM – Hohenbichler		0.00601671
SORM – Tvedt		0.00601671

Table 3.

Summary of SORM results – Example 1.0 Part B.

Part C

Compare the results of parts A and B with results obtained by Monte Carlo simulation:

The simulation was conducted using 2e5 and 1e6 simulation cycles as shown in Table 1. The results obtained by Breitung, Tvedt, and Hohenbichler methods are close to the Monte Carlo simulation results when using 1e6 simulation cycles.

3.2 Example 2

The performance function for a system has been formulated as:

$$Z\left(X\right)=24x_1^2+13x_2^2+x_3^2-100,x_1\sim N\left(2,0.2\right),x_2\sim N\left(2,0.32\right),x_3\sim N\left(4,0.4\right)$$

1. Determine the reliability index ($\beta$) and the most probable point of failure (MPPF) in the $X$ space and the$U$ space. Use a tolerance of 0.0001 for convergence.

2. Determine the probability of failure using SORM by applying Breitung, Hohenbichler, and Tvedt methods.

3. Compare the results of parts A and B with results obtained by Monte Carlo simulation. Monte Carlo results are shown in Table 4.

Number of simulation cycles	${\mathrm{P}}_{\mathrm{f}}$
2e5	0.0018250
1e6	0.0018360

Table 4.

MCS results – Example 2.0, Part C.

Solution

Part A

Determine $\beta$ and the MPPF in the $X$ space and the $U$ space. FORM is used to solve part A. The computation continues until the solution converges with a tolerance error of 0.0001. $\left(\frac{\left|{\beta }_i-{\beta }_{i-1}\right|}{{\beta }_{i-1}}<0.0001\right).$

Assume the initial MPPF as the given mean of each variable.

Transform the LSF $Z\left(\mathit{X}\right)$ to $Z\left(\mathit{U}\right)$:

$$Z\left(\mathit{U}\right)=24{\left(0.2u_1+2\right)}^2+13{\left(0.32u_2+2\right)}^2+{\left(0.4u_3+4\right)}^2-100$$

Compute the initial values of the design point in U space:

$$x_1={\mu }_{x_1}=2$$

$$x_2={\mu }_{x_2}=2$$

$$x_3={\mu }_{x_3}=4$$

$$u_1=\frac{x_1-{\mu }_{x_1}}{{{\sigma }_x}_1}=\frac{2-2}{0.2}=0$$

$$u_2=\frac{x_2-{\mu }_{x_2}}{{{\sigma }_x}_2}=\frac{2-2}{0.32}=0$$

$$u_3=\frac{x_3-{\mu }_{x_3}}{{{\sigma }_x}_3}=\frac{4-4}{0.4}=0$$

Compute the initial estimate of $Z\left(U\right)$:

$$Z\left(0,0,0\right)=64$$

Compute the partial derivatives of the LSF:

$$\frac{\partial Z\left(U^{\ast }\right)}{\partial u_1}=1.92u_1+19.2=19.2$$

$$\frac{\partial Z\left(U^{\ast }\right)}{\partial u_2}=2.6624u_2+16.64=16.64$$

$$\frac{\partial Z\left(U^{\ast }\right)}{\partial u_3}=0.32u_3+3.2=3.2$$

Compute the standard deviation of the LSF:

$${\sigma }_z=\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(U^{\ast }\right)}{{\partial u}_i}\right)}^2}=\sqrt{{\left(19.2\right)}^2+{\left(16.64\right)}^2+{\left(3.2\right)}^2}=25.6080$$

Compute the initial reliability index β:

$$\beta =\frac{{\mu }_z}{{\sigma }_z}=\frac{64}{25.6080}=2.4992$$

Compute the directional cosines ${\alpha }_i$:

$${\alpha }_1=-\frac{\left(\frac{\partial Z\left(U^{\ast }\right)}{\partial u_1}\right)}{\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(U^{\ast }\right)}{\partial u_i}\right)}^2}}=-\frac{19.2}{25.6080}=-0.7498$$

$${\alpha }_2=-\frac{\left(\frac{\partial Z\left(U^{\ast }\right)}{\partial u_2}\right)}{\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(U^{\ast }\right)}{{\partial u}_i}\right)}^2}}=-\frac{16.64}{25.6080}=-0.6498$$

$${\alpha }_3=-\frac{\left(\frac{\partial Z\left(U^{\ast }\right)}{\partial u_3}\right)}{\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(U^{\ast }\right)}{\partial u_i}\right)}^2}}=-\frac{3.2}{25.6080}=-0.1250$$

The initial calculation of the probability of failure is:

$$P_f=\mathrm{\Phi }\left(-\beta \right)=1-\mathrm{\Phi }\left(\beta \right)=1-\mathrm{\Phi }\left(2.4992\right)=0.00622$$

Iteration 1

Determine the new MPPF/design point in $U$ space and X space:

$$u_1=\beta {\alpha }_1=2.4992\ast \left(-0.7498\right)=-1.8738$$

$$u_2=\beta {\alpha }_2=2.4992\ast \left(-0.64980\right)=-1.6240$$

$$u_3=\beta {\alpha }_3=2.4992\ast \left(-0.1250\right)=-0.3123$$

$$x_1=\beta {\alpha }_1{\sigma }_{x_1}+{\mu }_{x_1}=u_1{\sigma }_{x_1}+{\mu }_{x_1}=-1.8738\ast 0.2+2=1.6252$$

$$x_2=\beta {\alpha }_2{\sigma }_{x_2}+{\mu }_{x_2}=u_2{\sigma }_{x_2}+{\mu }_{x_2}=-1.6240\ast 0.32+2=1.4803$$

$$x_3=\beta {\alpha }_3{\sigma }_{x_3}+{\mu }_{x_3}=u_3{\sigma }_{x_3}+{\mu }_{x_3}=-0.3123\ast 0.4+4=3.8751$$

Compute the LSF in terms of the new design points:

$$Z\left(U^{\ast }\right)=Z\left(-1.8738,-1.624,-0.3123\right)=6.8972$$

Compute the partial derivatives and the standard deviation of the LSF at the new design points:

The partial derivatives:

$$\frac{\partial Z\left(U^{\ast }\right)}{\partial u_1}=1.92u_1+19.2=15.6022$$

$$\frac{\partial Z\left(U^{\ast }\right)}{\partial u_2}=2.6624u_2+16.64=12.3163$$

$$\frac{\partial Z\left(U^{\ast }\right)}{\partial u_3}=0.32u_3+3.2=3.1001$$

The standard deviation of $Z$:

$${\sigma }_z=\sqrt{{\left(15.6022\right)}^2+{\left(12.3163\right)}^2+{\left(3.1001\right)}^2}=20.1179$$

Compute the new β in terms of the new design points:

$$\beta =\frac{Z\left({\mathit{U}}^{\ast }\right)-{\sum }_{i=1}^n\frac{\partial Z\left({\mathit{U}}^{\mathit{\ast }}\right)}{\partial u_i}\left(u_i^{\ast }\right)}{\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left({\mathit{U}}^{\ast }\right)}{{\partial u}_i}\right)}^2}}$$

$$=\frac{6.8972-\left(15.6022\ast \left(-1.8738\right)+12.3163\ast (-1.6240\right)+3.1001\ast \left(-0.3123\right))}{20.1179}$$

$$=2.8384$$

We examine the tolerance:

$$\frac{\left|{\beta }_i-{\beta }_{i-1}\right|}{{\beta }_{i-1}}=\frac{\left|2.8384-2.4992\right|}{2.4992}=0.1357$$

$$0.1357>0.0001$$

Since the error is greater than the established tolerance (0.0001), the computation should continue until convergence is obtained.

The next step is to compute the directional cosines:

$${\alpha }_1=-\frac{\left(78.011\ast 0.2\right)}{20.1179}=-0.7755$$

$${\alpha }_2=-\frac{\left(38.488\ast 0.32\right)}{20.1179}=-0.6122$$

$${\alpha }_3=-\frac{\left(7.7502\ast 0.4\right)}{20.1179}=-0.1541$$

The computation continues until β converges following the steps mentioned in Appendix A. Once convergence is obtained the probability of failure is calculated.

Table 5 shows the results for the next iterations. As the table indicates, the convergence occurred at the third iteration with an error less than the tolerance.

Iteration	0	1	2	3
u1	0.0000	−1.8738	−2.2013	−2.1911
u2	0.0000	−1.6240	−1.7377	−1.7580
u3	0.0000	−0.3123	−0.4374	−0.4478
x1	2.0000	1.6252	1.5597	1.5618
x2	2.0000	1.4803	1.4439	1.4375
x3	4.0000	3.8751	3.8250	3.8209
Z(U)	64.0000	6.8972	0.1227	0.0007
∂Z/∂u1	19.2000	15.6022	14.9735	14.9931
∂Z/∂u2	16.6400	12.3163	12.0136	11.9596
∂Z/∂u3	3.2000	3.1001	3.0600	3.0567
σ1	0.2000	0.2000	0.2000	0.2000
σ2	0.3200	0.3200	0.3200	0.3200
σ3	0.4000	0.4000	0.4000	0.4000
σZ	25.6080	20.1179	19.4396	19.4208
α1	−0.7498	−0.7755	−0.7703	−0.7720
α2	−0.6498	−0.6122	−0.6180	−0.6158
α3	−0.1250	−0.1541	−0.1574	−0.1574
β	2.49922	2.83840	2.84461	2.84463
Pf	0.006223	0.002267	0.002223	0.002223

Table 5.

Summary of the results – Example 2.0, Part A.

$$\frac{\left|{\beta }_i-{\beta }_{i-1}\right|}{{\beta }_{i-1}}=\frac{\left|2.84463-2.84461\right|}{2.84461}=0.000007031$$

$$0.000007031\ll 0.0001$$

$\beta$ is calculated to be 2.84463 and the MPPF is located at $U$(−2.1911, −1.7580, −0.4478) and$X$(1.5618, 1.4375, 3.8209).

The probability of failure is computed in terms of the final value of $\beta :$

$$P_f=\mathrm{\Phi }\left(-\beta \right)=1-\mathrm{\Phi }\left(\beta \right)=1-\mathrm{\Phi }\left(2.84463\right)=0.002223$$

Part B

Determine the probability of failure using SORM by applying Breitung, Hohenbichler, and Tvedt methods.

Determine the length of the gradient vector at the MPPF:

The partial derivatives at the third iteration are:

$$\frac{\partial Z\left(U^{\ast }\right)}{\partial u_1}=14.9931$$

$$\frac{\partial Z\left(U^{\ast }\right)}{\partial u_2}=11.9596$$

$$\frac{\partial Z\left(U^{\ast }\right)}{\partial u_3}=3.0567$$

$$\left|\nabla Z\left(U^{\ast }\right)\right|=\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(U^{\ast }\right)}{\partial u_i}\right)}^2}=\sqrt{{\left(14.9931\right)}^2+{\left(11.9569\right)}^2+{\left(3.0567\right)}^2}=19.4208$$

Construct the matrix ${\mathbf{R}}_{\mathbf{01}}$ using Eq. (11):

$${\mathbf{R}}_{\mathbf{01}}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ -\frac{\raisebox{1ex}{$\partial Z\left(U^{\ast }\right)$}\!\left/ \!\raisebox{-1ex}{${\partial u}_1$}\right.}{\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(U^{\ast }\right)}{\partial u_i}\right)}^2}}& -\frac{\raisebox{1ex}{$\partial Z\left(U^{\ast }\right)$}\!\left/ \!\raisebox{-1ex}{${\partial u}_2$}\right.}{\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(U^{\ast }\right)}{\partial u_i}\right)}^2}}& -\frac{\raisebox{1ex}{$\partial Z\left(U^{\ast }\right)$}\!\left/ \!\raisebox{-1ex}{${\partial u}_3$}\right.}{\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(U^{\ast }\right)}{\partial u_i}\right)}^2}}\end{array}\right]$$

$${\mathbf{R}}_{\mathbf{01}}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ -0.77201& -0.61581& -0.15739\end{array}\right]$$

The last row is the directional cosines of the reliability index, $\beta$, at the MPPF.

Perform Gram-Schmidt orthogonalization for the matrix ${\mathit{R}}_{\mathbf{01}}$ using Eqs. (12) and (13), and perform orthonormalization of each row vector to come up with the matrix $\mathit{R}$:

$$r_n=r_{0n}=r_{03}$$

$$r_i=r_{0i}-\sum _{j=i+1}^n\frac{r_j{r}_{0i}^t}{r_j{r}_j^t}{r}_j$$

i = 3

$$r_3=r_{03}=\left[\begin{array}{ccc}-0.77201& -0.61581& -0.15739\end{array}\right]$$

The remaining elements of the matrix are calculated as:

i = 2

$$r_2=r_{02}-\sum _{j=2}^n\frac{r_j{r}_{02}^t}{r_j{r}_j^t}{r}_j=\left[010\right]-\frac{r_3{r}_{01}^t}{r_3{r}_3^t}{r}_3$$

$$=\left[\begin{array}{ccc}0& 1& 0\end{array}\right]-\frac{\left[\begin{array}{ccc}-0.77201& -0.61581& -0.15739\end{array}\right]\left[\begin{array}{c}0\\ 1\\ 0\end{array}\right]}{\left[\begin{array}{ccc}-0.77201& -0.61581& -0.15739\end{array}\right]\left[\begin{array}{c}-0.77201\\ -0.61581\\ -0.15739\end{array}\right]}\left[\begin{array}{ccc}-0.77201& -0.61581& -0.15739\end{array}\right]$$

$$r_2=\left[\begin{array}{ccc}-0.47541& 0.62077& -0.096925\end{array}\right]$$

Normalizing the elements of the row vector, $r_2$ becomes:

$$r_2=\frac{\left[\begin{array}{ccc}-0.47541& 0.62077& -0.096925\end{array}\right]}{\sqrt{{\left(-0.47541\right)}^2+{0.62077}^2+{\left(-0.096925\right)}^2}}$$

$$r_2=\left[-0.60340.78789-0.12302\right]$$

Following the same procedures row vector $r_1$ becomes:

$$r_1=\left[0.199770-0.97984\right]$$

The matrix $\mathit{R}$ becomes:

$$\mathit{R}=\left[\begin{array}{ccc}0.19977& 0& -0.97984\\ -0.6034& 0.78789& -0.12302\\ -0.77201& -0.61581& -0.15739\end{array}\right]$$

Compute the second-order derivative of the LSF at the design point, $U^{\ast }$, using Equation (20) to obtain the $\mathit{H}$ matrix:

$$\mathit{H}\mathit{=}\frac{1}{\left|\nabla Z\left(U^{\ast }\right)\right|}\left[\begin{array}{ccc}\frac{{\partial }^{2}Z\left(U^{\ast }\right)}{{\partial u}_1^2}& \frac{{\partial }^{2}Z\left(U^{\ast }\right)}{\partial u_1\partial u_2}& \frac{{\partial }^{2}Z\left(U^{\ast }\right)}{\partial u_1\partial u_3}\\ \frac{{\partial }^{2}Z\left({\mathit{U}}^{\ast }\right)}{{\partial u}_2\partial u_1}& \frac{{\partial }^{2}Z\left(U^{\ast }\right)}{{\partial u}_2^2}& \frac{{\partial }^{2}Z\left(U^{\ast }\right)}{\partial u_2\partial u_3}\\ \frac{{\partial }^{2}Z\left(U^{\ast }\right)}{\partial u_3\partial u_1}& \frac{{\partial }^{2}Z\left(U^{\ast }\right)}{\partial u_3\partial u}& \frac{{\partial }^{2}Z\left(U^{\ast }\right)}{{\partial u}_3^2}\end{array}\right]$$

$$=\frac{1}{19.4208}\left[\begin{array}{ccc}1.9200& 0& 0\\ 0& 2.6624& 0\\ 0& 0& 0.3200\end{array}\right]=\left[\begin{array}{ccc}0.0989& 0& 0\\ 0& 0.1371& 0\\ 0& 0& 0.0165\end{array}\right]$$

Compute $\mathit{RH}{\mathit{R}}^{\mathit{t}}$ matrix:

$$\mathit{RH}{\mathit{R}}^{\mathit{t}}=\left[\begin{array}{ccc}0.0197650& -0.0099306& -0.0127056\\ -0.0099306& 0.1213464& -0.0201427\\ -0.0127056& -0.0201427& 0.1113186\end{array}\right]$$

Compute the eigenvalues of the matrix $\mathit{B}$ to obtain principal curvatures ($k_i$’s):

The principal curvatures, $k_i$’s, are computed by solving the eigenvalues of $\mathit{RH}{\mathit{R}}^{\mathit{t}}$. To do that, the last column and the last row are dropped from the matrix. The above matrix becomes:

$$\mathit{RH}{\mathit{R}}^{\mathit{t}}=\left[\begin{array}{cc}0.0197650& -0.0099306\\ -0.0099306& 0.1213464\end{array}\right]$$

Applying the eigenvalue method to obtain the principal curvatures:

$$k=\left[\begin{array}{cc}0.018803& 0.122308\end{array}\right]$$

Calculate the probability of failure $P_f$ using Eqs. (22)–(24):

Breitung method

$${P_f}_{\mathit{\text{Breitung}}}=\mathrm{\Phi }\left(-\beta \right)\prod _{i=1}^{n-1}\frac{1}{\sqrt{\left(1+\beta k_i\right)}}$$

$$=\mathrm{\Phi }\left(-2.8446\right)\frac{1}{\sqrt{\left(1+2.8446\left(0.018803\right)\right)\left(1+2.8446\left(0.122308\right)\right)}}$$

$${P_f}_{\mathit{\text{Breitung}}}=0.0018656$$

Hohenbichler method

$${P_f}_{\mathit{\text{Hohenbichler}}}=\mathrm{\Phi }\left(-\beta \right)\prod _{i=1}^{n-1}\frac{1}{\sqrt{\left(1+\frac{\varphi \left(-\beta \right)}{\mathrm{\Phi }\left(-\beta \right)}{k}_i\right)}}$$

$$=\frac{\mathrm{\Phi }\left(-2.8446\right)}{\sqrt{\left(1+\left(\frac{\varphi \left(-2.8446\right)}{\mathrm{\Phi }\left(-2.8446\right)}\right)\left(0.018803\right)\right)\left(1+\left(\frac{\varphi \left(-2.8446\right)}{\mathrm{\Phi }\left(-2.8446\right)}\right)\left(0.122308\right)\right)}}$$

$$=0.0018363$$

Tvedt method

$$f1=\prod _{i=1}^{n-1}\frac{1}{\sqrt{\left(1+\beta k_i\right)}}=0.83918$$

$$f2=\left(\beta \mathrm{\Phi }\left(-\beta \right)-\varphi \left(-\beta \right)\right)\left(f1-\left(\prod _{i=1}^{n-1}\frac{1}{\sqrt{\left((1+\left(1+\beta \right)k_i\right)}}\right)\right)$$

$$=\left(\left(2.8446\right)\left(0.0022231)-0.0069786\right)\right(0.83918-\left(0.9657)\left(0.82472\right)\right)$$

$$f2=-2.7977541E-05$$

$$f3=\left(1+\beta \right)\left(\beta \mathrm{\Phi }\left(-\beta \right)-\varphi \left(-\beta \right)\right)\left(f1-\mathit{\text{Real}}\left(\prod _{i=1}^{n-1}\frac{1}{\sqrt{\left(1+\left(j+\beta \right)k_i\right)}}\right)\right)$$

$$=\begin{array}{l}\left(1+2.8446\right)\left(2.8446\left(\mathrm{\Phi }\left(-2.8446\right)\right)-\varphi \left(-2.8446\right)\right)(0.83918\\ -\mathit{\text{Real}}\left(\frac{1}{\sqrt{\left(1+(j+2.8446\right)\left(0.018803\right)\left)\left(1+(j+2.8446\right)\left(0.122308\right)\right)}}\right))\end{array}$$

$$f3=-7.5839622E-06$$

$${P_f}_{\mathit{\text{Tvedt}}}=\mathrm{\Phi }\left(-\beta \right)f1+f2+f3$$

$=\left(0.0022231\right)\left(0.83918\right)+\left(-2.79775E-05\right)+\left(-7.58396E-06\right)$

$${P_f}_{\mathit{\text{Tvedt}}}\mathit{=}0.0018300$$

See Table 6.

Method	Reliability Index – $\mathrm{\beta }$	${\mathrm{P}}_{\mathrm{f}}$
FORM	2.844633	0.0022231
SORM–Breitung		0.0018656
SORM–Hohenbichler		0.0018363
SORM–Tvedt		0.0018300

Table 6.

Analytical methods results – Example 2.0, Part B.

Part C

Compare the results of parts A and B with results obtained by Monte Carlo simulation:

The simulation was conducted using 2e5 and 1e6 simulation cycles as shown in Table 4. The results obtained by Breitung, Tvedt, and Hohenbichler methods are close to the Monte Carlo simulation results when using 1e6 simulation cycles.

3.3 Example 3

Determine the second-order reliability index (${\beta }_{\mathit{\text{SORM}}}$) for the probability of failure calculated in example 2.

Solution

The reliability index for the second-order LSF can be calculated by taking the inverse of the cumulative distribution of the probability of failure as:

$${\beta }_{\mathit{\text{SORM}}}={\mathrm{\Phi }}^{-1}\left[1-P_f\right]$$

$${\beta }_{\mathit{\text{Breitung}}}={\mathrm{\Phi }}^{-1}\left[1-0.0018656\right]=2.90004$$

$${\beta }_{\mathit{\text{Tvedt}}}={\mathrm{\Phi }}^{-1}\left[1-0.0018300\right]=2.90606$$

$${\beta }_{\mathit{\text{Hohenbichler}}}={\mathrm{\Phi }}^{-1}\left[1-0.0018363\right]=2.90499$$

$${\beta }_{\mathit{MC}}={\mathrm{\Phi }}^{-1}\left[1-0.0018360\right]=2.90505$$

The absolute difference between Monte Carlo method and the other methods is calculated below:

$$\text{Breitung method}\mathbf{=}\left|2.90505-2.90004\right|=0.00500$$

$$\text{Tvedt method}\mathbf{=}\left|2.90505-2.90606\right|=0.00101$$

$$\text{Hohenbichler method}\mathbf{=}\left|2.90505-2.90499\right|=0.00006$$

In this example, the reliability index obtained by Hohenbichler is the closest to the reliability index obtained by Monte Carlo method.

4. Summary and concluding remarks

This chapter presented commonly used analytical probabilistic methods for determining the probability of failure for nonlinear and higher-order systems. They included Breitung, Hohenbichler, and Tvedt methods. Furthermore, Monte Carlo Simulation was used to estimate the probability of failure to validate the results obtained by the methods mentioned above. It is concluded that the results obtained by Breitung, Tvedt, and Hohenbichler methods are closer to the results obtained by Monte Carlo Simulation when using 1e6 simulation cycles. Future work should include the analysis of multi-dimensional and higher-order failure functions having correlated variables. Moreover, it is recommended that future work apply finite element method to gain more insight into how the methods compare with each other.

A. An overview of FORM

A.1 Failure modeling

The performance function, also called the limit state function (LSF), is formulated in terms of a system’s load and capacity. Both the load, $L$, and the capacity, $C$, have an impact on the performance of the system. The system fails when the load exceeds the capacity or when the LSF becomes less than zero. The LSF can be expressed as:

$$Z=C-L$$E25

The capacity $C$ and the load $L$ are formulated in terms of random variables,$x_1,x_{2,}\dots x_n$:

$$Z\left(\mathit{X}\right)=Z\left(x_1,x_2,\dots \dots \dots \dots ,x_n\right)$$E26

These variables are assumed to be statistically independent random variables having a normal distribution. The probability of failure can be expressed as:

$$P_f=P\left(Z<0\right)$$E27

This is the probability that the LSF becomes less than zero. In Figure A1, this is the region that exists above the straight line for linear LSF and above the curve for nonlinear LSF. The safe region or nonfailure region is below the curve, and this is the region where the LSF becomes greater than zero. The LSF becomes at the limit state when the LSF becomes equal to zero as Figure A1 indicates.

Figure A1.

Representation of linear and nonlinear LSF – original coordinates (X domain).

Also, the figure shows the reliability index, $\beta$, which is the shortest distance from the origin to the surface. It is calculated as the ratio of the mean, ${\mu }_Z$, to the standard deviation, ${\sigma }_Z$, of the LSF.

$$\beta =\frac{{\mu }_Z}{{\sigma }_Z}$$E28

$${\mu }_Z\approx Z\left({\mu }_{X_1,}{\mu }_{X_2},\dots \dots \dots \dots ..,{\mu }_{X_n}\right)$$E29

$${\sigma }_z=\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(X^{\ast }\right)}{{\partial x}_i}\right)}^2{\sigma }_{x_i}^2}$$E30

Then, the probability of failure is expressed in terms of $\beta$ as:

$$P_f=\varphi \left(-\beta \right)=1-\Phi \left(\beta \right)$$E31

A.2 First Order Reliability Method (FORM)

This section will use the Hasofer and Lind method to calculate the probability of failure by FORM [9]. The first step is to transform the normal random variables into standard normal variables as:

$$u_i=\frac{x_i-{\mu }_{x_i}}{{\sigma }_{x_i}},i=1,2,.\dots ,n$$E32

$$x_i=u_i\ast {\sigma }_{x_i}+{\mu }_{x_i},i=1,2,.\dots ,n$$E33

where ${\mu }_{x_i}$ and ${\sigma }_{x_i}$ are the mean and the standard deviation of the random variable $x_i$ respectively, and $u_i$ is the transformed standard normal variable. The LSF is then formulated in terms of the standard normal variables as:

$$Z\left(U\right)=Z\left(u_1\ast {\sigma }_{x_1}+{\mu }_{x_1},u_2\ast {\sigma }_{x_2}+{\mu }_{x_2},\dots ,u_n\ast {\sigma }_{x_n}+{\mu }_{x_n}\right)=0$$E34

Figure A2 shows the representation of linear and nonlinear LSF and the most probable point of failure, $(\mathit{\text{MPPF}}:u_1^{\ast }$, $u_2^{\ast })$, or sometimes referred to as the design point. The next subsection outlines the steps for calculating the probability of failure using FORM.

Figure A2.

Representation of linear and nonlinear LSF – transformed coordinates (U domain).

A.3 FORM calculation steps

The calculation steps for this method are listed below [10].

1. Formulate the LSF/performance function in terms of the original random variables, $x_i$:

   $$Z\left(\mathit{X}\right)=Z\left(x_1,x_2,\dots \dots \dots \dots .x_n\right)$$

2. Assume the initial MPPF/design point as the given mean of each variable.

3. Transform the LSF $Z\left(\mathit{X}\right)$ to $Z\left(\mathit{U}\right)$ using Eq. (34).

4. Compute the initial values of the design point in $U$ space using Eq. (32). It should be noted that the initial values will be zero.

5. Compute the initial estimate of $Z\left(U\right)$ in terms of the initial values of design point $U$.

6. Compute the partial derivatives of the LSF:

   $$\frac{\partial Z\left({\mathit{U}}^{\ast }\right)}{\partial u_i}$$E35

7. Compute the standard deviation of the LSF using the following Equation:

   $${\sigma }_z=\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(U^{\ast }\right)}{{\partial u}_i}\right)}^2}$$E36

8. Compute the initial reliability index β:

   $$\beta =\frac{{\mu }_Z}{{\sigma }_Z}$$E37

9. Compute the directional cosines ${\alpha }_i$ using the following Equation:

   $${\alpha }_i=-\frac{\left(\frac{\partial Z\left({\mathit{U}}^{\ast }\right)}{\partial u_i}\right)}{\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(U^{\ast }\right)}{\partial u_i}\right)}^2}}$$E38

10. Determine the new MPPF/design point in $U$ space using the following Equation and in the $X$ space using Eq. (33):

    $$u_i^{\ast }=\beta {\alpha }_i$$E39

11. Compute the LSF in terms of new design points.

12. Compute the partial derivatives, Eq. (35), and the standard deviation of the LSF, Eq. (36), at the new design points.

13. Compute the new β using the following Equation:

    $$\beta =\frac{Z\left({\mathit{U}}^{\ast }\right)-{\sum }_{i=1}^n\frac{\partial Z\left({\mathit{U}}^{\mathbf{\ast }}\right)}{\partial u_i}\left(u_i^{\ast }\right)}{\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left({\mathit{U}}^{\ast }\right)}{{\partial u}_i}\right)}^2}}$$E40

14. Repeat steps 9 through 13 until β converges. A tolerance, ε, of 0.001 is usually used:

    $$\epsilon =\frac{\left|{\beta }_i-{\beta }_{i-1}\right|}{{\beta }_{i-1}}<0.001$$E41

15. Use Eq. (31) to calculate the probability of failure.