4.2.1 Dual Sampling Technique

If U represents a set of samples uniformly distributed in [0, 1], and the corresponding basic random variable is X(U), where X obeys the distribution of the probability density function (PDF) f(x₁, x₂, … x_n), there exist I-U and X(I-U) as well, which are negatively correlated with U and X(U). Therefore, the simulation estimation of Equation (4.4) is

(4.18)

Apparently, Equation (4.18) is the unbiased estimation of P_f, and the variance of simulation estimation can be expressed as

(4.19)

(4.20)

Correspondingly, the variance of the simulation estimation is

in is a variable expected to be estimated. So

(4.21)

(4.22)

Therefore, not only does the conditional expectation sampling technique reduce the variance of the sampling simulation, but it is also very beneficial to the calculation of truncated distribution probability in Equation (4.20).

4.2.3 Importance Sampling Technique

If there exists a sampling density function h(X), and the following relation is satisfied

(4.23)

then Equation (4.2) can be written in the form of importance sampling

(4.24)

Then the unbiased estimation of Equation (4.24) is

(4.25)

(4.26)

When the sampling density function is

(4.27)

the simulated variance of Equation (4.26) reaches a minimum. It should be said that Equation (4.27) just provides a way to select h(X), but h(X) is very difficult to select because this depends on the distribution form of the random variables, the LSF and the sampling simulation accuracy ^[4-3][4-4]. However, Equation (4.26) has both upper and lower bounds, which can be obtained by means of the Cauchy-Schwary inequality, i.e.,

(4.28)

Equation (4.28) can also be expressed as the ratio of the PDF f(X) of random variables to the sampling probability density function h(X), i.e.,

(4.29)

This expression provides the boundaries of the sampling function h(X). Obviously, Equation (4.27) is the upper boundary of the sampling function h(X), requiring that h(X) should satisfy the ratio in the failure zone D_f. Moreover, all samples should fall within D_f (see Equation (4.23)). At the same time, FU ^[4-5] has proved the upper boundary of Equation (4.29), as follows

(4.30)

where X^⋆ is the maximum likelihood point on the LSF. Considering that the ratio of Equation (4.26) is always greater than or equal to , and that f(X) can always find such a point in the gradient direction of X^⋆ to meet the conditions of Equation (4.27), this point is selected as the subdomain center of the failure zone so that at a given confidence level, the sampling mean obtained from this subdomain meets the conditions of Equation (4.27). This observation is very useful for the construction of sampling function h(X), implying that the sampling center may be in the gradient direction of X^⋆ and near X^⋆ in the failure zone.

The above observation has been partially reflected in the existing importance sampling methods ^{[4-6][4-7][4-8][4-9][4-10][4-11]}. However, there is still the tricky problem of how to effectively determine the importance sampling density function (type and parameters).

4.2.4 Stratified Sampling Method

The concept of stratified sampling ^[4-12][4-13] is analogous to that of importance sampling because both of them require the samples which contribute a lot to P_f. However, stratified sampling does not change the original density function, but simply divides the sampling interval into further subintervals, keeping different numbers of sampling points in the subintervals so that more samples should be taken from those subintervals which make greater contributions.

The domain of integration G(X)<0 is divided into M mutually disjointed subintervals L_j, and from each subinterval, N_j uniform random number vectors r uniformly distributed within this subinterval are taken. Here, N_j represents not only the number of uniform random number vectors generated in the j^-th sub-interval, but also the frequency of the uniform random number vectors falling into this subinterval from [0, 1]. In this way, the simulation results for stratified sampling can be written as:

(4.31)

(4.32)

(4.33)

In the simulated variance Equation (4.33), it can be proved that if

(4.33)

(4.34)

(4.35)

Suppose Equation (4.2) can be divided into two parts ^[4-2]:

(4.36)

(4.37)

(4.38)

The estimation formula for simulated variance is:

(4.39)

As can be seen from Equation (4.38), if y(X) is very close to f(X), i.e., , then . So, y(X) is called the control variate of f(X).

The significance of the control variates method is that if part of the known analytical solution is contained in the problem to be solved, then the variance of sampling simulation can be significantly reduced by removing this part and calculating the difference using the simulation method.

4.2.6 Correlated Sampling Method

Usually, one of the main purposes of simulation is to ascertain the effect of slight changes that take place in the system. Thus, it is necessary to carry out simulation repeatedly. If two tests are independent of each other, then the variance of the difference between the simulation results is the sum of the variances of each simulation; if the same random number is used in the two tests, the test results are highly correlated with each other, thus making it possible to reduce the test variance of the difference between the test results ^[4-2].

Consider the following two integrals:

(4.40)

(4.41)

(4.42)

(4.43)

(4.44)

(4.45)

(4.46)

(4.47)

(4.48)

(4.49)

so

(4.50)

(4.51)

where is the probability value of the random variable x_k at ; correspondingly, the PDF in Equations (4.52) and (4.53) is

(4.52)

or

(4.53)

If there exists an importance sampling density function h(X), as follows

(4.54)

Then the failure probability ^[4-4] in Equations (4.2) and (4.24) can be expressed as

(4.55)

or

(4.56)

Equations (4.55) and (4.56) are the expression (4.52) of conditional expectation sampling. According to the results of Equation (4.53), Equations (4.55) and (4.56) also have good sampling efficiency.

If a normal distribution probability density function is chosen as the sampling density function of x_i(i≠ k) in Equation (4.54), then Equations (4.55) and (4.56) will have higher sampling efficiency; during the selection of parameters for the sampling distribution probability density function , the sampling function h(X) is required to satisfy at least the first and second moments of the random variable in the failure zone, i.e.,

(4.57)

(4.58)

or

(4.59)

(4.60)

For the random variable x_i in X ∈ D_f, this can be selected according to Equation (4.49) or by the Taylor expansion of Z = G(X); here, x_i is allowed to take its value near the boundary domain of D_f, and the accumulated deviation will be eliminated in Equation (4.49); therefore, for the selection of x_k, it is necessary to consider the univariate form in LSF, as well as the variable with high discreteness in probability distribution.

After Equations (4.59) and (4.60) are determined, the importance sampling density function becomes . Therefore, Equation (4.55) can be rewritten as

(4.61)

where represents the sample generated by the importance sampling density function; φ[·] is a standard normally distributed probability density function.

4.3.2 Composite Important Sampling

For the extremal distribution function, can be obtained by the deterministic method. For other types of PDF, can be obtained by sampling. An effective importance sampling method will be given here.

If the importance sampling density function is set to a truncated distribution function ^[4-2], as follows

(4.62)

where d represents the truncated value, as shown in Equation (4.49); σ₁ is a parameter that varies with the probability distribution of the basic variable x. When the basic variable obeys the normal distribution N(μ, σ²), σ₁ can be expressed as

(4.63)

Figure 4.1 shows the form of the truncated distribution probability density function; Table 4.1 shows the change of σ₁ with σ when G(X)=3.0 – x; as can be seen, the results obtained in Equation (4.63) are obviously better than other results; when the basic variable obeys lognormal distribution, it can be transformed into normal distribution first, before the value of σ₁ can be worked out using Equation (4.63). It can then be inverted to lognormal distribution. Therefore, the sampling simulation of and can be expressed as

(4.64)

where is a sub-sample generated by h_T(x); M represents the number of samples. Table 4.2 shows a comparison of the results achieved by different sampling simulation methods. ISM is the important sampling method, whose sampling function obeys normal distribution; its mean is the “design checking point” and the mean-square error of sampling is set to unit value. IFM is an iteratively faster Monte-Carlo method; IISM is an improved numerical simulation method; AIISM is the IISM with dual sampling technique. As can be seen, IISM reflects the simulated sampling results of Equation (4.64).

Thus, the sampling simulation of Equation (4.61) can be expressed as

(4.65)

To sum up, in order to reduce the number of sampling simulations in the Monte-Carlo method and improve sampling efficiency, conditional expectation sampling can be combined with importance sampling to perform a composite sampling simulation to ensure the effectiveness of sampling every time. This also avoids the complexity caused by the use of an optimization method to determine sampling parameters. As a new way, sampling parameters can be determined by the moment method; the effectiveness of sampling simulation can be improved by adopting a normal distribution function for importance sampling and a truncated distribution function for conditional expectation variables; various examples confirm the effectiveness and wide applicability of this improved simulation method.

4.3.3 Calculation Steps

1. Select the basic variable x_k of conditional expectation from the LSF expression of the structure;
   
2. Determine parameters for the importance sampling density function according to Equations (4.57) and (4.58);
   
3. Generate simulation sample from , and calculate the truncated value d_i through Equation (4.49);
   
4. Determine sampling parameter σ₁ by Equation (4.63) according to the probability distribution of d_i and x_k;
   
5. Generate simulation sample from Equation (4.62), and calculate ;
   
6. Repeat step (5) M times, and calculate using Equation (4.64);
   
7. Calculate and ;
   
8. Repeating steps (3) to (7) N times;
   
9. Calculate the failure probability of the structure based on Equation (4.65).

4.4.1 V Space

If X space is an original basic random variable space, U space, a standard normally distributed variable space, can be obtained by means of the Rosenblatt transformation ^[4-14][4-15] x=T_xu (u). In U space, the joint PDF f(u) of random variables is monotonous, and u^⋆ is the maximum likelihood point on the failure surface, i.e., the minimum distance from the origin of coordinates, which can be obtained by FORM/SORM. This also means that the curvature at u^⋆ is always greater than or equal to −1/β ^[4-16]. P is the reliability index, and β=(u^⋆T u^⋆)^1/2.

To make full use of the quadratic effect on the failure surface, the failure function is expanded by Taylor series at u^⋆, with the quadratic term taken. So

(4.66)

where G_u(u^⋆) and G_uu (u^⋆) are first and second derivatives of G(u) at u^⋆. Equation (4.66) is orthogonally transformed, meaning

(4.67)

so that z_n is parallel to the direction cosines of -G_u(u^⋆)/|G_u(u^⋆)| at u^⋆. H can be obtained by the standard Gram-Schmidt Equation [4-17]. Thus, Equations (4.66) can be expressed as

(4.68)

(4.69)

κ is the eigenvector of matrix , i.e., the diagonal matrix of principal curvature. is the orthogonal eigenvector matrix of , while is a (n-1)×(n-1)-order matrix obtained from the matrix A with the n^th row and n^th column deleted; therefore, the new coordinate system for these variables corresponds to the main curvature coordinates. So

(4.70)

(4.71)

where T_{u, x} = P · H is the orthogonal transformation matrix for U space and V space. Therefore, V space is only the linear transformation of U space, and is also a standard normally distributed variable space ^[4-15].

By substituting Equation (4.71) into Equation (4.66), we can establish a failure function for V space:

(4.72)

where a_nn is an element on the n^th row and n^th column of the matrix A, and A_n is a vector on the n^th row or n^th column of the matrix A that does not contain a_nn. It should be said that Equation (4.72) is a general expression of the failure function with a quadratic effect in V space. For simplification, an approximate parabolic function is usually adopted to replace Equation (4.72). This type of approximate parabolic function is effective in describing the quadratic effect of the surface at u^⋆, as follows

(4.73)

(4.74)

Figure 4.2 shows the approximation of the parabolic surface in V space at u^⋆.

4.4.2 Importance Sampling Area

The importance sampling zone refers to the fact that the samples in this zone have an important contribution to the estimation of P_f, and it is composed of the sampling center and a region with a certain confidence level. In V space, such an importance sampling zone can be represented as an ellipse, as shown in Figure 4.3. According to the symmetry of the approximate parabolic surface and the observation of the importance sampling method, the sampling center is located on the v_n axis and within the failure region D_f, while the semi-axis length of the ellipse is related to the variance of sampling variables and changes within the geometric properties of the failure surface.

(1) Sampling center

As shown in Figure 4.3, the position of the sampling ellipse center changes with the principal curvature κ_i, while the sampling center is located on the v_n axis. Its distance from the maximum likelihood point v^⋆ on the failure surface is k_ib_i (i.e., in Figure 4.3), where b_i represents the semi-axis length of the sampling ellipse on the v_n axis, and parameter k_i should satisfy the following conditions:

• ① When κ_i= 0, k_i = k₀; this means that, when the failure surface is a plane, the sampling center is near v^⋆;
  
• ② When κ_i= -1/β, k_i = 0; this shows that, when the failure surface is a sphere in V space or a circle with the radius of β in the v_i – v_n plane, the sampling center is located on v^⋆;
  
• ③ When κ_i= +∞, k_i = 1. This is a corner case, i.e., the failure surface is a very convex surface, and the distance from the sampling center to the maximum likelihood point v^⋆ is equal to the semi-axis length of the sampling ellipse.

Therefore, the parameter k_i can be constructed into the following function:

(4.75)

where k₀ is related to the ratio δ₀ of the effective sample region A_eff to the whole sample region A_whole when κ_i=0. This can be obtained by the following equation:

(4.76)

Table 4.3 shows the numerical values of k₀ corresponding to different δ₀ values.

(2) Axial length of sampling ellipse

As can be seen from Figure 4.3, the semi-axis length of the sampling ellipse on the v_n axis can be expressed as

(4.77)

where Δβ is a parameter related to the simulation accuracy. Reference [4.20] provides a suitable range of Δβ, i.e., 0.7∼ 1.4. The length of the other half axis of the sampling ellipse can be obtained depending on the geometric properties of the failure surface:

• ① For the convex surface, i.e., when –β^-1 ≤ κ_i≤ 0, this is
  
  (4.78a)
  
  
• ② For the flat surface, i.e., when 0 ≤ κ_i ≤ κ_1,cr, this is
  
  (4.78b)
  
  
• ③ For the concave surface, i.e., when κ_i,cr < κ_i, this is
  
  (4.78c)
  
  

  where κ_i,cr can be obtained by equating Equation (4.78b) and (4.78c). Figure 4.4 shows the relationship between the principal curvature κ of the failure surface and the elliptic parameters a, b and k.

4.4.3 Importance Sampling Function

Choosing an appropriate sampling function is one of the main parts of importance sampling. Because V space is standard normally distributed variable space, an n-dimensional normal distribution probability density function can be selected as the importance sampling function h(v). The mean of the sampling function is

(4.79)

The standard variance of the sampling function involves the confidence level of the samples generated in the importance sampling zone. Let α be the confidence coefficient, and the double-axis length of the sampling ellipse be a confidence interval. The standard variance of sampling will be

(4.80)

where α_S = 1/Φ^-1 [(1+α)/2], and Φ^-1[·] is the reversal of the standard normal probability function. Figure 4.5 shows the effect of different confidence coefficients α on the simulation results. As can be seen, an appropriate α value ranges from 0.90 to 0.999, but it should be noted that α also depends on the geometric characteristics of the failure surface.

Therefore, the importance sampling expression of V space can be expressed as

(4.81)

where is the sample vector generated by the sampling function h(v).

4.4.4 Simulation Procedure

We can now make a comprehensive comment on the above process of constructing V space by the importance sampling method. As can be seen, the construction of h(v) is related not only to the geometric parameters (maximum likelihood point, gradient and curvature) of the failure surface, but also to parameters of the sampling function (effective sampling zone ratio δ₀, effective sampling zone Δβ, confidence coefficient α). The former directly affects the accuracy and efficiency of sampling, but it can be expressed by V-space sampling; the latter plays an important role in the process of constructing h(v), but it seems to be less sensitive to the simulation results given suitable parameters. When parameters are set to δ₀ = 0.8, Δβ = 1.0 and α = 0.95 in the following examples, ideal V-space sampling results can be achieved under all failure function conditions.

Based on the contents outlined in the previous sections, the simulation process of V-space importance sampling (ISM-V) can be described as follows:

1. Use FORM/SORM to determine the maximum likelihood point u^⋆ on the failure surface in U space;
   
2. Calculate the principal curvature κ at u^⋆ (Equation (4.69) and the transformation matrix for V space and U space (Equation (4.71);
   
3. Determine the center and axial length of the sampling ellipse in the v_i ∼ v_n plane (Equations (4.75), (4.77) and (4.78);
   
4. Calculate the mean and standard variance of the importance sampling function h (v) (Equations (4.79) and (4.80);
   
5. Generate a set of samples from h (v) (Equation (4.81);
   
6. Transform into and by T_uv and Rosenblatt transformation;
   
7. Check whether the sample is within the failure zone D_f, and work out the ratio;
   
8. Repeat step (5)∼(7) N times;
   
9. Calculate the failure probability using Equation (4.81).

4.4.5 Evaluation

A random variable space (U space) can be transformed into another random variable space (V space) through a linear orthogonal transformation matrix so as to reflect the quadratic effect of the failure surface at the maximum likelihood point more accurately. Thus, an importance sampling method (ISM) is established for V space. This type of ISM fully considers the geometric properties of the failure surface (e.g., maximum likelihood point, gradient, curvature, etc.), thus ensuring the validity of samples and improving calculation efficiency. With a wide variety of applications, it is not only used for the convex failure surface, but also for flat or concave failure surfaces; moreover, with increased simulation times, the simulation results show very stable convergence, close to the exact value.