Investigation on a novel reliability analysis approach integrating adaptive space division and direction sampling

Abstract

Multitudes of algorithms have been proposed to evaluate the failure probability of components and systems. Among all these algorithms, direction sampling is a promising one. However, a major computational effort involved in direction sampling is that for each direction sample a root-searching process needs to be conducted to obtain the distance between the origin and the failure surface, which can be computationally expensive. In addition, the failure probability along a direction is hard to obtain when the direction vector intersects the failure surface multiple times. This paper proposes a novel approach to obtain densely populated points locating on the failure surface and utilizing these points to conduct reliability analysis with direction sampling. The process of obtaining the population consists in dividing the initial hypercube space into multiple lattice grids, identifying the grids crossing the failure surface, and dividing these grids into smaller ones. The iterative division process stops when the grids cut across by the failure surface are small enough so that the centers of these grids can be considered located on the failure surface. Several numerical tests are performed to validate the applicability of the proposed approach. Through these tests, it can be concluded that for the situation where the intersection between the direction vector and the failure surface is unique, the proposed approach will lead, with high accuracy, to a failure probability. For the situation where the intersection is not unique, satisfactory estimation of failure probabilities can also be achieved using approximation methods. For high-dimensional reliability problems, clustering techniques can be utilized to reduce the computational cost.

Keywords

reliability analysis failure surface direction sampling grids division cross angle

Introduction

Structural reliability analysis is a classic research problem with a long history. The analysis result provides a sound estimation of the structural performance under load effects and possibly multiple deterioration mechanisms, which serves as the basis to lay out intervention strategies during the life-cycle of structures. By definition, the failure probability of a structure is

P_{f} = \int_{G (x) < 0} f (x) d x

(1)

where

x = {[x_{1}, x_{2}, \dots, x_{m}]}^{T}

is a vector of random variables (m is the number of random variables), f( x ) is joint probability density function of x , and G( x ) is the performance function that characterizes the limit state. Theoretically the failure probability can be directly obtained through an integration process provided that f( x ) and G( x ) are known. However, in most practical problems, G( x ) are nonlinear, rendering it difficult to pinpoint the failure region (G( x ) may even be implicit in some cases). In addition, when the number of random variables involved is large, multi-dimensional integration is computationally prohibitive. As a result, numerical integration method as presented in equation (1) is rarely carried out in real-world reliability analysis.

Different methodologies have been proposed to circumvent the integration in equation (1). Among these methodologies, approximation approach such as the First Order Reliability Method (FORM) and Second Order Reliability Method (SORM) are widely adopted, especially FORM. In FORM, the performance function is approximated at the most probable point (MPP) with first order Taylor series (i.e., linear function). MPP is a point located on the failure surface (i.e., G( x )=0) with the minimum distance from the origin in the space of standard normal variables (referred to as the U-space). In these two methods, the original random variables are converted to independent standard normal variates using Rosenblatt or Nataf transformation (Melchers, 1987). The rational for this linear approximation is that the probability densities of the random variables located in the region affected by this approximation is small. HL-RF algorithm (Fiessler et al., 1979; Hasofer and Lind, 1974) is widely used to search the MPP. SORM takes into account the curvature at the MPP, which leads to an improvement of the accuracy of reliability analysis results to some extent. In recent years new algorithms based on FORM have been developed, such as Zhong et al. (2020). However, the accuracy of these approximation approaches is dubious for the highly nonlinear and multimodal performance functions (in which multiple local MPPs exist). In addition, in the process of searching MPP, unstable phenomena such as bifurcation and chaos may occur (Lopez et al., 2015).

Sampling methods are another school-of-thought to tackle the reliability problem. The most basic sampling method is crude Monte Carlo (MC) simulation. Large numbers of samples of random variables are generated to evaluate the associated values of performance functions. The failure probability is determined as the ratio of the number of samples resulting in negative values of performance function (indicating the failure of the structure) over the total number of generated samples. For sampling methods, different trials using the same number of samples can result in different probabilities of failure. The coefficient of variation of the failure probability obtained using crude MC simulation, COV(P_f), is expressed as (Melchers, 1987)

C O V (P_{f}) = \sqrt{\frac{1 - P_{f}}{P_{f} n_{M C}}}

(2)

where n_MC is the number of generated samples. It can be seen that when P_f is small (which is the case for most of civil engineering structures), very large number of samples needs to be generated to achieve an acceptable COV. Therefore, the crude MC simulation (also referred to as brutal-force MC simulation), is not an economic approach to solve the reliability problem for structures with a small failure probability.

Quantities of smart sampling techniques have been proposed to reduce the number of samples needed to achieve a satisfactory level of accuracy, such as stratified sampling (a well-known type of stratified sampling is Latin hypercube sampling (McKay et al., 1979)), importance sampling (Tokdar and Kass, 2010), and subset simulation (Au and Beck, 2001). Song and Kawai (2023) proposed an adaptive stratified sampling method to conduct structural reliability analysis. Zuniga et al., (2021) integrated dimension reduction techniques into importance sampling approach to conduct structural reliability assessment. Although Latin hypercube sampling can reduce the number of samples, the number of samples needed are still large. For importance sampling, the region of importance needs to be identified which may be a tricky issue in itself. Subset sampling involves an iterative process which can be time-consuming although the number of samples needed is reduced. In addition, subset sampling has some limitations in terms of applicability (Breitung, 2017).

Direction sampling (Ditlevsen and Bjerager, 1986) and line sampling (Koutsourelakis et al., 2004; Schuëller et al., 2004) are two sampling methods taking the geometry of the failure surface into consideration. These two methods are also conducted in the U-space. Direction sampling approach involves uniformly generating direction samples on the unit hypersphere. The one-dimensional reliability index along each direction is associated with the distance from the origin to the intersection point between the direction vector and the failure surface. The distance can be determined by searching the root of a unit-variate equation. Line sampling, on the other hand, identifies an important direction and generates samples on the hyperplane perpendicular to the important direction. An issue with using line sampling is how to identify this importance direction. Given this issue, more robust line sampling techniques have been proposed, such as the adaptive line sampling (Pradlwarter et al., 2007) and combined line sampling (Papaioannou and Straub, 2021), among others. Recent studies have proposed new algorithms based on the direction sampling method and line sampling method, such as Cheng et al., (2023) and Jafari-Asl et al. (2021), among others. Overall, the direction sampling approach is preferred over the line sampling approach as the former skips the direction determination process which can be complicated. However, a weakness shared by both direction sampling and line sampling lies in the root-searching process. The associated one-variate equation can be highly nonlinear, and the roots can be multiple. Due to this weakness, direction sampling is widely applied only for reliability problems associated with simple-connected failure regions.

It is worth noting that using surrogate models to conduct reliability analysis is another type of approach. Surrogate models are used to represent the original performance functions that is very complex or implicit. Sampling methods are then utilized to obtain the failure probability associated with the surrogate model. Two major research aspects on the surrogate model approach are the metamodels used to fit the data and the learning function. Metamodels considered in the surrogate models includes polynomial response surface (Bucher and Bourgund, 1990; Faravelli, 1989), polynomial chaos expansions (Ghanem and Spanos, 1991), and Kriging (Vahedi et al., 2018; Wang et al., 2022), among others. Learning functions to obtain the design of experiment (the point on the exact failure surface, abbreviated as DoE) include expected feasibility function (EFF) (Bichon et al., 2008), U-function (Echard et al., 2011), and least improvement function (LIF) (Sun et al., 2017), among others. It is worth mentioning that several studies have integrated surrogate models with smart sampling approach, such as Zhang et al., (2020) and Nguyen et al., (2022), among others. The philosophy of surrogate models is to obtain an accurate reliability analysis result while minimizing the number of DoEs used to conduct the parameter fitting. As the number of DoEs is limited, there is some discrepancy between the function used in the surrogate model and the original performance function, which is the source of error associated with this type of approach. In addition, an iterative method is also needed in the process of fitting the surrogate model, which may be time-consuming.

Vast amounts of research effort have been devoted to the calculation of failure probability of structures. It should be noted that information on the configuration of failure surface is also beneficial. A primary benefit of knowing the configuration of the failure surface is that when the values of some random variables are identified through structural health monitoring (SHM) techniques, the failure probability can be updated with ease without carrying out new reliability analyses. For instance, for a bridge structure subjected to dead and live loads, the uncertainties associated with the dead loads, live loads, and resistances are involved in the reliability analysis processes. If the failure surface is identified, the time-variant failure probability is known instantly if live load effects are measured using weigh-in-motion (WIM) systems. For a steel structure subjected to corrosion-enhanced fatigue, the fatigue reliability involves the uncertainties associated with both corrosion modelling and fatigue modelling. The fatigue reliability can be updated with ease when the stress histogram is obtained using SHM if the failure surface of the structure subjected to corrosion-enhanced fatigue is known. In addition, failure surface can help pinpointing the critical regions of values of random variables that lead to the failure of the structure, such as the region around the most probable point (MPP) which is sought in FORM.

This paper proposes a novel approach to tackle the reliability problem, which is to obtain densely populated points on the failure surface and use them to conduct reliability analysis with the direction sampling technique. For each generated direction sample, the point to which the direction from the origin in the U-space is closest to the direction sample is selected. The distance from the point to the origin in the U-space is considered as the reliability index along the direction of the sample if the intersection point is unique. For the case where a direction vector has multiple intersection points with the failure surface, an approximation method is proposed to obtain all the intersection points. In addition, clustering techniques can be utilized to reduce the computational expense when the dimensionality of the reliability problem is high. The novelty of this approach lies in using the densely populated points to infer the intersections between direction vectors and failure surfaces. Using this approach to solve the reliability problem associated with multiple intersections between the direction vectors and the failure surface as well as integrating clustering techniques into this approach further enhances the novelty of this study. Several examples are investigated to verify the validity of this approach. The flowchart of the proposed approach is shown in Figure 1.

Figure 1.

Flowchart of the proposed reliability approach.

Theoretical background of direction sampling

In the space of the independent standard normal random variables (i.e., the U-space), the shape of the joint PDF ϕ_u( u ) is rotationally symmetric with respect to the origin. The PDF of ϕ_u( u ) is largest at the origin and decreases drastically along each direction. This characteristic renders the application of the polar (spherical) coordinate system in the failure probability analysis. In fact, using polar coordinate system may be more convenient than Cartesian coordinate system in solving the reliability problem.

An m -dimensional independent standard normal vector U in the U-space can be expressed as U = R A , where R is the Euclidean norm of the vector (i.e., R=‖ U ‖) and A is the random unit direction vector (i.e., A = U /‖ U ‖).

The failure probability in equation (1) can be rewritten as

P_{f} = \int_{Ω_{A}} \int_{0}^{+ \infty} I_{f} [T^{- 1} (r α)] f_{R} (r) f_{A} (α) d r d α

(3)

where T⁻¹(⋅) is the inverse normal transformation (transforming the independent normal random variables in the U-space into original random variables X ); Ω_A is an m−1 dimensional hypersphere with a unit radius; I_f( X ) is an indicator function equal to one when G( X )≤0 and 0 when G( X )>0.

If Monte Carlo simulation is conducted using equation (3), the computation effort will be similar to that using samples generated based on Cartesian coordinates. An alternative methodology based on the concept of conditional probability was proposed in Ditlevsen and Bjerager (1989), which is expressed as

P_{f} = \int_{α {\in Ω}_{A}} P (F | α) f_{A} (α) d α = \int_{α {\in Ω}_{A}} {1 - χ_{m}^{2} [r {(α)}^{2}]} f_{A} (α) d α

(4)

In which

P (F | α)

is the probability of failure in the direction of α ;r( α ) is the distance between the failure surface and the origin in the direction of α , and

χ_{m}^{2} [\cdot]

is the chi-square distribution function of m degrees of freedom. Equation (4) is based on the fact that the distance from origin r conforms to the chi-square distribution with m degrees of freedom in the m-dimensional U-space. Assuming uniform sampling on an m−1 dimensional sphere, Monte Carlo simulation with samples of α can be conducted using equation (4) to obtain failure probability:

P_{f, d i r} = \frac{\sum_{i = 1}^{N_{d}} {1 - χ_{m}^{2} [r {(α_{i})}^{2}]}}{N_{d}}

(5)

where N_d is the number of direction samples and α _i is the i th direction sample. The method proposed in Muller (1959) and Marsaglia (1972) can be adopted to uniformly generate samples on the surface of an m−1 dimensional sphere, in which samples are generated by dividing an m-dimensional vector by its norm. It should also be noted that equations (4) and (5) only apply for the case where r( α _i) is unique. Assuming α _i intersects the failure surface twice, the distances associated with two intersection points are referred to as r_i1 and r_i2, respectively, in that case

P (F | α) = χ_{n}^{2} (r_{i 2}^{2}) - χ_{n}^{2} (r_{i 1}^{2})

(assuming r_i2>r_i1 ). Therefore, when the intersection points between the direction vector and the failure surface are multiple, all the intersection points need to be found out to accurately predict the failure probability in this direction. However, even in the case where only one intersection point exists, finding r( α _i) is not an easy task, as the equation having r( α _i) can be highly nonlinear. The arduous process of finding r( α _i) is a major drawback of the direction sampling method.

Novel approach to conduct reliability analysis

The key concept of the novel approach to conduct reliability analysis is to obtain densely populated points on the failure surface in the U -space. When the points are obtained, for each direction vector pointing from the origin, it is possible to find the point very close to vector and consider this point as the intersection between the direction vector and the failure surface, thereby eliminating the root-searching process. The random variables that are not standard normal can be converted into standard normal random variables using Rosenblatt or Nataf transformation (Melchers, 1987). The approach is divided into two parts: (i) obtaining the lattice grids crossing the failure surface that are small enough so that the center of gravity of these grids can be considered as locating on the failure surfaces, and (ii) utilizing direction sampling to obtain the failure probability.

Division of lattice grids

The procedure of obtaining the small lattice grids crossing the failure space is as follows:

• Step 1: Identify a hypercube space in the m -dimension U -space. The hypercube space is defined as $(U_{l 1}, U_{u 1}) \times \dots (U_{l i}, U_{u i}) \dots \times (U_{l m}, U_{u m})$ , where U_li and U_bi (i=1,…,m) are the lower and upper bounds associated with the i th dimension, respectively. Generate lattice points within the initial hypercube. Lattice points are uniformly spaced at an interval of U_int,i in the i th dimension. An illustration in this regard is provided in Figure 2. Using the lattice points, the initial hypercube space is divided into multiple lattice grids.

• Step 2: Determine if a grid crosses the failure surface. This process is conducted by using the values of the performance function at the vertices of the grid and randomly generated points within the grid.

• Step 3: Ignore the grids that do not cross the failure surface. Divide all the grids crossing the failure surface into smaller grids. The division is performed by dividing the segment of a grid into two parts in each dimension. After the division, each new grid is associated with one-half of the length of the divided grid in each dimension (i.e., U_int,in = U_int,i/2, where U_int,in is the new grid interval in the i th dimension).

• Step 4: Return to step 2 unless all the grids crossing the failure surface are small enough.

Figure 2.

Illustration of a three-dimension cube divided using lattice points.

The lower and upper bounds of the hypercube in each dimension (U_li and U_bi, i = 1,…,m) and the length of lattice grid in each dimension (U_int,i; i = 1,…,m) are user-specified. Theoretically, a larger initial hypercube can help capture the failure surface more thoroughly, especially the regions in the failure surface that are distant from the origin. However, as the contribution of the far-away regions in the failure surface to the failure probability is small, capturing these regions may not improve the accuracy of the calculated failure probability significantly. In addition, a large initial space may incur a heavy burden of computational cost. Therefore, a balance between the accuracy and the computational cost should be made by the users regarding the initial hypercube space. The initial U_int,i needs not to be small, although a smaller U_int,i can also help characterizing the failure surface more accurately.

It is self-evident that a grid crossing the failure surface is split into two parts by this surface. Therefore, points inside the grid crossing the failure surface can render the value of the performance function both positive and negative. On the other hand, for a grid not crossing the failure surface, the points inside it can only make the value of the performance function positive or negative. A grid can be considered as crossing the failure surface when the points within it are associated with both positive and negative values of the performance function. A simple approach to determine if a grid crosses the failure surface is to evaluate the performance function at the vertices of the performance functions. If both positive and negative values are obtained, then it can be determined that the grid crosses the failure surface. However, it should be noted that if all the values of the performance function at vertices are positive/negative, it is still possible that the failure surface crosses the grid. Therefore, sample points inside the grid need to be generated to determine if a grid crosses the failure surface. It is acknowledged that adding points inside the grid may still fail to determine if a grid crosses the failure surface correctly, although this possibility decreases with the increase of number of sample points. The number of samples should be related to the size of the grid. The larger the size of a grid is, the larger the number of samples should be adopted. For the value of the number of samples generated within a grid N_s, a balance should be achieved between the probability of accurate detection and the acceptable computational cost.

Apparently, a specific criterion should be set up to define the smallest grid. For a m-dimension grid, the criteria should be related to the length of the grid in each dimension (i.e., U_int,i;i = 1,…,m). Assuming the lengths of grid in each dimension are equal (U_int,i = U_int,j = U_int, i≠j), the iterative process stops when U_int≤ϵ. The smaller ϵ is, the more accurately the failure surface can be represented using the centers of the grids. As a smaller ϵ also implies a larger computational expense, a tradeoff between the accuracy in portraying the failure surface and computational cost should be achieved. In this paper, when the lengths of a grid in each dimension are equal, the number of samples N_s used to determine if a grid intersects the failure surface is expressed as

N_{s} = \max (N_{r} \cdot U_{int}, 1)

(6)

where N_r is the number of samples per unit length of a grid in each dimension. Equation (6) ensures that at least one sample is generated even if the length of a grid in each dimension is very small.

It is worth noting that for system reliability analysis involving multiple performance functions, an equivalent performance function can be derived based on the system model. Assuming n_g performance functions $g_{1}, g_{2}, \dots, g_{n_{g}}$ are involved, the equivalent performance function for a series system is

g_{e, s} = \min_{i} g_{i}, i = 1, \dots, n_{g}

(7a)

For a parallel system, the equivalent performance function is

g_{e, p} = \max_{i} g_{i}, i = 1, \dots, n_{g}

(7b)

Direction sampling using the center of grids

After the small grids crossing the failure surface are obtained, the centers of these grids are used in the direction sampling process. Depending upon the number of intersections between a direction vector and the failure surface, the following two cases are discussed separately.

• Case A: The failure surface and the direction vector only intersect at one point.

In this case, as the centers of grids are densely populated on the failure surface, one of these centers can be considered as the intersection point between the direction vector and the failure surface. For a center of grid, a vector v can be formed by pointing from the origin to w . For a specific direction vector, the center of gravity, w _r, the vector of which having the smallest cross angle with this direction vector is selected. The cross angle θ between the two vectors, v ₁ and v ₂, can be determined as

θ = \arccos (\frac{v_{1} \cdot v_{2}}{‖ v_{1} ‖ \cdot ‖ v_{2} ‖})

(8)

Let the smallest cross angle obtained be θ_s. If θ_s is larger than a specific threshold (referred to as θ_s,th ), the direction vector can be considered as not crossing the failure surface and the failure probability along this direction is zero. Otherwise, the failure probability along the direction vector is (Melchers, 1987)

P_{f, d i r} = {\begin{cases} 1 - χ_{m}^{2} ({‖ w_{r} ‖}^{2}) if \underset{r \to 0}{g} > 0 \\ χ_{m}^{2} ({‖ w_{r} ‖}^{2}) if \underset{r \to 0}{g} < 0 \end{cases}

(9)

where

\underset{r \to 0}{g}

is the value of the performance function when the distance r between an intersection point and the origin in the U-space is close to 0. The value of θ_s,th should be neither too small nor too large. A very small θ_s,th may result in a failure probability of zero when the direction vector actually intersects the failure surface. A very large θ_s,th may lead to a severe enlargement of the real failure surface. In fact, the value of θ_s,th is dependent upon the value ϵ. The smaller ϵ is, the smaller the value of θ_s,th can be. Before setting up θ_s,th, the cross angle between the vectors associated with two centers of adjacent grids can be calculated. θ_s,th can be set a value several times larger than this cross angle.

• Case B: The failure surface and the direction vector intersect at multiple points.

The problem becomes much more complex when multiple intersection points exist between the direction vector and the failure surface. The reason for the existence of multiple intersection points is that either the failure region is enclosed or there are holes in the failure region. An illustration of this scenario is shown in Figure 3(a). Depending upon the sign of $\underset{r \to 0}{g}$ and $\underset{r \to \infty}{g}$ (the value of the performance function when r is approaching infinity), a total of four scenarios are involved:

Figure 3.

(a) Illustration of multiple intersection points between a direction vector and failure surface; (b) illustration of two close intersection points between a direction vector and failure surface.

Scenario B1: $\underset{r \to 0}{g} > 0$ and $\underset{r \to \infty}{g} < 0$

In this scenario, the intersection points between a direction vector and the failure surface are odd-numbered. Assuming the distances from the origin to the failure surface along the direction vector are r₁,r₂,…,r_2k+1(k = 1,2,…) in the ascending order, the failure probability along the direction vector is

P_{f, d i r} = 1 - χ_{m}^{2} (r_{2 k + 1}^{2}) + \sum_{i = 1}^{k} χ_{m}^{2} (r_{2 i}^{2} - r_{2 i - 1}^{2})

(10a)

Scenario B2: $\underset{r \to 0}{g} > 0$ and $\underset{r \to \infty}{g} > 0$

In this scenario, the intersection points between a direction vector and the failure surface are even-numbered. Assuming the distances from the origin to the failure surface along the direction vector are r₁,r₂,…,r_2k(k =1,2,…) in the ascending order, the failure probability along the direction vector is

P_{f, d i r} = \sum_{i = 1}^{k} χ_{m}^{2} (r_{2 i}^{2} - r_{2 i - 1}^{2})

(10b)

Scenario B3: $\underset{r \to 0}{g} < 0$ and $\underset{r \to \infty}{g} > 0$

P_{f, d i r} = χ_{m}^{2} (r_{1}^{2}) + \sum_{i = 1}^{k} χ_{m}^{2} (r_{2 i + 1}^{2} - r_{2 i}^{2})

(10c)

Scenario B4: $\underset{r \to 0}{g} < 0$ and $\underset{r \to \infty}{g} < 0$

In this scenario, the intersection points between a direction vector and the failure surface are even-numbered. Assuming the distances from the origin to the failure surface along the direction vector are r₁,r₂,…,r_2k(k = 1,2,…) in the ascending order, the failure probability along the direction vector is

P_{f, d i r} = χ_{m}^{2} (r_{1}^{2}) + 1 - χ_{m}^{2} (r_{2 k}^{2}) + \sum_{i = 1}^{k - 1} χ_{m}^{2} (r_{2 i + 1}^{2} - r_{2 i}^{2})

(10d)

The methodology of finding multiple r is detailed as follows:

• Step 1: Identify the center of grid associated with the smallest cross angle with the direction vector. This center of grid, w _cr, can be considered as one of the intersections between the direction vector and the failure surface. The method of obtaining this center of gravity has been indicated in Case A.

• Step 2: For all the vectors v _j (j =1,…,N_c−1) associated with the rest of centers of grids w _j (j =1,…,N_c−1), identify all the v _j with a cross angle associated with the direction vector smaller than θ_th. The set of centers of grids associated with these v _j is denoted as $W$ . All the centers of grids in $W$ can be approximately considered in the same direction as w _cr.

• Step 3: Among all the centers of grids in $W$ , find all the centers of grids w _f ( f =1,…,N_f, where N_f is the number of the centers that meets the criteria) that satisfy ‖ w _f− w _cr‖>r_th. The rationale is that among all the centers of gravity in $W$ with a distance from w _cr larger than a certain value, at least one of the points must be another intersection point between the direction vector and the failure surface.

• Step 4: Consider these w _f as the new $W$ . The center of grids in $W$ having the smallest distance from the origin is selected as the next intersection point and named as the new w _cr. Remove w _cr from $W$ .

• Step 5: Repeat step 3 until w _f is an empty set.

• Step 6: Sort all the intersection points identified based on their distances from the origin. The reason to conduct the sorting is that the first identified w _cr may not be the closest intersection point from the origin in the U-space.

The value of θ_th should be small enough so that all the centers of grids that are approximately in the same direction as the direction vector can be identified. The value of θ_th is also dependent upon the value of ϵ. In fact, θ_th can be set equal to θ_s,th. The value of r_th should be large enough so that centers of gravity that are near the selected intersection point and yet not the real intersection points can be filtered out. Based on numerical experiences, r_th should be several times larger than θ_s,th. It should be noted that there are cases where the direction vector is almost tangent to the failure surface (an illustration in this regard is provided in Figure 3(b)). In this case, some intersection points may fail to be detected by the methodology indicated previously. Two different approaches are proposed when this scenario occurs. The first, denoted as approximation approach I, considers the failure probability associated with the direction vector of which the intersection points fail to be detected in full as zero. The other, denoted as approximation approach II, considers the failure probability associated with the original direction vector to be the same as the failure probability associated with a substitute direction vector. The substitute direction vector is selected among all the direction vectors of which the intersection points are successfully identified as the one with the smallest cross angle with the original direction vector. In Figure 3(b), assuming the intersection points between direction vector v ₁ and the failure surface fail to be identified in full while all the intersection points between direction vector v ₂ and the failure surface are successfully identified, the failure probability associated with v ₂ may be used as the failure probability associated with v ₁.

Criteria can be proposed regarding if the intersection between a direction vector and the failure surface is unique or not. For a series system reliability problem, it can be asserted that the intersection between a direction vector and the failure surface is unique. For a component reliability problem of which the performance functions involve nonperiodic functions, such as polynomial performance functions, it is likely that intersection between a direction vector and the failure surface is unique. For a component reliability problem associated with performance functions involving periodic functions such as trigonometric functions and parallel system reliability problem, the intersection between a direction vector and the failure surface may not be unique. Monotonicity analysis of the performance functions can be carried out to determine if the intersection point is unique or not.

Using clustering algorithms to reduce the computational cost

It is acknowledged that as the dimension of the performance function increases, the size of the matrix containing the centroids of the grids can increase substantially. Storing such matrices when the grid size is small requires large computer memory and therefore represents a major obstacle for the implementation using the proposed algorithm. For the case where the intersection between a direction vector and the failure surface is unique, clustering techniques can be used to reduce the computational cost and expedite the calculation process. By using clustering techniques, the centroids of grids crossing the failure surface are divided into multiple clusters and only the clusters of grids for which the direction is close to the generated direction samples are selected and continue to be divided.

The procedure of using clustering algorithms to accelerate the calculation process is as follows:

• Step 1: Divide the grids crossing the failure surface until the length of the grid U_int = ϵ₁ (assuming the length of the grid in each dimension is identical).

• Step 2: Separate the grids into N_k clusters. For each cluster, calculate the coordinate of each centroid. This coordinate can be considered as a vector, with the origin and the centroid considered as the initial point and terminal point of a vector, respectively. The set of the vectors representing the centroids of all the clusters is denoted as $C$ .

• Step 3: For a specific direction sample, calculate cross angles of this direction sample with each vector in $C$ . Sort the vectors consisting of cross angles.

• Step 4: Select N_ks clusters of grids based on the sorting results, i.e., only the clusters associated the first N_ks cross angles in the ascending order are selected.

• Step 5: Continue dividing the girds in the selected clusters until the final size threshold ϵ is reached.

• Step 6: Compute the failure probability along this direction using equations (8) and (9).

• Step 7: Repeat Steps 1–6 for all direction samples.

A flowchart of this procedure is shown in Figure 4. Clustering techniques aim to divide a cloud of spatial points into multiple clusters, with the purpose of minimizing the sum of variances of clusters (the variance of a cluster is the sum of square of point-to-cluster-centroid distances) (Du et al., 1999). K-means clustering algorithm (Hartigan and Wong, 1979), a widely used machine learning algorithm, is adopted herein to implement the clustering process. Theoretically, the larger the number of selected clusters N_ks is, the higher the accuracy of the proposed approach is. However, a larger N_ks also indicates a larger computational cost. Therefore, the value of N_ks is user-specified, depending upon the acceptable computational cost as well as the required accuracy of the users. The number of total clusters N_k is also user-specified. When N_ks is fixed, a larger N_k signifies a lower computational cost. However, a larger N_k corresponds to a longer duration to complete the clustering process.

Figure 4.

Flowchart of the proposed approach integrating clustering techniques.

Illustrative examples

To check the applicability of the proposed reliability analysis approach, six example problems are investigated herein. All the example problems investigated herein were considered in previous studies, where reliable results of reliability analysis were provided. To verify the results obtained using this approach, four different methods are considered herein, namely

• Method A: The approach proposed herein with 10⁴ direction samples.

• Method B: Crude Monte Carlo simulation with 10⁶ samples.

• Method C: First Order Reliability Method (FORM)

• Method D: Adaptive Kriging Monte Carlo Method (AK-MCS) (Echard et al., 2011)

Both Methods B and C are widely-used approaches to obtain failure probability. Method D is a metamodel-based approach combining Monte Carlo simulation with Kriging metamodels. Comparison between the results of Method A with those associated with Methods B, C, and D can indicate whether the result of the approach proposed herein is reliable.

It should be noted that FORM can only obtain the reliability index of individual components. For system reliability, FORM needs to be integrated with the cumulative distribution function (CDF) of multi-variate normal distribution. For a series system of γ failure modes, the failure probability is calculated as (Hohenbichler and Rackwitz, 1982)

P_{f, s y s} = 1 - \int_{- \infty}^{β_{1}} \dots \int_{- \infty}^{β_{γ}} \frac{1}{\sqrt{{(2 π)}^{γ} | Σ |}} \exp (- \frac{1}{2} z^{T} Σ^{- 1} z) d z_{1} \dots d z_{γ}

(11a)

For a parallel system of γ failure modes, the failure probability is calculated as (Hohenbichler and Rackwitz, 1982)

P_{f, s y s} = \int_{β_{1}}^{\infty} \dots \int_{β_{γ}}^{\infty} \frac{1}{\sqrt{{(2 π)}^{γ} | Σ |}} \exp (- \frac{1}{2} z^{T} Σ^{- 1} z) d z_{1} \dots d z_{γ}

(11b)

where β_h(h=1,…,γ) is the reliability index associated with the h th failure mode;

z = {[z_{1}, \dots, z_{γ}]}^{T}

is the γ -dimensional vector of standard normal random variables; Σ is the correlation matrix indicating the correlation among different failure modes. Detailed information is provided in Hohenbichler and Rackwitz (1982).

For Method D, the learning function and the convergence criteria are based on Echard et al., (2011). The number of samples in the initial DoE and the number of samples added to the experimental design are 10 and 1000, respectively. The number of Monte Carlo samples in Method D is 10⁵. In the last example, reliability analysis integrating clustering techniques is presented.

Example 1: Component reliability problem with linear performance functions

The linear performance functions investigated in Katsuki and Frangopol (1994) are adopted as the first example. These performance functions are:

g_{1} (u_{1}, u_{2}) = - u_{1} - u_{2} + 3 \sqrt{2}

(12a)

g_{2} (u_{1}, u_{2}, u_{3}) = - u_{1} - u_{2} - u_{3} + 3 \sqrt{3}

(12b)

g_{3} (u_{1}, u_{2}, u_{3}, u_{4}) = - u_{1} - u_{2} - u_{3} - u_{4} + 3 \sqrt{4}

(12c)

g_{4} (u_{1}, u_{2}, u_{3}, u_{4}, u_{5}) = - u_{1} - u_{2} - u_{3} - u_{4} - u_{5} + 3 \sqrt{5}

(12d)

where u₁,u₂,u₃,u₄, and u₅ are independent standard normal random variables. These four performance functions are associated with the same failure probability, Φ(−3)=1.350×10⁻³, where Φ(⋅) is the CDF of standard normal distribution.

The three-dimension (3D) case (i.e., g₂ ) is used to illustrate the procedures of calculating failure probability in Method A. For each dimension of the hypercube, the lower bound U_l, upper bound U_b, and the initial size the of grid U_int are set as 0, 4, and 2, respectively. N_r is set as one while ϵ is 0.1. The initial grid plot, the grids crossing the failure surface after the first grid division, and the centers of final grids crossing the failure surface are shown in Figure 5. It can be seen that the linear failure surface can be represented by densely populated centers of grids crossing the failure surface. Apparently, a direction vector can only intersect the failure surface once for this problem. Setting θ_s,th as 0.1, the failure probabilities using Method A associated with these four performance functions and the relative errors (defined as $\frac{| P_{f A} - Φ (- 3) |}{Φ (- 3)}$ , where P_fA is the failure probability calculated using Method A) are shown in Table 1. It can be seen that the relative error obtained using Method A is generally small (less than 10% in this case), which indicates the analysis using the approach proposed herein is reliable.

Figure 5.

Illustrations of grid division process for the 3D case in Example 1. (a) Initial grid plot. (b) Grids crossing failure surface after first division. (c) Centers of final grids crossing failure surface.

Table 1.

Probability of failure associated with Example 1.

Performance function	Failure probability P_fA (×10⁻³ )	Relative error (%)
g ₁	1.430	5.9
g ₂	1.411	4.5
g ₃	1.320	2.2
g ₄	1.321	2.1

Note: Relative error is defined as |p_fA−Φ(−3)|/Φ(−3).

An analysis is carried out using the 2D case to illustrate the effect of number of samples to conduct the simulation associated with Methods A and B on the failure probability results. The relationship between the number of samples used and the obtained failure probability is shown in Figure 6. For this specific case, it can be seen that the dispersion of failure probability is small after the number of samples reaches 10⁴ in Method A, which can prove that the number of samples used herein for Method A (i.e., 10⁴) is sufficient to obtain an accurate failure probability. On the other hand, if Method B is used, the failure probability stabilizes after the number of samples reaches 2 × 10⁵. It can be observed that for the problem associated with a small failure probability, using direction sampling leads to a reduced number of samples compared with using crude MC sampling.

Figure 6.

Effect of the number of samples used in Methods A and B on the failure probability. (a) Method A, (b) Method B.

Example 2: Component reliability problem with nonlinear performance function

The second example is associated with a nonlinear performance function, which is expressed as

G (u_{1}, u_{2}) = \exp (0.4 u_{1} + 7) - \exp (0.3 u_{2} + 5) - 200

(13)

where u₁ and u₂ are independent standard normal random variables. This problem was investigated in Roussouly et al., (2013) using the response surface method. Based on the analysis result in Roussouly et al., (2013), the failure probability for this problem is around 3.6 × 10⁻³. Method B is expected to provide an accurate result for this problem.

For this problem, the initial space is set as [−4, 4]×[−4, 4]. U_int and ϵ are set as 2.5 and 0.1, respectively. N_r and θ_s,th are set as 1 and 0.1, respectively. The grids crossing the failure surface after the first division and the centers of the final grids crossing the failure surface are plotted in Figure 7. The failure and safe regions associated with this performance function are indicated in Figure 7(b). In this case, each direction vector has only one intersection with the failure surface. The failure probability associated with Method A is 3.674 × 10⁻³. The failure probabilities associated with Methods B, C, and D are 3.672 × 10⁻³, 3.365 × 10⁻³, and 3.77 × 10⁻³, respectively. It can be seen that the failure results obtained by these Methods A and B are close to each other. The result obtained by Method A is more accurate than those obtained by Methods C and D.

Figure 7.

Grids after the first division and the center of final grids associated with Example 2.

A parametric analysis using Method A is conducted on θ_s,th with this problem. The result is shown in Figure 8. It can be seen that for the case where the intersection point between a direction vector and failure surface is unique, if θ_s,th is too small (less than 10⁻² in this case), the intersection points of a direction vector cannot be identified. As a result, the failure probability obtained is smaller than the true failure probability. If θ_s,th is larger than a threshold (10⁻² in this case), accurate failure results can be obtained. Although a large θ_s,th (e.g., larger than 0.1) may result in a nonzero failure probability associated with a direction vector not crossing the failure surface, as these nonzero failure probabilities are very small, the influence of further increasing θ_s,th when it surpasses a threshold (10⁻² in this case) on the overall failure probability considering all the direction vectors is minimal.

Figure 8.

Influence of θ_s,th on the failure probability associated with Example 2 using Method A.

Example 3: System reliability for a series system

The third example is associated with a series system consisting of four components. The equivalent performance function is expressed as

G (x_{1}, x_{2}) = \min {\begin{array}{c} 3 + 0.1 {(x_{1} - x_{2})}^{2} - \frac{x_{1} + x_{2}}{\sqrt{2}} \\ 3 + 0.1 {(x_{1} - x_{2})}^{2} + \frac{x_{1} + x_{2}}{\sqrt{2}} \\ (x_{1} - x_{2}) + \frac{6}{\sqrt{2}} \\ (x_{2} - x_{1}) + \frac{6}{\sqrt{2}} \end{array}}

(14)

where x₁ and x₂ are standard normal random variables. Three different correlation cases are explored to shed light on the influence of correlation among random variables on the system reliability. In the first correlation case, x₁ and x₂ are independent (i.e., correlation coefficient ρ is 0); in the second and third correlation cases, x₁ and x₂ are correlated with a correlation coefficient ρ of 0.5 and 0.9, respectively. The problem associated with the first correlation case was investigated in Sun et al., (2017). According to Sun et al., (2017), the failure probability associated with considering independent random variables is about 4.40 × 10⁻³.

As indicated previously, for a series system, a direction vector can only intersect the failure surface once. The initial space and the values of parameter of interest are the same as those in Example 2. Scatter plots of the centers of the final grids crossing the failure surface associated with the three correlation cases in the U-space are shown in Figure 9, where failure and safe regions are indicated. It can be seen that correlation among random variables makes a difference in the shape of failure surface. The failure probabilities associated with four different methods are presented in Table 2. For the first correlation case (ρ = 0), the failure probabilities associated with Methods A and B match the result given in Sun et al., (2017). Considering all the three correlation cases, it can be seen that Methods A, B, and D give similar reliability analysis results. When the failure probability is not very small (larger than 10⁻⁴), using Monte Carlo simulation is expected to provide an accurate failure probability. Therefore, Methods A, B, and D can provide an accurate result for this problem. As near perfectly negative correlations exist among safety margins of some components in this series system, the failure probability obtained using Method C has a considerable error. It can be seen that the correlation among random variables in the performance function does not affect the applicability of the approach proposed herein to the system reliability analysis for series systems.

Figure 9.

Centers of final grids crossing the failure surface associated with different correlation coefficients ρ between the random variables in Example 3. (a) ρ = 0, (b) ρ = 0.5, (c) ρ = 0.9.

Table 2.

Probability of failure associated with Example 3.

Method	Failure probability P_f (×10⁻³)
Method	ρ = 0.0	ρ = 0.5	ρ = 0.9
A	4.632	11.800	28.579
B	4.491	11.944	28.777
C	5.402	14.293	29.524
D	4.37	12.05	28.49

Note: ρ is the correlation coefficient among the two random variables.

Example 4: System reliability for a parallel system

The fourth example is associated with a parallel system consisting of three components. The equivalent performance function is expressed as

G (x_{1}, x_{2}) = \max {\begin{array}{c} \frac{(x_{1}^{2} + 4) (x_{2} - 1)}{20} - \sin (\frac{5 x_{1}}{2}) - 2 \\ {(x_{1} + 2)}^{4} - x_{2} + 4 \\ {(x_{1} - 4)}^{3} - x_{2} + 2 \end{array}}

(15)

where x₁ has a normal distribution with a mean of two and a standard deviation of 1; x₂ has also a normal distribution with a mean of five and a standard deviation of 1. For this problem, it is possible that a direction vector intersects the failure surface more than once. The three correlation cases assumed in the previous example are also considered for this problem, in which a correlation coefficient ρ of 0, 0.5, and 0.9 is used to characterize the correlation between x₁ and x₂. This problem with ρ = 0 was investigated in Bichon et al., (2011) using surrogate models. According to Bichon et al. (2011), the failure probability for this problem associated with ρ = 0 is about 1.15 × 10⁻³.

For the first and second correlation cases (i.e., ρ = 0 and = 0.5), the initial space is set as [−5, 5]×[−5, 5]. For the third correlation case (i.e., = 0.9), as the failure probability is very small, a larger initial space needs to be set to capture the failure surface. Through a trial-and-error process, the initial space for the third correlation case is determined as [−10, 10]×[−10, 10]. For all the three correlation cases, ϵ is set as 0.1. Taking the first correlation case as an example, three different values of N_r, 1, 5, and 10, are set to investigate the influence of the number of generated samples on the failure probability. The centers of the final grids crossing the failure surface associated with these three N_r values are shown in Figure 10. It can be seen in Figure 10(a) that a small N_r may result in a failure to identify part of the failure surface. Therefore, for performance functions associated with an enclosed failure region, increasing N_r may help capture the failure surface more accurately. When N_r increases to a certain level, further increasing N_r may not improve the characterization of failure surface significantly, as can be seen from Figure 10(b) and Figure 10(c) (in which the failure and safe regions are indicated).

Figure 10.

Centers of final grids crossing the failure surface associated with different numbers of generated samples in a grid N_r in Example 4 with ρ = 0. (a) N_r = 1, (b) N_r = 5, (c) N_r = 10.

In the direction sampling process, θ_s,th and θ_th are both set as 0.01 while r_th is set as 0.1. A smaller θ_s,th is set for this problem compared with the previous examples to pinpoint all the centers of grids that are approximately in the same direction as the direction vector in a more accurate manner. It can be determined in this case that the number of intersection points between a direction vector and failure surface is even-numbered. The failure probability results associated with three correlation cases using two approximation approaches are shown in Table 3. For the results associated with the first correlation case (ρ = 0), the failure probabilities associated with Methods A, B, and D are close to that given in Bichon et al. (2011), and therefore, can be considered as acceptable results. As FORM fails to determine the first performance function involving trigonometric functions in equation (15), the failure probability associated with Method C is off the mark. For the second correlation case (ρ = 0.5), the failure probability associated with Method B can still be considered as an accurate result. For this case, the failure probabilities associated with Methods A and B are deemed as acceptable. For the third correlation case (ρ = 0.9), as the failure probability is extremely small, Methods B, C, and D fail to give a reasonable result. To verify the results given by Method A, a fifth method (Method E) using importance sampling is adopted to obtain the failure probability associated with a correlation coefficient = 0.9. For the importance sampling, the design point in the U-space is (−4, 8), with 10⁵ samples in total and 10⁴ samples in a batch. The failure probability associated with Method E is calculated as 1.116 × 10⁻¹⁰, which is closer than that given by Method A. Therefore, the correlation coefficient among random variables does not influence the applicability of the approach proposed herein to system reliability analysis for parallel systems.

Table 3.

Probability of failure associated with Example 4.

Method	Failure probability P_f
Method	ρ = 0.0	ρ = 0.5	ρ = 0.9
A (approximation approach I)	8.957 × 10⁻⁴	1.457 × 10⁻⁴	6.346 × 10⁻¹¹
A (approximation approach II)	1.138 × 10⁻³	1.689 × 10⁻⁴	1.123 × 10⁻¹⁰
B	1.152 × 10⁻³	1.770 × 10⁻⁴	0
C	1.110 × 10⁻¹⁴	0	0
D	1.13 × 10⁻³	2.0 × 10⁻⁴	0

Note: ρ is the correlation coefficient among random variables; for the correlation case ρ = 0.9, the number of samples used in Methods B, C, and D is insufficient to generate a sample of failure, therefore, the failure probability associated with these three methods is zero.

A parametric analysis on θ_s,th in Method A is carried out for this problem associated with the first correlation case (ρ = 0). θ_th is considered equal to θ_s,th and r_th is equal to 10θ_th. The relationship between failure probability given by Method A and θ_s,th is shown in Figure 11. It can be seen that a very small θ_s,th (smaller than 10⁻² in this case) results in an underestimation of the failure probability while a large θ_s,th (larger than 0.05) results in an overestimation of the failure probability. Failure probabilities associated with a θ_s,th of 0.01 (the value used herein), is025, and 0.1 are indicated in Figure 11. The value of θ_s,th can be chosen from a range associated with the transition of the slope of the curve indicating the relationship between the failure probability and θ_s,th. In this case, θ_s,th between 0.01 and 0.025 is considered as acceptable.

Figure 11.

Influence of θ_s,th in Method A on the failure probability associated with Example 4 with ρ = 0.

Example 5: Reliability analysis associated with multiple disconnected failure regions

For direction sample technique, the trickiest problem is the one associated with multiple disconnected failure regions (e.g., a performance function with multiple “holes” in its failure region). In this regard, reliability analysis associated with modified Rastrigin function (Echard et al., 2011) is investigated herein. The performance function is

G (u_{1}, u_{2}) = 10 - \sum_{i = 1}^{2} (u_{i}^{2} - 5 \cos (2 π u_{i}))

(16)

where u₁ and u₂ are independent standard normal variables. This problem was used in Echard et al., (2011) to test the applicability of several reliability algorithms. Based on Echard et al., (2011), the failure probability associated with this problem is around 7.3 × 10⁻².

A 3D surface plot on the values of the performance function is shown in Figure 12. It can be seen that the function is associated with multiple local maxima and minima and intersects with the plane G(u₁,u₂) = 0 at multiple regions. As the modified Rastrigin function is symmetric with respect to both u₁-axis and u₂-axis, the failure surface in the first quadrant of the U-space is investigated and the direction sampling is conducted with the centers of grids in the first quadrant. The initial space is set as [0, 5]×[0, 5]. The centers of the final grids crossing the failure surface using N_r = 1, 5, and 50 samples associated with ϵ = 0.1 are shown in Figure 13. It can be seen that to fully capture the failure surface with multiple small holes, a relatively large N_r is needed. The influence of ϵ on the failure surface characterization is also investigated using this example. Using N_r = 50, the centers of grids crossing the failure surface associated with ϵ = 0.05 and ϵ = 0.01 are shown in Figure 14(a) and 14(b), respectively. It can be seen that a smaller ϵ leads to a more accurate characterization of the failure surface. The safe region associated with this performance function is indicated in Figure 14(b).

Figure 12.

3D surface plot of the modified Rastrigin function.

Figure 13.

Centers of final grids crossing the failure surface associated with different numbers of generated samples N_r with an ϵ of 0.1 in Example 5. (a) N_r = 1, (b) N_r = 5, (c) N_r = 50.

Figure 14.

Centers of final grids crossing the failure surface associated with different final grid sizes ϵ with a N_r of 50 in Example 5. (a) ϵ = 0.05 , (b) ϵ = 0.01.

The centers of final grids crossing the failure surface associated with ϵ = 0.01 are used to obtain failure probabilities. θ_s,th and θ_s,th are both set as 0.005 while r_th is set as 0.05. In this case the number of intersection points between a direction vector and the failure surface should be odd-numbered. For Method A, the failure probabilities associated with approximation approaches I and II are calculated as 5.490 × 10⁻² and 7.170 × 10⁻², respectively. The failure probabilities associated with Methods B, C, and D are calculated as 7.302 × 10⁻², 6.848 × 10⁻⁶, and 4.029 × 10⁻², respectively. For this problem, the results provided by Method A with approximation approach II and Method B are acceptable, while the results provided by Methods C and D are off the mark. It should be noted that although it is possible to obtain a satisfactory result using Method A, a relatively large N_r and a relatively small ϵ are entailed, which implies a large computational cost. Therefore, for the reliability problem associated with multiple disconnected regions, using direction sampling may not be the optimal option.

Example 6: Reliability analysis integrating clustering techniques

The integration of clustering techniques into the proposed approach is investigated herein. The example problem is associated with a nonlinear oscillator subjected to a rectangular pulse. An illustration of the oscillator and the external force is shown in Figure 15. This example problem was investigated in Zhang et al. (2021). The performance function is

G (c_{1}, c_{2}, m, r, t_{1}, F_{1}) = 3 r - | \frac{2 F}{m ω_{0}^{2}} \sin (\frac{ω_{0} t_{1}}{2}) |

(17)

with

ω_{0} = \sqrt{\frac{c_{1} + c_{2}}{m}}

. All the random variables involved herein have a normal distribution. The mean values of m, c₁, c₂, r, F₁, and t₁ are 1, 1, 0.1, 0.5, 0.6, and 1, respectively, while the standard deviations of these variables are 0.05, 0.1, 0.01, 0.05, 0.1, and 0.2, respectively. The failure probability was calculated in Zhang et al. (2021) using several algorithms. The failure probability is about 9.1 × 10⁻⁶.

Figure 15.

Illustration of a nonlinear oscillator subjected rectangular pulse (adapted from Zhang et al., (2021)). (a) Nonlinear oscillator. (b) Rectangular pulse.

Although the performance function only involves six random variables, simulating the entire 6-dimensional failure surface with a fine discretization entails large memory so that a typical desktop with a memory of 16 gigabytes cannot be used. Therefore, clustering technique needs to be adopted in order to solve the problem using an ordinary desktop. As the focus of this example is to check the accuracy of Method A using clustering techniques, only the failure probability associated with Method A is calculated. For a failure probability as small as 9.1 × 10⁻⁶, the number of samples used in Methods B and D is insufficient. When absolute value is involved in a performance function, Method C also cannot provide an accurate failure probability.

The lower bound U_l, upper bound U_b, and the initial length the of grid in each dimension U_int are set as −4, 4, and 2, respectively. The values of N_r, ϵ, and θ_s,th are 1, 0.15, and 0.2, respectively. Clustering process is conducted for the grids crossing the failure surface after the first division (i.e., ϵ₁=1). A total of 2¹¹ clusters are generated using K-means clustering algorithms. The number of grids within each cluster before further division is less than 10. For each generated direction samples, 12 clusters are selected based on the cross angles between the vector associated with the cluster centroid and the direction vector. To expedite the calculation process, it is assumed that if the minimum cross angle between the vector associated with the cluster centroid and the direction vector is larger than 0.3 radian, the direction vector does not intersect the failure surface and, therefore, the failure probability in that direction is zero. This is valid based on the cross angles between the vectors associated with the centroids of adjacent clusters. For the selected clusters, grid division and failure surface identification continue until U_int is less than ϵ.

A total of 2 × 10⁴ samples is used to make sure the convergence of failure probability when only a part of failure surface is finely discretized. The relationship between number of samples and the failure probability is shown in Figure 16. It can be observed that for this example problem, the failure probability stabilizes when the number of samples reaches 10⁴. The failure probability using 10⁴ samples in Method A is 9.83 × 10⁻⁶. When the number of samples further increases, some small-scale fluctuations exist and the peak of the fluctuation is 1.02 × 10⁻⁵. As only part of grid on the failure surface is divided in a refined manner for each direction sample, error between the calculated result and the accurate result can be observed, However, even for the maximum value of the fluctuation when the number of samples is larger than 10⁴, the relative error is within a reasonable range (about 12% in this case). Given the drastic reduction of the memory required and the acceleration of the computation speed, incentive exists on incorporating a clustering process in the entire lattice grid division process when the requirement on the accuracy of the reliability analysis is not very stringent.

Figure 16.

Relationship between number of samples and failure probability using Method A for Example 6.

Conclusions

A novel approach to conduct reliability analysis is investigated herein. This approach is based on densely populate points on the failure surface and uses direction sampling to obtain failure probability. An iterative space division process is conducted to obtain small lattice grids crossing the failure surface. For a generated direction vector sample, the intersection points between it and the failure surface are selected among these centers of grids, thereby avoiding a root-searching process. The novelty of the proposed approach is enhanced by extending the scope of application to problems associated with multiple intersections between the direction vector and the failure surface. The applicability of the proposed approach can be further enhanced by using clustering techniques.

The following conclusions are drawn:

1. When the intersection between the direction vector and the failure surface is unique, the approach proposed herein can be adopted to solve the associated reliability problem.

2. When a direction vector intersects the failure surface multiple times, the proposed approach yields satisfactory results by using approximation methods. Furthermore, the proposed approach is applicable to reliability problems with connected failure regions. A significant advantage of the proposed approach is its ability to calculate small failure probabilities.

3. For reliability problems associated with multiple disconnected regions, it is still possible for the proposed approach to yield satisfactory results. However, as the number of samples increases and the final grid size decreases, the associated computational cost is expected to increase significantly. This is a limitation of the proposed approach.

4. For reliability problems with a high dimensionality, using clustering techniques can significantly reduce the computational cost.

5. It should be acknowledged that the curse of dimensionality is still the most daunting challenge for the proposed approach to achieve wide application. In addition to clustering techniques, directional importance sampling and dimension reduction techniques (Breitung, 2015) can also be used to enhance the applicability of the proposed approach.

Footnotes

Acknowledgments

The authors are grateful for the support provided by the National Science Foundation Award CMMI-1537926. The opinions and conclusions presented this paper are those of the authors and do not necessarily reflect the views of the sponsoring organization.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Science Foundation Award CMMI-1537926.

ORCID iDs

Xu Han

Dan M. Frangopol

References

S-K

Beck

(2001) Estimation of small failure probabilities in high dimensions by subset simulation. Probabilistic Engineering Mechanics 16(4): 263–277. DOI: 10.1016/S0266-8920(01)00019-4

Bichon

Eldred

Swiler

, et al. (2008) Efficient global reliability analysis for nonlinear implicit performance functions. AIAA Journal 46(10): 2459–2468.

Bichon

Mcfarland

Mahadevan

(2011) Efficient surrogate models for reliability analysis of systems with multiple failure modes. Reliability Engineering & System Safety 96(10): 1386–1395. DOI: 10.1016/j.ress.2011.05.008

Breitung

(2015) Dimension reduction methods for reliability problems. In: 12th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP 2015, Vancouver, Canada, July 2015. The University of British Columbia. DOI: 10.14288/1.0076113

Breitung

(2017) Extrapolation, invariance, geometry and subset sampling. In: Proceedings of 14th International Probabilistic Workshop, Cham, Switzerland, December 2017, pp. 33–43. Springer International Publishing. DOI: 10.1007/978-3-319-47886-9

Bucher

Bourgund

(1990) A fast and efficient response surface approach for structural reliability problems. Structural Safety 7(1): 57–66. DOI: 10.1016/0167-4730(90)90012-E

Cheng

Papaioannou

, et al. (2023) Rare event estimation with sequential directional importance sampling. Structural Safety 100: 102291. DOI: 10.1016/j.strusafe.2022.102291

Ditlevsen

Bjerager

(1986) Methods of structural systems reliability. Structural Safety 3(3–4): 195–229. DOI: 10.1016/0167-4730(86)90004-4

Ditlevsen

Bjerager

(1989) Plastic reliability analysis by directional simulation. Journal of Engineering Mechanics 115(6): 1347–1362. DOI: 10.1061/(ASCE)0733-9399

10.

Faber

Gunzburger

(1999) Centroidal Voronoi tessellations: applications and algorithms. SIAM Review 41(4): 637–676. DOI: 10.1137/S0036144599352836

11.

Echard

Gayton

Lemaire

(2011) AK-MCS: an active learning reliability method combining Kriging and Monte Carlo Simulation. Structural Safety 33(2): 145–154. DOI: 10.1016/j.strusafe.2011.01.002

12.

Faravelli

(1989) Response-surface approach for reliability analysis. Journal of Engineering Mechanics 115(12): 2763–2781. DOI: 10.1061/(ASCE)0733-9399

13.

Fiessler

Rackwitz

Neumann

(1979) Quadratic limit states in structural reliability. Journal of the Engineering Mechanics Division 105(4): 661–676. DOI: 10.1061/JMCEA3.0002512

14.

Ghanem

Spanos

(1991) Spectral stochastic finite-element formulation for reliability analysis. Journal of Engineering Mechanics 117(10): 2351–2372. DOI: 10.1061/(ASCE)0733-9399

15.

Hartigan

Wong

(1979) Algorithm AS 136: a K-means clustering algorithm. Applied Statistics 28(1): 100–108. DOI: 10.2307/2346830

16.

Hasofer

Lind

(1974) Exact and invariant second-moment code format. Journal of the Engineering Mechanics division 100(1): 111–121. DOI: 10.1061/JMCEA3.0001848

17.

Hohenbichler

Rackwitz

(1982) First-order concepts in system reliability. Structural Safety 1(3): 177–188. DOI: 10.1016/0167-4730(82)90024-8

18.

Jafari-Asl

Ohadi

Ben Seghier

MEA

, et al. (2021) Accurate structural reliability analysis using an improved line-sampling-method-based slime mold algorithm. ASCE-ASME Journal of Risk and Uncertainty in Engineering Systems, Part A: Civil Engineering 7(2): 04021015. DOI: 10.1061/ajrua6.0001129

19.

Katsuki

Frangopol

(1994) Hyperspace division method for structural reliability. Journal of Engineering Mechanics 120(11): 2405–2427. DOI: 10.1061/(ASCE)0733-9399

20.

Koutsourelakis

Pradlwarter

Schuëller

(2004) Reliability of structures in high dimensions, part I: algorithms and applications. Probabilistic Engineering Mechanics 19(4): 409–417. DOI: 10.1016/j.probengmech.2004.05.001

21.

Lopez

Torii

Miguel

LFF

, et al. (2015) Overcoming the drawbacks of the FORM using a full characterization method. Structural Safety 54: 57–63. DOI: 10.1016/j.strusafe.2015.02.003

22.

Marsaglia

(1972) Choosing a point from the surface of a sphere. The Annals of Mathematical Statistics 43(2): 645–646. DOI: 10.1214/aoms/1177692644

23.

McKay

Beckman

Conover

(2000) A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 42(2): 55–61. DOI: 10.1080/00401706.2000.10485979

24.

Melchers

(1987) Structural Reliability: Analysis and Prediction. Hertfordshire, England: Ellis Horwood Ltd.

25.

Muller

(1959) A note on a method for generating points uniformly on N-dimensional spheres. Communications of the ACM 2(4): 19–20. DOI: 10.1145/377939.377946

26.

Nguyen

T-T

M-H

Nguyen

T-P

, et al. (2022) Seismic reliability analysis of building structures using subset simulation coupled with deep learning-based surrogate model. Advances in Structural Engineering 25(11): 2301–2318. DOI: 10.1177/13694332221092677

27.

Papaioannou

Straub

(2021) Combination line sampling for structural reliability analysis. Structural Safety 88: 102025. DOI: 10.1016/j.strusafe.2020.102025

28.

Pradlwarter

Schueller

Koutsourelakis

, et al. (2007) Application of line sampling simulation method to reliability benchmark problems. Structural Safety 29: 208–221. DOI: 10.1016/j.strusafe.2006.07.009

29.

Roussouly

Petitjean

Salaun

(2013) A new adaptive response surface method for reliability analysis. Probabilistic Engineering Mechanics 32: 103–115. DOI: 10.1016/j.probengmech.2012.10.001

30.

Schuëller

Pradlwarter

Koutsourelakis

(2004) A critical appraisal of reliability estimation procedures for high dimensions. Probabilistic Engineering Mechanics 19(4): 463–474. DOI: 10.1016/j.probengmech.2004.05.004

31.

Song

Kawai

(2023) Adaptive stratified sampling for structural reliability analysis. Structural Safety 101: 102292. DOI: 10.1016/j.strusafe.2022.102292

32.

Sun

Wang

, et al. (2017) LIF: a new Kriging based learning function and its application to structural reliability analysis. Reliability Engineering & System Safety 157: 152–165. DOI: 10.1016/j.ress.2016.09.003

33.

Tokdar

Kass

(2010) Importance sampling: a review. WIREs Computational Statistics 2(1): 54–60. DOI: 10.1002/wics.56

34.

Vahedi

Ghasemi

Miri

(2018) Structural reliability assessment using an enhanced adaptive Kriging method. Structural Engineering and Mechanics 66(6): 677–691. DOI: 10.12989/sem.2018.66.6.677

35.

Wang

, et al. (2022) AKSE: a novel adaptive Kriging method combining sampling region scheme and error-based stopping criterion for structural reliability analysis. Reliability Engineering & System Safety 219: 108214. DOI: 10.1016/j.ress.2021.108214

36.

Zhang

Wang

Sørensen

(2020) AKOIS: an adaptive Kriging oriented importance sampling method for structural system reliability analysis. Structural Safety 82: 101876. DOI: 10.1016/j.strusafe.2019.101876

37.

Zhang

Cheng

(2021) AK-DS: an adaptive Kriging-based directional sampling method for reliability analysis. Mechanical Systems and Signal Processing 156: 107610. DOI: 10.1016/j.ymssp.2021.107610

38.

Zhong

Wang

Dang

, et al. (2020) Structural reliability assessment by salp swarm algorithm–based FORM. Quality and Reliability Engineering International 36(4): 1224–1244. DOI: 10.1002/qre.2626

39.

Zuniga

Murangira

Perdrizet

(2021) Structural reliability assessment through surrogate based importance sampling with dimension reduction. Reliability Engineering & System Safety 207: 107289. DOI: 10.1016/j.ress.2020.107289