Probabilistic Modeling of Failure

Abstract

Failure of a system or a component of a system is and has been a major concern to systems’ operators and owners. Failure could be traced back to different causes and may take different forms and shapes. It may result from software malfunction, hardware degraded performance, human error, sabotage, environmental as well as other external factors. There are various techniques found in the literature that can assist in the analysis of failure. These techniques comprise deterministic and probabilistic techniques. Deterministic techniques ignore the variability and uncertainties of the variables in the analysis which may lead to unsatisfactory and inaccurate results. While probabilistic techniques produce accurate and an all-inclusive result because they incorporate the variabilities and uncertainties in the analysis. The focus of this chapter is to present commonly used probabilistic failure analysis techniques and their mathematical derivations. Examples to enhance the understanding of the concept of failure analysis are also presented.

1. Introduction

Traditionally, failure analysis is conducted using deterministic techniques to assess the operability and integrity of industrial systems. These techniques lack the ability to report or predict the probabilistic nature of the systems’ behavior. Moreover, they ignore the probabilistic and random nature of the external factors that have direct impact on the performance of the systems. Implementing these techniques may produce inadequate assessment and eventually, may lead to wrong decisions concerning the integrity and reliability of the evaluated systems. To make an informed and reliable decision about the reliability and operability of such systems, probabilistic failure analysis should be adopted as an alternative analysis technique. This technique should be made as an integral part of the decision process as well as be a part of the overall organization’s risk control.

This chapter presents techniques that assist in the analysis of failure using engineering probabilistic methods. They include simulation as well as analytical methods. Simulation methods can be conducted using Monte Carlo simulation technique. Two different Monte Carlo simulation approaches are presented in this chapter. These are, the counting approach and sample statistics approach. The main drawback with simulation is that it takes great deal of time to perform and may require an extensive processing power. However, it is an essential step in the analysis to validate the results obtained by the analytical methods.

Some of the analytical methods include first order reliability method (FORM) and second order reliability method (SORM). FORM involves two approaches to calculate probability of failure, these are first order and second moment (FOSM) and advanced FOSM.

The focus of this chapter will be only on the FORM with the assumption that all random variables are uncorrelated. Analysis that require the use of FORM for correlated random variables is beyond the scope of this chapter. Likewise, analysis requiring the application of SORM to analyze limit state functions involving second order representation is beyond the scope of this chapter.

2. Failure modeling

Failure can be defined as the inability of industrial systems or subsystems either partially or totally to satisfy operational requirements as set forth by design specifications. Failure of a system could be partial or complete; in either case the consequences of failure may result in adverse consequences. Interruption of services, degraded performance, system shutdown, environmental damage and customers’ dissatisfaction are some of the consequences. Such consequences may lead to financial losses, liabilities and destroyed image of the operating company. As an example, if failure involves leak detection system to detect oil and gas leakage from subsea oil and gas pipelines, consequences could be so severe. Pollution of the ocean, damage to the fishery and tourism industries are some of the major consequences.

The system fails when the imposed demand or load on the system exceeds its capacity or resistance. The strength or the capacity of the system is a design parameter that specifies the maximum load the system can endure or the maximum demand the system can satisfy. The variabilities of the system’s capacity to satisfy the demand or load imposed on it are mainly attributed to the inherent uncertainties of the operation characteristics of the system’s components as well as external environmental factors. Therefore, the capacity of the system is assumed to be probabilistic in nature that varies from time to time due to the reasons mentioned above. Likewise, using same argument the load or demand imposed on the system are considered probabilistic in nature due to the effect of the varying environmental conditions.

Considering the above, the performance function of the system or sometime is called limit state function can be formulated as the difference between the system’s capacity and the load or demand imposed on it. The same argument can be used for production facility, the performance is the difference between supply and demand, supply being the capacity, or the strength and demand is the load. If the two parameters are the same then it can be said that the system is at a limit state, if the system cannot meet the demand then the system is at a failure state, and if the system capacity exceeds the load imposed on it, the system is at a satisfactory state.

3. Reliability analysis methods

Knowing in advance when the system is going to fail or degrade in performance is an essential step in the failure analysis. Under this step, the probability of failure is calculated in terms of the random variables affecting the performance of the system. There are several approaches found in the literature that can be used to evaluate the probability of failure either analytically or by simulation. Analytical methods approximate the probability of failure by using first order reliability method (FORM) or second order reliability method (SORM). The FORM uses two approximation techniques that evaluates the probability of failure, these are the first order and second moment (FOSM) and advanced first order second moment (AFOSM) techniques.

Calculating the probability of failure based on the methods mentioned above can be used to predict the ability of the system to satisfy operational as well as safety requirements during its life cycle. Combining this analysis with risk analysis, the consequences of failure can be easily determined. First order reliability method consists of two techniques namely:

3.1 Analytical methods

3.1.1 First order second moment (FOSM)

FOSM makes use of second moment statistics (mean and variance) and ignores higher moments (skewness and kurtosis) of the random variables. It evaluates the performance function by using the first order Taylor series expansion of the limit state function (LSF) at the mean value. This method is used when the performance function is linear having statistically independent, normally distributed and noncorrelated random variables ${\mathrm{X}}_{\mathrm{i}}^{\prime }\mathrm{s}$.

Performance function can be defined as [3, 5, 10]:

$$\mathrm{Z}=\mathrm{C}-\mathrm{D}$$E1

where C is the capacity and D is the demand are statistically independent random variables and are assumed to be normally distributed. Failure occurs when:

$$Z<0$$E2

Then the probability of failure ($P_f$) can be computed as:

$$P_f=(P\left(C<D,)\right)$$E3

Or

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

Figure 1 shows the probability density function (PDF) of the random variable Z, as it can be depicted from the figure that the probability of failure comprises the shaded area where Z < 0.

Figure 1.

Probability density function of the performance function.

The probability of failure is expressed as [3, 5, 10]:

$$P_f=\left(P\left(Z<0\right)\right)=\underset{-\infty }{\overset{0}{\int }}{f}_{z\left(Z\right)\mathit{dz}}$$E5

Alternatively, the performance function $Z$ can be formulated in terms of many random variables designated by a vector X as [3, 5, 10]:

$$Z\left(\mathit{X}\right)=Z\left(X_1,X_{2,}\dots .X_n\right)$$E6

$X_1,\dots .X_n$ are the random variables in the performance function.

The integration of the performance function as indicated in Eq. (5) is performed for the region where $\mathrm{Z}<0$, this type of integration is difficult to solve, alternatively Taylor series expansion is used. The first order Taylor series approximation about the mean of the random variables is shown in Eq. (7). The expansion is truncated at the linear terms to obtain the first order approximation of the performance function [7, 8, 9].

$$Z=Z\left({\mu }_X\right)+\sum _{i=1}^n\frac{\partial Z}{\partial X_i}\left(X_i-{\mu }_{X_i}\right)+\frac{1}{2!}\sum _{i=1}^n \sum _{j=1}^n\frac{{\partial }^{2}Z}{\partial X_i\partial X_j}\left(X_i-{\mu }_{X_i}\right)\left(X_j-{\mu }_{X_j}\right)$$E7

$$Z=Z\left({\mu }_X\right)+\sum _{i=1}^n\frac{\partial Z}{\partial X_i}\left(X_i-{\mu }_{X_i}\right)$$E7.1

Then the mean and variance are given as:

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

$$\mathit{Var}\left[Z\right]=\mathit{Var}\left[Z\left({\mu }_X\right)+\sum _{i=1}^n\frac{\partial Z}{\partial X_i}\left(X_i-{\mu }_{X_i}\right)\right]$$E9

$$\mathit{Var}\left[Z\right]=\mathit{Var}\left[Z\left({\mu }_X\right)\right]+\mathit{Var}\left[\sum _{i=1}^n\frac{\partial Z}{\partial X_i}\left(X_i\right)\right]-\mathit{Var}{\left[\sum _{i=1}^n\frac{\partial Z}{\partial X_i}\left({\mu }_{X_i}\right)\right]}_{X_i}$$E10

$$\mathit{Var}\left[Z\left({\mu }_X\right)\right]=0$$E10.1

$$\mathit{Var}{\left[\sum _{i=1}^n\frac{\partial Z}{\partial X_i}\left({\mu }_{X_i}\right)\right]}_{X_i}=0$$E10.2

$$\mathit{Var}\left[Z\right]={\left[\sum _{i=1}^n\frac{\partial Z}{\partial X_i}\right]}^2\mathit{Var}\left(X_i\right)$$E10.3

$${\sigma }_Z=\sqrt{{\left[\sum _{i=1}^n\frac{\partial Z}{\partial X_i}\right]}^2\mathit{Var}\left(X_i\right)}$$E10.4

The reliability index (β) is taken as ratio of the mean to the standard deviation of the performance function.

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

The reliability index is computed for every failure mode, where the probability of failure is expressed as:

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

This method is simple to use and assumes that the random variables are normally distributed. All is needed for the calculation is the knowledge of the mean and the standard deviation and it is not necessary to know the distribution of the random variables. The downside of this method is that it can cause error in the final results if the function is nonlinear or if the tail of the distribution cannot be approximated by normal distribution. Moreover, if the function is nonlinear it will be provided different answer than that if it is linear. Advanced FOSM is used to deal with the limitations of the FOSM mentioned above.

3.1.2 Advanced first order second moment method (AFOSM)

AFOSM provides solution for linear and nonlinear performance function by determining the shortest distance from the origin to the failure surface. This method is also called Hasofer-Lind method. It evaluates the probability of failure for the limit state function or the performance function by determining the most probable failure point instead of the mean. Hasofer and Lind developed this advanced method in 1974 which is called Hasofer-Lind method and is abbreviated as H-L method. As stated above, the main objective of this method is to estimate the failure point which is the shortest distance from the origin to the failure surface that separates the failure region from the safe region. This can be clearly shown in Figure 2. The failure point is sometimes called in the literature design point or check point, but in this chapter, it will be referred to as the most probable point of failure (MPPF). Let us consider a limit state function/performance function with normally distributed and independent random variables X as:

Figure 2.

(a) Nonlinear limit state function—original coordinates. (b) Nonlinear limit state function—transformed coordinates.

$$Z\left(X\right)=Z\left(x_1,x_2,\cdots \cdots \cdots \cdots \cdots \cdots ,x_n\right)=0$$E13

This method transforms the random variables into reduced form as:

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

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

The performance function is then formulated in terms of the reduced random 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 \dots \dots \dots .u_n\ast {\sigma }_{x_n}+{\mu }_{x_n}\right)=0$$E16

Figure 2 shows the plot of the limit state function in the original as well as the transformed coordinates. It shows that the MPPF is the tangent point on the curve $Z\left(X\right)=0$ and the reliability index $\beta$ as the shortest distance from the origin to the limit surface.

To find the MPPF ${\mathrm{x}}_{\mathrm{i}}^{\prime }$ on the limit surface under the condition that $\mathrm{Z}\left(\mathbf{X}\right)=0$, Taylor series expansion is used around the MPPF, considering the first order terms only, this gives:

$$Z\left(U\right)=Z\left(U^{\ast }\right)+\sum _{i=1}^n\frac{\partial Z\left(U^{\ast }\right)}{{\partial U}_i}\left(u_i-u_i^{\ast }\right)$$E17

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

Using chain rule for derivative and considering that relationship between U and X as:

$$Z\mathit{\text{function of}}x\to Z=g\left(x\right)$$

$$X\mathit{\text{function of}}u\to X=f\left(u\right)$$

and using Eqs. (14) and (15) the partial derivative $\frac{\partial g\left(U^{\ast }\right)}{\partial u_i}$ becomes:

$$\frac{\partial Z}{{\partial U}_i}=\frac{\partial Z}{\partial X_i}\frac{\partial X}{{\partial U}_i}=\left(\frac{\partial Z}{\partial X_i}\right)\left(\frac{\partial X\left(u_i^{\ast }\ast {\sigma }_{x_i}+{\mu }_{x_i}\right)}{{\partial U}_i}\right)=\frac{\partial Z}{\partial X_i}{\sigma }_{x_i}$$E19

Substituting Eq. (19) into Eq. (17) gives:

$$Z\left(\mathit{U}\right)=Z\left({\mathit{U}}^{\ast }\right)+\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}{\sigma }_{x_i}\left(u_i-u_i^{\ast }\right)$$E20

The mean of $Z\left(U\right)$ is:

$$u_z=E\left[Z\left({\mathit{U}}^{\ast }\right)\right]+E\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}{\sigma }_{x_i}{u}_i\right]-E\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}{\sigma }_{x_i}{u}_i^{\ast }\right]$$E21

$$E\left[Z\left(U^{\ast }\right)\right]=Z\left(U^{\ast }\right)$$E22

$$E\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}{\sigma }_{x_i}{u}_i\right]=E\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}\left(x_i-u_{x_i}\right)\right]=E\left[x_i\right]\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}\right]-E\left[{\mu }_{x_i}\right]\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}\right]$$E23

$${\mu }_{x_i}\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}\right]-{\mu }_{x_i}\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}\right]=0$$E24

$$E\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}{\sigma }_{x_i}{u}_i^{\ast }\right]=\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}{\sigma }_{x_i}{u}_i^{\ast }$$E25

$$u_z=Z\left({\mathit{U}}^{\ast }\right)-\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\mathbf{\ast }}\right)}{{\partial x}_i}{\sigma }_{x_i}{u}_i^{\ast }$$E26

The variance is expressed as:

$${\sigma }_Z^2=\mathit{Var}(Z\left(\mathit{U}\right)=\mathit{Var}\left[Z\left({\mathit{U}}^{\ast }\right)\right]+\mathit{Var}\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}{\sigma }_{x_i}\left(u_i-u_i^{\ast }\right)\right]$$E27

It must be noted that constants have no variance, their variances equal to zero. The first term of the Taylor expansion in Eq. (27) is constant; therefore, its variance equals to zero. Similarly, the variance at the mean value is zero.

$$\mathit{Var}\left[Z\left({\mathit{U}}^{\ast }\right)\right]=0$$E28

$$\mathit{Var}(Z\left(\mathit{U}\right)=\mathit{Var}\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}{\sigma }_{x_i}\left(u_i\right)\right]\mathbf{-}\mathit{Var}\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}{\sigma }_{x_i}\left(u_i^{\ast }\right)\right]$$E29

$$\mathit{Var}\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}{\sigma }_{x_i}\left(u_i^{\ast }\right)\right]=0$$E30

$$\mathit{Var}(Z\left(\mathit{U}\right)=\mathit{Var}\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}{\sigma }_{x_i}\left(u_i\right)\right]\mathbf{=}\mathit{Var}\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}\left(x_i-u_{x_i}\right)\right]$$E31

$$\mathit{Var}\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}\left(x_i-u_{x_i}\right)\right]\mathit{=}\mathit{Var}\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}{x}_i\right]-\mathit{Var}\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}{u}_{x_i}\right]$$E31.1

$$\mathit{Var}\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}{u}_{x_i}\right]\mathbf{=}\mathbf{0}$$E32

$$\mathit{Var}(Z\left(\mathit{U}\right)\mathbf{=}\sum _{i=1}^n{\left(\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}\right)}^2\mathit{Var}\left[x_i\right]$$E33

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

Alternatively, Eq. (31) can be calculated as [6]:

$$\begin{array}{cc}\mathit{Var}\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}\left(x_i-u_{x_i}\right)\right]& =\mathit{E}{\left[\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}\left(x_i-u_{x_i}\right)\right]}^{\mathbf{2}}\\ & =\sum _{i=1}^n{\left(\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}\right)}^{2}E{\left[\left(x_i-u_{x_i}\right)\right]}^2\\ & =\sum _{i=1}^n{\left(\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}\right)}^2\mathit{Var}\left[x_i\right]\end{array}$$E34.1

The reliability index is calculated as:

$$\beta =\frac{{\mu }_Z}{{\sigma }_Z}=\frac{Z\left({\mathit{X}}^{\ast }\right)-\sum _{i=1}^n\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{\partial x_i}{\sigma }_{x_i}\left(u_i^{\ast }\right)}{\sqrt{\sum _{i=1}^n{\left(\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{{\partial x}_i}{\sigma }_{x_i}\right)}^2}}$$E35

The directional cosine ${\alpha }_i$ along the coordinate axes is computed as:

$${\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}}=-\frac{\left(\frac{\partial Z\left({\mathit{X}}^{\ast }\right)}{\partial x_i}{\sigma }_{x_i}\right)}{\sqrt{\sum _{i=1}^n{\left(\frac{\partial Z\left(\mathit{X}\right)}{\partial x_i}{\sigma }_{x_i}\right)}^2}}$$E36

It can be shown from Figure 2 that:

$$u_i^{\ast }=\beta \ast \mathit{Cos}\left({\theta }_i\right)\mathit{\text{or}}{u}_i^{\ast }={\alpha }_{x_i}\beta$$E37

Using Eqs. (14), (15) and (36) we can determine the design point in the original coordinates as:

$$x_i^{\ast }=\beta {\alpha }_{x_i}{\sigma }_{x_i}+{\mu }_{x_i}=u_i^{\ast }{\sigma }_{x_i}+{\mu }_{x_i}i=1,2,\dots \dots \dots \dots \dots ,n$$E38

The probability of failure, $P_f$ can be computed as:

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

The steps to estimate the reliability index are:

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

   $$Z\left(x_1,x_2,\cdots \cdots \cdots \cdots \cdots \cdots ,x_n\right)$$

2. Assume the initial design points as the given mean of each variables.

3. Compute the initial reliability index β in terms of the mean values of the random variables using Eq. (11).

4. Compute the partial derivatives of the limit state function (LSF)/performance function in terms of mean value of the random variables.

5. Compute the directional cosines ${\alpha }_i$using Eq. (36).

6. Compute the new design point using Eq. (38).

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

8. Compute the partial derivatives at the new design points.

9. Compute the new reliability index β using Eq. (35).

   1. Another alternative is to compute β from the limit state function $\left(Z\left(X\right)=Z\left(\beta {\alpha }_{x_i}{\sigma }_{x_i}+{\mu }_{x_i}\right)=0\right)$using the newly determined design points in step 7 and solve for β.

10. Compute the directional cosine.

11. Repeat steps six through nine until β converges to a pre-established tolerance level.

12. Use Eq. (39)to calculate the probability of failure.

The steps mentioned above are used with the assumption that the random variables are normally distributed. For non-normally distributed variables additional steps are needed to determine the mean and standard deviation of the equivalent normal distribution as listed below. These steps should be carried out after step number two to determine mean and standard deviation of the equivalent normal distribution. Assuming the random variables are statistically independent and non-normally distributed:

1. Determine the distribution parameters.

2. Compute the cumulative distribution function cdf, $F\left(x_i\right)$, the probability density function pdf, $f\left(x_i\right)$ and the inverse cdf, ${\mathrm{\Phi }}^{-1}\left[F\left(x_i\right)\right]$ of the original non-normal random variables at the initial design point.

3. Compute the values of the standard deviation, ${\sigma }_{x_i}$ and the mean, ${\mu }_{x_i}$ of the equivalent normal distribution as:

   $${\sigma }_{x_i}=\frac{\varphi \left({\mathrm{\Phi }}^{-1}\left[F\left(x_i\right)\right]\right)}{f\left(x_i\right)}$$E40

   It must be noted that $f\left(x_i\right)$ refers to the pdf of the original non-normal random variable and $\varphi \left(\right)$ refers to the pdf of the equivalent standard normal random variable.

   $${\mu }_{x_i}=x_i-{\mathrm{\Phi }}^{-1}\left[F\left(x_i\right)\right]{\sigma }_{x_i}$$E41

4. Compute the standard normal variable ($u_i)$ of $x_i$ by:

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

For a log normally distributed random variable, distribution parameters ${\mu }_{\mathit{ln}\left(x\right)}$, ${\sigma }_{\mathit{ln}\left(x\right)}$ are defined using the following equations:

$${\sigma }_{ln\left(x\right)}^2=\mathit{ln}\left({\left(\frac{{\sigma }_x}{{\mu }_x}\right)}^2+1\right),{\sigma }_{ln\left(x\right)}=\sqrt{{\sigma }_{ln\left(x\right)}^2}$$E43

$${\mu }_{\mathit{ln}\left(x\right)}=ln\left({\mu }_x\right)-0.5{\sigma }_{ln\left(x\right)}^2$$E44

The pdf and cdf are defined as:

$$f\left(x\right)=\frac{1}{\sqrt{2\pi }\left(x\right){\sigma }_{ln\left(x\right)}}exp\left(-0.5{\left(\frac{ln\left(x\right)-{\mu }_{\mathit{ln}\left(x\right)}}{{\sigma }_{\mathit{ln}\left(x\right)}}\right)}^2\right)$$E45

$$F\left(x\right)=\mathrm{\Phi }\left(\frac{ln\left(x\right)-{\mu }_{ln\left(x\right)}}{{\sigma }_{ln\left(x\right)}}\right)$$E46

For other distribution types the readers are referred to Refs. [4, 7, 8].

3.2 Simulation methods

Alternatively, the probability of failure is computed using Monte Carlo simulation method. Two methods are considered in this chapter, the counting and sample statistics methods. The simulation is conducted using computer programs such MATLAB, C++ or MINITAB or any other simulation programming packages.

3.2.1 Monte Carlo counting method

The counting method is formulated by dividing the number of simulation cycles at the events when the Z function becomes less than 0 (${\mathrm{N}}_{\mathrm{f}}<0)$ by the total number of simulation cycles (N).

$$P_f=\frac{N_f}{N}$$E47

The steps for Monte Carlo simulation counting method are listed below:

Counting method

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

   $$Z\left(x_1,x_2,\cdots \cdots \cdots \cdots \cdots \cdots ,x_n\right)$$
   .

2. Determine the distribution and its parameters for each random variable.

3. Assign N to be the number of simulation cycles.

4. Assign M to be the number of calculation times.

5. Assign Nf to be the number of simulation cycles when the Z function becomes less than 0.

6. Initialize:

   1. Set M to zero

   2. Set Nf to zero

   7. Generate random values from the given distribution. With the determined distribution parameters for each variable.

   a. As an example, if the random variable is normally distributed use the following Matlab function to generate the random values:

   $$x=\mathit{\text{normrnd}}\left(\mathit{mu},\mathit{var},n,1\right)$$E48

   N is number of simulation cycles, mu and var. are the distribution parameters (the mean and standard deviation of the random variable).

   b. If the random variable is log normally distributed with a mean, ${\mu }_{x_i}$and standard deviation, ${\sigma }_{x_i}$:

   i. Determine the distribution parameters; Eqs. (43)and (44)can be used use to calculate these two parameters.

   ii. Use the following Matlab function to generate the random values:

   $$x=\mathit{\text{lognrnd}}\left({\mu }_{\mathit{ln}\left(x\right)},{\sigma }_{\mathit{ln}\left(x\right)},N,1\right)$$E49

7. Calculate Pf:

   1. Is Z < 0?

      1. If yes:

         • Nf = Nf +1,

         • M = M + 1,

      2. If no:

         • Is M = N?

           1. If no—go to step 7

           2. If yes—calculate $\mathit{Pf}=\frac{\mathit{Nf}}{N}$

Alternative approach to random number generation is to use the following steps:

• Generate uniformly distributed random numbers $u_i$between 0 and 1; this can be accomplished by using software packages such as MATLAB, excel or C++ and other software packages.

• Equate the inverse cumulative distribution function cdf of the random variable to the generated random numbers, $u_i$by using the following equation:

$$u_i=F_x\left(x_i\right)$$E49.1

$$x_i=F_x^{-1}\left(u_i\right)$$E49.2

As an example, if $x_i$follows normal distribution with a mean ${\mu }_{x_i}$, and standard deviation ${\sigma }_{x_i}$then its random number becomes:

$$x_i={\mu }_{x_i}+{\sigma }_{x_i}\left(F_x^{-1}\left(u_i\right)\right)$$E49.3

$\left(F_x^{-1}\left(u_i\right)\right)$can be determined using the cdf of the standard normal distribution, these tables are included in most probability and statistics books.

The random numbers can be generated N times using the following MATALB function:

$$u_i=\mathit{rand}\left(N,1\right)$$E49.4

3.2.2 Monte Carlo sample statistics method

The Monte Carlo sample statistics methods considers the mean ${\mu }_z$and the standard deviation ${\sigma }_z$in computing the reliability index.

$$\beta =\frac{{\mu }_z}{{\sigma }_z}$$E50

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

The steps for Monte Carlo simulation sample statistics method are listed below:

Sampling method

1. Formulate the performance function in terms of the original random variables, $x_i$: $Z\left(x_1,x_2,··\cdots \cdots \cdots \cdots \cdots \cdots ,x_n\right)$.

2. Determine the distribution and its parameters for each random variable.

3. Assign N to be the number of simulation times.

4. Assign M to be the number of calculation times.

5. Initialize:

   1. Set M to zero

   2. Set Nf to zero

6. Generate random values from the given distribution with the determined distribution parameters for each variable.

7. Calculate Z function.

8. Calculate Pf:

   c. Calculate $\beta =$ $\frac{{\mu }_z}{{\sigma }_z}$

   d. Calculate $P_f=1-\Phi \left(\beta \right)$

   e. Is M = N?

   i.If no—go to step 6

   ii.If yes—stop

It must be noted that Monte Carlo sample statistics method can be used only for linear functions having uncorrelated normal random variables.

Example 1

The performance function for a system has been formulated as: $z=2.5x_1-x_2$

where

$$x_1\sim N\left(2,0.2\right),x_2\sim N\left(3.1,0.32\right)$$

1. Determine the probability of failure using simulation methods: Monte Carlo (MCS) sample statistics and Monte Carlo counting simulation methods.

2. Determine the probability of failure using the analytical methods: FOSM and AFOSM (Hasofer-Lind) methods.

Solution:

1. Monte Carlo simulation was conducted for both methods, the sample statistics and counting methods. The number of simulation cycles used in the analysis are 2e5 and 1e6 cycles (Table 1).

2. See Table 2.

	MCS—sample statistics method		MCS—counting method
Number of simulation cycles	Reliability index—β	Pf	Reliability index—β	Pf
2e5	3.1993	0.000689	3.1967	0.000695
1e6	3.2030	0.00068	3.2005	0.000686

Table 1.

Simulation methods results—Example 1.

FOSM		AFOSM
Reliability index—β	Pf	Reliability index—β	Pf
3.2006	0.000686	3.2006	0.000686

Table 2.

Analytical methods results—Example 1.

^**

This example adopted from [9].

Example 2

The performance function for a leak detection system has been formulated as [2]:

$z=2.8252x^{-1}-x_2$, $x_1\sim N\left(1.1,0.02\right),x_2\sim N\left(2.37,0.04\right)$

1. Determine the probability of failure using Monte Carlo counting simulation method.

2. Determine the probability of failure using AFOSM method.

This example is adopted from reference [1, 2], with some modifications.

Solution

Part a

Monte Carlo simulation was conducted for both methods, the sample statistics and counting methods as indicated in the table.

Table 3 indicates that the probability of failure converges to 0.000476.

	MCS—counting method
Number of simulation cycles	Reliability index—β	Pf
1e5	3.2962	0.000490
1e6	3.3043	0.000476

Table 3.

Simulation methods results—Example 2.

Part b

First the Z function is computed in terms of the mean values of the random variables:

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

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

$$u_1=\frac{x_1-{\mu }_{x_1}}{{{\sigma }_x}_1}=\frac{1.1-1.1}{0.02}=0$$

$$u_2=\frac{x_2-{\mu }_{x_2}}{{{\sigma }_x}_2}=\frac{2.37-2.37}{0.04}=0$$

$$Z\left(X\right)=Z\left(\mathit{1.1},\mathit{2.37}\right)=0.1984$$

The partial derivatives:

$$\frac{\partial z}{\partial x_1}=\frac{-2.8252}{{x_1}^2}=-2.3349$$

$$\frac{\partial z}{\partial x_2}=-1$$

$${\sigma }_z=\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(X\right)}{\partial x_i}{\sigma }_{x_i}\right)}^2}=\sqrt{{\left(\mathit{-}\mathbf{2.3349}\mathit{\ast }\mathbf{0.02}\right)}^2+{\left(\mathit{-}\mathbf{1}\mathit{\ast }\mathbf{0.04}\right)}^2}=0.061487$$

$${\alpha }_1=-\frac{\left(\frac{\partial Z\left(X\right)}{\partial x_1}{\sigma }_{x_1}\right)}{\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(X\right)}{\partial x_i}{\sigma }_{\mathit{xi}}\right)}^2}}=-\frac{\left(-2.3349\ast 0.02\right)}{0.061487}=0.75947$$

$${\alpha }_2=-\frac{\left(\frac{\partial Z\left(X\right)}{\partial x_2}{\sigma }_{x_2}\right)}{\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial g\left(X\right)}{\partial x_i}{\sigma }_{\mathit{xi}}\right)}^2}}=-\frac{\left(-1\ast 0.04\right)}{0.061487}=0.65054$$

$$\beta =\frac{{\mu }_z}{{\sigma }_z}=\frac{0.19836}{0.061487}=3.2261$$

$$P_f=\mathrm{\Phi }\left(-\beta \right)=1-\mathrm{\Phi }\left(\beta \right)=1-\mathrm{\Phi }\left(3.2261\right)=0.000627$$

Iteration 1

$$u_1=\beta {\alpha }_{x_1}=3.2261\ast 0.75947=2.45013$$

$$u_2=\beta {\alpha }_{x_2}=3.2261\ast 0.65054=2.09872$$

$$x_1=\beta {\alpha }_1{\sigma }_{x_1}+{\mu }_{x_1}=u_1{\sigma }_{x_1}+{\mu }_{x_1}=2.45013\ast 0.02+1.1=1.149003$$

$$x_2=\beta {\alpha }_2{\sigma }_{x_2}+{\mu }_{x_2}=u_2{\sigma }_{x_2}+{\mu }_{x_2}=2.09872\ast 0.04+2.37=2.453949$$

$$Z\left(X\right)=Z\left(\mathit{1.149003},\mathit{2.453949}\right)=0.004880$$

$$\frac{\partial Z}{\partial x_1}=\frac{-2.8252}{{x_1}^2}=\frac{-2.8252}{{\left(1.149003\right)}^2}=-2.13997$$

$$\frac{\partial Z}{\partial x_2}=-1$$

$${\alpha }_1=-\frac{\left(\frac{\partial Z\left(X\right)}{\partial x_1}{\sigma }_{x_1}\right)}{\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(X\right)}{\partial x_i}{\sigma }_{\mathit{xi}}\right)}^2}}=-\frac{\left(-2.13997\ast 0.02\right)}{\sqrt{{\left(-2.13997\ast 0.02\right)}^2+{\left(-1\ast 0.04\right)}^2}}=\frac{0.042799}{0.05858}=0.7306$$

$${\alpha }_2=-\frac{\left(\frac{\partial Z\left(X\right)}{\partial x_2}{\sigma }_{x_2}\right)}{\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(X\right)}{\partial x_i}{\sigma }_{\mathit{xi}}\right)}^2}}=-\frac{\left(-1\ast 0.04\right)}{\sqrt{{\left(-2.13997\ast 0.02\right)}^2+{\left(-1\ast 0.04\right)}^2}}=0.6828$$

$$\begin{array}{cc}{\beta }_{\mathit{HL}}& =\frac{Z\left(X\right)-{\sum }_{i=1}^n\frac{\partial Z\left(X\right)}{\partial x_i}{\sigma }_{x_i}\left(u_i\right)}{\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z}{\partial x_i}{\sigma }_{x_i}\right)}^2}}\\ & =\frac{0.004880-(\left(-2.13997\ast 0.02\ast 2.45013\right)+\left(-1\ast 0.04\ast 2.09872\right)}{0.05858}\\ & =\frac{0.19369}{0.05858}=3.3064\end{array}$$

$$P_f=\mathrm{\Phi }\left(-\beta \right)=1-\mathrm{\Phi }\left(\beta \right)=1-\mathrm{\Phi }\left(3.3064\right)=0.0004726$$

The probability of failure obtained by counting method is very close to that obtained by AFOSM (Table 4).

Iteration	x1	x2	u1	u2	Z(X)	∂z/∂x1	∂z/∂x2	α1	α2	β	pf
1	1.100000	2.370000	0.000000	0.000000	0.198364	−2.334876	−1	0.759469	0.650543	3.226104	0.000627
2	1.149003	2.453949	2.450126	2.098721	0.004880	−2.139968	−1	0.730596	0.682810	3.306375	0.000473
3	1.148312	2.460305	2.415624	2.257626	0.000001	−2.142540	−1	0.731005	0.682372	3.306389	0.000473
4	1.148340	2.460248	2.416987	2.256188	0.000000	−2.142439	−1	0.730989	0.682389	3.306389	0.000473

Table 4.

Summary of the analytical methods results—Example 2.

Example 3

A pipeline segment is suffering corrosion that grows annually at steady rate. The extent of the initial growth has been estimated to be 4.7 mm and it is assumed that it follows log normal distribution with standard deviation of 1.1.

The corrosion annual growth follows log normal distribution with a mean and standard deviation values of 0.2 and 0.01. The pipeline wall thickness follows normal distribution with a mean and standard deviation values of 14 mm and 4.7 respectively. The critical pipeline wall thickness has been determined to be 80% of the wall thickness [3]. A summary of the relevant information pertaining to the pipeline corrosion is presented in the net table, Table 5.

	Wall thickness	Initial corrosion	Corrosion growth
Distribution type	Normal distribution	Log normal distribution	Log normal distribution
Mean	14 mm	4.7 mm	0.2 mm
Standard deviation	0.07	1.1	0.01

Table 5.

Pipeline corrosion data—Example 3.

The owner of the pipeline decides not to repair the corrosion and wants to know if the pipeline can survive for the next 14 years without causing a leakage. It has been decided that in order to be in the safe side the maximum acceptable probability of failure has been set to 1e4 [3].

Solve this problem using analytical method as well as Monte Carlo (MCS) simulation method.

Solution

Monte Carlo simulation

$x_1$: Wall thickness, $x_1\sim N\left(1.1,0.02\right)$

$x_2$: Extent of the initial corrosion, $x_2\sim \mathit{LGN}\left(4.7,1.1\right)$

$x_3$: Corrosion annual growth,$x_3\sim \mathit{LGN}\left(0.2,0.01\right)$

The limit state function can be formulated as indicated below:

The capacity portion is formulated as: $0.8x_1$

The demand portion is formulated as: $x_2+14x_3$

$$Z\left(X\right)=0.8x_1-\left(x_2+14x_3\right)$$

Pipeline wall thickness:

$${\mu }_{x_1}=14,{\sigma }_{x_1}=0.07$$

$$x_1=\mathit{\text{normrnd}}\left({\mu }_{x_1},{\sigma }_{x_1},N,1\right)$$

Initial corrosion:

$${\mu }_{x_2}=4.7,{\sigma }_{x_2}=1.1$$

$${\sigma }_{ln\left(x_2\right)}^2=\mathit{ln}\left({\left(\frac{{\sigma }_{x_2}}{{\mu }_{x_2}}\right)}^2+1\right)=0.0533$$

$${\sigma }_{ln\left(x_2\right)}=\sqrt{{\sigma }_{ln\left(x_2\right)}^2}=0.2309$$

$${\mu }_{\mathit{ln}\left(x_2\right)}=ln\left({\mu }_{x_2}\right)-0.5{\sigma }_{ln\left(x_2\right)}^2=1.5209$$

$$x_2=\mathit{\text{lognrnd}}\left({\mu }_{ln\left(x_2\right)},{\sigma }_{ln\left(x_2\right)},N,1\right)$$

Corrosion growth rate:

$${\mathrm{\mu }}_{{\mathrm{x}}_3}=0.2,{\mathrm{\sigma }}_{{\mathrm{x}}_3}=0.01$$

$${\sigma }_{ln\left(x_3\right)}^2=\mathit{ln}\left({\left(\frac{{\sigma }_{x_3}}{{\mu }_{x_3}}\right)}^2+1\right)=0.002497$$

$${\sigma }_{ln\left(x_3\right)}=\sqrt{{\sigma }_{ln\left(x_3\right)}^2}=0.04997$$

$${\mu }_{\mathit{ln}\left(3\right)}=ln\left({\mu }_{x_3}\right)-0.5{\sigma }_{ln\left(x_3\right)}^2=-1.61069$$

$$x_3=\mathit{\text{lognrnd}}\left({\mu }_{ln\left(x_3\right)},{\sigma }_{ln\left(x_3\right)},N,1\right)$$

Monte Carlo simulation produced the following results as shown in Table 6. The probability of failure and beta converge to the following values:

	MCS—counting method
Number of simulation cycles	Reliability index—β	Pf
1e6	2.6249	0.004334
1e7	2.6206	0.004389

Table 6.

Simulation methods results—Example 3.

$$\mathit{Pf}=\frac{\mathit{Nf}}{N}=\frac{4396}{1e6}=0.004389$$

$$\mathit{\text{Beta}}={\Phi }^{-1}\left(\mathit{Pf}\right)=2.6206$$

Here only the counting method is used because the sampling method produces different results. The sampling method produces accurate results for linear and normal limit state/performance function only.

Analytical Solution

Assume the initial value for each random variable to be its mean.

$$Z\left(X\right)=Z\left(\mathrm{14,4.7,0.2}\right)=3.7$$

$$\frac{\partial Z}{\partial x_1}=0.8$$

$$\frac{\partial Z}{\partial x_2}=-1$$

$$\frac{\partial Z}{\partial x_3}=-14$$

For non-normal variables, the standard deviation and mean values of the equivalent normal variables are calculated using Eqs. (40) and (41).

Pipeline wall thickness, ${\mathit{x}}_{\mathbf{1}}$:

$${\mu }_{x_1}=14,{\sigma }_{x_1}=0.07$$

Initial corrosion, ${\mathit{x}}_{\mathbf{2}}$:

$${\mu }_{x_2}=4.7,{\sigma }_{x_2}=1.1$$

$${\sigma }_{ln\left(x_2\right)}^2=\mathit{ln}\left({\left(\frac{{\sigma }_{x_2}}{{\mu }_{x_2}}\right)}^2+1\right)=0.0533$$

$${\sigma }_{ln\left(x_2\right)}=\sqrt{{\sigma }_{ln\left(x_2\right)}^2}=0.2309$$

$${\mu }_{\mathit{ln}\left(x_2\right)}=ln\left({\mu }_{x_2}\right)-0.5{\sigma }_{ln\left(x_2\right)}^2=1.5209$$

Compute the pdf of the original non-normal variable (log normal distribution) using Eq. (45):

$$\begin{array}{cc}f\left(x\right)& =\frac{1}{\sqrt{2\pi }\left(x\right){\sigma }_{ln\left(x\right)}}exp\left(-0.5{\left(\frac{ln\left(x\right)-{\mu }_{ln\left(x\right)}}{{\sigma }_{ln\left(x\right)}}\right)}^2\right)\\ & =\frac{1}{\sqrt{2\pi }\left(4.7\right)\left(0.2309\right)}exp\left(-0.5{\left(\frac{ln\left(4.7\right)-1.5209}{0.2309}\right)}^2\right)=0.36512\end{array}$$

Compute the cdf of the original non-normal variable (log normal distribution) using Eq. (46):

$$F\left(x\right)=\mathrm{\Phi }\left(\frac{ln\left(x_2\right)-{\mu }_{\mathit{ln}\left(x\right)}}{{\sigma }_{\mathit{ln}\left(x\right)}}\right)=\mathrm{\Phi }\left(\frac{ln\left(4.7\right)-1.5209}{0.2309}\right)$$

$${\mathrm{\Phi }}^{-1}\left[F\left(x_i\right)\right]=\left(\frac{ln\left(4.7\right)-1.5209}{0.2309}\right)=0.115466$$

Compute the mean and standard deviation of the equivalent normal variable at the design point using Eqs. (40) and (41):

$${\sigma }_{x_2}=\frac{\varphi \left({\mathrm{\Phi }}^{-1}\left[F\left(x_2\right)\right]\right)}{f\left(x_2\right)}=\frac{\varphi \left(0.115466\right)}{0.36512},$$

$\varphi \left(\right)$ refers to the equivalent standard normal random variable and is calculated as:

$$\begin{array}{cc}\varphi \left(0.115466\right)& =\frac{1}{\sqrt{2\pi }}exp\left(-0.5{\left(\frac{ln\left(x_2\right)-{\mu }_{ln\left(x\right)}}{{\sigma }_{ln\left(x\right)}}\right)}^2\right)\\ & =\frac{1}{\sqrt{2\pi }}exp\left(-0.5{\left(0.115466\right)}^2\right)=0.3963\end{array}$$

$${\sigma }_{x_2}=\frac{0.3963}{0.36512}=1.0854$$

$${\mu }_{x_2}=x_2-{\mathrm{\Phi }}^{-1}\left[F\left(x_2\right)\right]{\sigma }_{x_2}=4.7-\mathit{0.115466}\left(\mathit{1.0854}\right)=\mathit{4.5747}$$

Corrosion growth, ${\mathit{x}}_{\mathbf{3}}$

The same procedures outline for initial corrosion extent, ${\mathrm{x}}_2$: to produces the following results:

$$f\left(x_3\right)=39.907,{\mathrm{\Phi }}^{-1}\left[F\left(x_3\right)\right]=0.02526,\varphi \left({\mathrm{\Phi }}^{-1}\left[F\left(x_3\right)\right]\right)=0.398815,$$

$${\sigma }_{x_3}=0.0099936\mathit{\text{and}}{\mu }_{x_3}=0.19975$$

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

$$\beta =\frac{{\mu }_z}{{\sigma }_z}=\frac{3.7}{1.0958}=3.3766$$

$$P_f=\mathrm{\Phi }\left(-\beta \right)=0.00037$$

$${\mathrm{\alpha }}_1=-\frac{\left(\frac{\mathrm{\partial Z}\left(\mathrm{X}\right)}{\partial {\mathrm{x}}_1}{\mathrm{\sigma }}_{{\mathrm{x}}_1}\right)}{\sqrt{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left(\frac{\mathrm{\partial Z}\left(\mathrm{X}\right)}{\partial {\mathrm{x}}_1}{\mathrm{\sigma }}_{\mathrm{x}1}\right)}^2}}=-0.0511$$

$${\alpha }_2=-\frac{\left(\frac{\partial Z\left(X\right)}{\partial x_2}{\sigma }_{x_2}\right)}{\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(X\right)}{\partial x_2}{\sigma }_{x2}\right)}^2}}=0.9905$$

$${\alpha }_3=-\frac{\left(\frac{\partial Z\left(X\right)}{\partial x_3}{\sigma }_{x_3}\right)}{\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z\left(X\right)}{\partial x_3}{\sigma }_{x3}\right)}^2}}=0.1277$$

Iteration 1

$$u_1=\beta {\alpha }_{x_1}=-0.1013$$

$$u_2=\beta {\alpha }_{x_2}=3.3445$$

$$u_3=\beta {\alpha }_{x_3}=0.4311$$

$$x_1=\beta {\alpha }_1{\sigma }_{x_1}+{\mu }_{x_1}=u_1{\sigma }_{x_1}+{\mu }_{x_1}=13.9879$$

$$x_2=\beta {\alpha }_2{\sigma }_{x_2}+{\mu }_{x_2}=u_2{\sigma }_{x_2}+{\mu }_{x_2}=8.2047$$

$$x_3=\beta {\alpha }_3{\sigma }_{x_3}+{\mu }_{x_3}=u_2{\sigma }_{x_3}+{\mu }_{x_3}=0.20406$$

$$\frac{\partial g}{\partial x_1}=0.8,\frac{\partial g}{\partial x_2}=-1,\frac{\partial g}{\partial x_3}=-14$$

Mean and standard deviation for normal variable ${\mathrm{x}}_1$ remain the same:

$${\mu }_{x_1}=14,{\sigma }_{x_1}=0.07$$

Mean and standard deviation for the equivalent normal distribution for the non-normal variables, ${\mathrm{x}}_2\text{and}{\mathrm{x}}_3$ (Table 7):

Iteration	1	2	3	4	5	6
x1	14.0000	13.9879	13.9929	13.9953	13.9947	13.9947
x2	4.7000	8.2047	9.9091	8.5641	8.4595	8.3734
x3	0.2000	0.2041	0.2023	0.2015	0.2017	0.2016
u1	0.0000	−0.1726	−0.1013	−0.0675	−0.0753	−0.0750
u2	0.0000	3.3445	3.4276	2.7576	2.6609	2.6164
u3	0.0000	0.4311	0.2582	0.1706	0.1896	0.1889
mu1	14.0000	14.0000	14.0000	14.0000	14.0000	14.0000
mu2	4.5747	3.4147	2.2538	3.1972	3.2620	3.3145
mu3	0.1998	0.1997	0.1997	0.1997	0.1997	0.1997
σ1	0.0700	0.0700	0.0700	0.0700	0.0700	0.0700
σ2	1.0854	1.8947	2.2883	1.9777	1.9536	1.9337
σ3	0.0100	0.0102	0.0101	0.0101	0.0101	0.0101
Z(x)	3.7000	0.1288	−1.5475	−0.1882	−0.0869	−0.0006
∂z/∂x1	0.8000	0.8000	0.8000	0.8000	0.8000	0.8000
∂z/∂x2	−1.0000	−1.0000	−1.0000	−1.0000	−1.0000	−1.0000
∂z/∂x3	−14.0000	−14.0000	−14.0000	−14.0000	−14.0000	−14.0000
σ	1.0958	1.9009	2.2933	1.9835	1.9594	1.9396
α1	−0.0511	−0.0295	−0.0244	−0.0282	−0.0286	−0.0289
α2	0.9905	0.9967	0.9978	0.9971	0.9970	0.9969
α3	0.1277	0.0751	0.0617	0.0711	0.0720	0.0727
β	3.3766	3.4388	2.7637	2.6687	2.6243	2.6240
Pf	0.0004	0.0003	0.0029	0.0038	0.0043	0.0043

Table 7.

Analytical methods results—Example 3.

$${\mu }_{x_2}=3.4147,{\sigma }_{x_2}=1.8947,{\mu }_{x_3}=0.1997,{\sigma }_{x_3}=0.010197$$

$${\alpha }_1=-0.0295,{\alpha }_2=0.9967,{\alpha }_3=0.0751$$

${\sigma }_z=\sqrt{{\left(0.8\ast 0.07\right)}^2+{\left(-1\ast 1.8947\right)}^2+{\left(-14\ast 0.010197\right)}^2}$ $=1.9009$

$$Z\left(X\right)=Z\left(13.9879,8.2047,0.20406\right)=0.1288$$

$${\beta }_{\mathit{HL}}=\frac{Z\left(X\right)-{\sum }_{i=1}^n\frac{\partial Z\left(X\right)}{\partial x_i}{\sigma }_{x_i}\left(u_i\right)}{\sqrt{{\sum }_{i=1}^n{\left(\frac{\partial Z}{\partial x_i}{\sigma }_{x_i}\right)}^2}}=\frac{0.1288-\left(-6.4080\right)}{1.9009}=3.4388$$

$$P_f=\mathrm{\Phi }\left(-\beta \right)=0.00037$$

Table 8 shows that the obtained values for beta and the probability of failure are so close to each other.

Analytical method		MCS—counting method
Reliability index—β	Pf	Reliability index—β	Pf
2.6240	0.004345	2.6206	0.004389

Table 8.

Comparison of results—Example 3.

The calculated probability of failure exceeds the target probability of failure ${10}^{-4}$.

Example 4

Two systems, system A and system B, each system has three main components and each component with it; the probability distribution type and its parameters are shown in Table 9. It is required to determine the probability of failure of each component using FOSM, AFOSM and Monte Carlo simulation [9].

System A	Distribution type and its parameters
Electric motor	Capacity	Demand
System energy	N∼(0.52, 0.05)	N∼(0.46, 0.08)
Number of parts	N∼(2, 0.5)	N∼(1, 0.3)
Weight	N∼(0.42 0.05)	N∼(0.400, 0.06)
System B	Distribution type and its parameters
Pneumatic cylinder	Capacity	Demand
System energy	N∼(0.52, 0.05)	N∼(0.4, 0.04)
Number of parts	N∼(2, 0.5)	N∼(1, 0.3)
Weight	N∼(0.420, 0.05)	N∼(0.37, 0.04)

Table 9.

Capacity and demand variables—Example 4.

Solution

Formulate a general LSF for each component as:

$$Z\left(X\right)=x_1-x_2$$

Formulate a general LSF for each component see the (Table 10).

System A	MCS—sample statistics method		MCS—counting method
Number of simulation cycles	Reliability index—β	Pf	Reliability index—β	Pf
Electric motor—system energy
1e5	0.6334	0.2632	0.6326	0.2635
1e6	0.6363	0.2623	0.6363	0.2623
Electric motor—number of parts
1e5	1.7207	0.0427	1.7197	0.0427
1e6	1.7168	0.0430	1.7183	0.0429
Electric motor—weight
1e5	0.2556	0.3991	0.2572	0.3985
1e6	0.2558	0.3991	0.2564	0.3988
System B	MCS—sample statistics method		MCS—counting method
Number of simulation cycles	Reliability index—β	Pf	Reliability index—β	Pf
Pneumatic cylinder—system energy
1e5	1.8808	0.0300	1.8927	0.0292
1e6	1.8748	0.0304	1.8770	0.0303
Pneumatic cylinder—number of parts
1e5	1.7088	0.0437	1.7036	0.0442
1e6	1.7150	0.0432	1.7156	0.0431
Pneumatic cylinder—weight
1e5	0.7808	0.2175	0.7801	0.2177
1e6	0.781443467	0.2172709	0.78127631	0.21732

Table 10.

Simulation methods results systems A and B—Example 4.

Analytical Solution

Analytical solution see (Table 11).

Electric motor—system energy
Iteration	x1	x2	u1	u2	σ1	σ2	Z(X)	∂z/∂x1	∂z/∂x2	α1	α2	σ_z	β	pf
1	0.5200	0.4600	0.0000	0.0000	0.0500	0.0800	0.0600	1.0000	−1.0000	−0.5300	0.8480	0.0943	0.6360	0.2624
2	0.5031	0.5031	−0.3371	0.5393	0.0500	0.0800	0.0000	1.0000	−1.0000	−0.5300	0.8480	0.0943	0.6360	0.2624
Electric motor—number of parts
Iteration	x1	x2	u1	u2	σ1	σ2	Z(X)	∂z/∂x1	∂z/∂x2	α1	α2	σ_z	β	pf
1	2.0000	1.0000	0.0000	0.0000	0.5000	0.3000	1.0000	1.0000	−1.0000	−0.8575	0.5145	0.5831	1.7150	0.0432
2	1.2647	1.2647	−1.4706	0.8824	0.5000	0.3000	0.0000	1.0000	−1.0000	−0.8575	0.5145	0.5831	1.7150	0.0432
Electric motor—weight
Iteration	x1	x2	u1	u2	σ1	σ2	Z(X)	∂z/∂x1	∂z/∂x2	α1	α2	σ_z	β	pf
1	0.4200	0.4000	0.0000	0.0000	0.0500	0.0600	0.0200	1.0000	−1.0000	−0.6402	0.7682	0.0781	0.2561	0.3989
2	0.4118	0.4118	−0.1639	0.1967	0.0500	0.0600	0.0000	1.0000	−1.0000	−0.6402	0.7682	0.0781	0.2561	0.3989
Pneumatic cylinder—system energy
Iteration	x1	x2	u1	u2	σ1	σ2	Z(X)	∂z/∂x1	∂z/∂x2	α1	α2	σ_z	β	pf
1	0.5200	0.4000	0.0000	0.0000	0.0500	0.0400	0.1200	1.0000	−1.0000	−0.7809	0.6247	0.0640	1.8741	0.0305
2	0.4468	0.4468	−1.4634	1.1707	0.0500	0.0400	0.0000	1.0000	−1.0000	−0.7809	0.6247	0.0640	1.8741	0.0305
Pneumatic cylinder—number of parts
Iteration	x1	x2	u1	u2	σ1	σ2	Z(X)	∂z/∂x1	∂z/∂x2	α1	α2	σ_z	β	pf
1	2.0000	1.0000	0.0000	0.0000	0.5000	0.3000	1.0000	1.0000	−1.0000	−0.8575	0.5145	0.5831	1.7150	0.0432
2	1.2647	1.2647	−1.4706	0.8824	0.5000	0.3000	0.0000	1.0000	−1.0000	−0.8575	0.5145	0.5831	1.7150	0.0432
Pneumatic cylinder—weight
Iteration	x1	x2	u1	u2	σ1	σ2	Z(X)	∂z/∂x1	∂z/∂x2	α1	α2	σ_z	β	pf
1	0.4200	0.3700	0.0000	0.0000	0.0500	0.0400	0.0500	1.0000	−1.0000	−0.7809	0.6247	0.0640	0.7809	0.2174
2	0.3895	0.3895	−0.6098	0.4878	0.0500	0.0400	0.0000	1.0000	−1.0000	−0.7809	0.6247	0.0640	0.7809	0.2174

Table 11.

Analytical method results—Example 4.