https://doi.org/10.5194/ms-14-439-2023
https://doi.org/10.5194/ms-14-439-2023 # Modified control variates method based on second-order saddle-point approximation for practical reliability analysis

Xinong En, Yimin Zhang, and Xianzhen Huang
Abstract

A novel method is presented for efficiently analyzing the reliability of engineering components and systems with highly nonlinear complex limit state functions. The proposed method begins by transforming the integral of the limit state function into an integral of a highly correlated limit state function using the control variates method. The second-order reliability method is then employed within the control variates framework to approximate the highly correlated limit state function as a quadratic polynomial. Subsequently, the probability of failure is obtained through the estimation of the saddle-point approximation and a small number of samples generated by Latin hypercube sampling. To demonstrate the effectiveness of the proposed method, four examples involving mathematical functions and mechanical problems are solved. The results are compared with those obtained using the second-order reliability method (SORM), control variates based on Monte Carlo simulation (CVMCS) with second-order saddle-point approximation (SOSPA), importance sampling (IS) and Monte Carlo simulation (MCS). The findings demonstrate that, while maintaining high-precision reliability results, the proposed method significantly reduces the number of evaluations of the limit state function through a small number of initial samples. The method is capable of efficiently and accurately solving complex practical engineering reliability problems.

Share
Dates
1 Introduction

Safety evaluation of engineering components and systems is crucial for ensuring the economic and efficient operation of industrial equipment. Uncertainties in geometric and material properties of engineering structures arise due to environmental factors and operational influences. In light of these uncertainties, reliability analysis provides a valuable approach to quantifying the ability to fulfill specified system functions within specified conditions and time (Helton et al., 2006; Zhang et al., 2015; Zhao et al., 2017). The probability of failure (PF), based on the limit state function (LSF), is typically expressed as a multidimensional integral of the joint probability density function (PDF) over the random variables. However, directly calculating this multidimensional integral for engineering reliability is challenging. Consequently, various methods for reliability analysis have been developed in recent years. These methods can be broadly categorized into three categories: approximate analytical methods (Hasofer and Lind, 1974; Zhao and Ono, 1999; Meng et al., 2018a), simulation methods (Du, 2008a; Zhang et al., 2018; Meng et al., 2018b) and surrogate model methods (Zhang et al., 2017; Echard et al., 2011; Dai et al., 2015).

In practical engineering problems, the presence of a relatively large number of random variables can result in low computational efficiency for certain methods, such as Monte Carlo simulation (MCS). Therefore, it is preferable to use methods that can provide sufficiently accurate results with a minimum number of evaluations of the limit state function (Lemaire, 2013). In this regard, the approximate analytical methods, such as the first-order reliability method (FORM) and the second-order reliability method (SORM), are advantageous as they can approximate the limit state function with relatively low computational costs. The FORM is widely recognized as one of the most mainstream approximate analytical reliability methods. In recent years, significant improvements have been developed, as evidenced by studies conducted by Chiralaksanakul and Mahadevan (2005), Shayanfar et al. (2017) and Keshtegar and Chakraborty (2018). However, in highly nonlinear complex problems, SORM offers a more accurate approximation of the LSF by using a quadratic polynomial surface. Recent studies have focused on improving SORM, with notable contributions from researchers such as Breitung (1984), Tvedt (1990), Der Kiureghian et al. (1987) and Zhang and Du (2010). Although they have advantages in terms of computation, their accuracy decreases when dealing with reliability problems that involve non-normal random variables and nonlinear limit state functions. Therefore, applying these methods may not be reasonable in such cases. To improve accuracy, some methods based on FORM and SORM have been developed.

Among these methods, saddle-point approximation (SPA) stands out with its advantages of simple calculation, operational feasibility and estimation of the PF without the need to approximate the LSF. Due to the efficiency and accuracy of SPA, scholars have further explored its application in reliability engineering (Xiao et al., 2012; Du, 2008b; Guo, 2014) and integrated it into approximate analytical methods (Du and Sudjianto, 2004; Huang and Du, 2008; Huang et al., 2018). Du and Sudjianto (2004) proposed the first-order saddle-point approximation (FOSA) method, which involves expanding the cumulant-generating function of the LSF using a first-order Taylor series expansion around the most probable point (MPP) and estimating the PF using saddle-point approximation. Huang and Du (2008) introduced the modified variance function-based saddle-point approximation method (MVFOSA), which transforms highly nonlinear LSFs into approximate linear functions at the mean point of each random variable, resulting in a more efficient approach. Huang et al. (2018) developed a new second-order reliability method based on saddle-point approximation, which obtains a second-order approximation of the LSF by performing a second-order Taylor series expansion around the MPP. However, it is important to acknowledge the existence of approximation errors and the potential accuracy disparity between approximate analytical methods and MCS in engineering applications.

Since approximate analytical methods do not introduce any errors into probability estimation, they cannot quantify the inherent approximation errors of these methods. Therefore, there exists a difference in accuracy between approximate analytical methods and simulation methods. MCS (Metropolis and Ulam, 1949; Zio, 2013) and its variants, such as importance sampling (IS) (Papaioannou et al., 2019) and subset simulation (SS) (Au and Beck, 2001), are more robust for complex limit state functions. MCS is most accurate as a standard for testing other theoretical methods. Despite these advantages, in order to achieve an acceptable level of accuracy, the number of random samples required by MCS grows exponentially for practical engineering problems with small failure probabilities. If the failure probability is on the order of 10k, MCS requires 10k+210k+3 samples to achieve the corresponding accuracy (Binder, 1987). Because the variance of MCS estimates is inversely proportional to the failure probability, some studies (Melchers, 1989; Au and Beck, 2001; Rubinstein and Kroese, 2016; Oates et al., 2016) have explored the use of variance reduction techniques in simulation methods to reduce the variance of MCS.

Among these variance reduction techniques, the control variable (CV) method (Law et al., 2000) stands out for its relatively simple operation and broad applicability, as it has almost no restrictions on its range of application. The CV method has been employed recently in reliability analysis. For instance, Kawai (2020) utilized both importance sampling and control variates to construct and investigate an adaptive variance reduction framework. Portier and Segers (2019) successfully reduced the variance of Monte Carlo integration using control variates in a multiple linear regression model. Dimarco and Pareschi (2020) introduced a multiscale control variates strategy to enhance the efficiency of numerical calculations for kinetic equations with stochastic parameters, particularly in high-dimensional problems. Additionally, in engineering practice, where reliability problems often involve high dimensionality and small failure probabilities, the CV method is frequently combined with other techniques. Rashki (2018) introduced a new formulation for approximating small failure probabilities, leveraging the adaptive cooperation of SS and the CV method. To minimize the number of evaluations of the LSF, Ameryan et al. (2022) proposed a novel approach that integrates an active learning kriging metamodel (AK-MCS) and the sequential space conversion method (SESC), where the SESC formulation is derived using the control variates method. Moreover, to enhance the efficiency of addressing high-dimensional problems, Mehni and Mehni (2023) presented a method for improving the performance of cross-entropy-based Gaussian mixture importance sampling, incorporating a control variates scheme to further enhance its effectiveness.

In this study, the control variates method is combined with variance reduction sampling techniques to enhance the efficiency of probability estimation. These techniques, with significantly smaller sample sizes compared to MCS-based CVs, often result in a substantial reduction in the number of evaluations of the LSF. Among these techniques, Latin hypercube sampling (LHS), initially proposed by Mckay et al. (1979), is a widely used multidimensional stratified sampling method in engineering practice. Zou et al. (2008) employed an optimal symmetric LHS technique to construct the response surface of the indicator function for multiple limit state and multiple design point problems. Yang et al. (2005) utilized the optimal LHS method to uniformly distribute sampling points across the entire design space, thereby improving the accuracy of a real-world frontal impact metamodel. Christensen et al. (2017) employed a novel constrained nested orthogonal-maximin Latin hypercube approach for sampling to determine variable-shape mold configurations used in the training and test sets. LHS exhibits characteristics such as quick convergence and uniformly distributed samples. Therefore, using LHS to generate samples for the control variates method can significantly improve computational efficiency.

To address the reliability issues of engineering components and systems characterized by high nonlinearity and complexity with limited samples, this study proposes a method that combines second-order saddle-point approximation (SOSPA) with control variable Latin hypercube sampling (CVLHS). The proposed method leverages the benefits of both variance reduction techniques and approximate analytical methods. By doing so, it achieves high-precision reliability results while significantly reducing the number of evaluations of the limit state function. The remainder of the paper is organized as follows. Section 2 discusses the method combining CVLHS and SOSPA. Section 3 outlines the implementation procedure of the proposed method in real-world applications. Section 4 presents two numerical examples and two engineering cases to demonstrate the efficiency and accuracy of the proposed method. Finally, Sect. 5 provides conclusions.

2 CVLHS with SOSPA

To increase the efficiency of the MCS, based on the theory of the CV method, the index function of the highly correlated limit state function (LSF) as the CV is used to convert the original index function in the PF multidimensional integral. By incorporating LHS into the CV method, only a small number of samples need to be generated to remove the errors raised from the highly correlated LSF. Combining the expectation of the CV estimated by SOSPA, PF is estimated accurately without any simplification and assumptions.

## 2.1 CV for the failure probability function

The PF is expressed as the integral of the joint PDF fX(x) over the n-dimensional random variable vector $\mathbit{x}=\left[{x}_{\mathrm{1}},{x}_{\mathrm{2}},\mathrm{\dots },{x}_{n}\right]$:

$\begin{array}{}\text{(1)}& {P}_{f}=\underset{F}{\int }{f}_{X}\left(\mathbit{x}\right)\mathrm{d}\mathbit{x}={\int }_{X}{I}_{F}\left(\mathbit{x}\right){f}_{X}\left(\mathbit{x}\right)\mathrm{d}x,\end{array}$

where F is the failure domain, $F=\left\{\mathbit{x}:g\left(\mathbit{x}\right)\le \mathrm{0}\right\}$, g(x) is the LSF and IF is the index function:

This integral can be described as an expectation of IF(x), and with the N randomly sampled variable vector x, it is approximated as follows:

$\begin{array}{}\text{(3)}& {P}_{f}={\int }_{X}{I}_{F}\left(\mathbit{x}\right){f}_{X}\left(\mathbit{x}\right)\mathrm{d}\mathbit{x}=E\left[{I}_{F}\left(\mathbit{x}\right)\right]\approx \frac{\mathrm{1}}{N}\sum _{i=\mathrm{1}}^{N}{I}_{F}\left({\mathbit{x}}_{i}\right).\end{array}$

If gc(x) is the highly correlated LSF of g(x), the corresponding failure domain can be denoted as ${F}_{\mathrm{c}}=\left\{\mathbit{x}:{g}_{\mathrm{c}}\left(\mathbit{x}\right)\le \mathrm{0}\right\}$. Moreover, ${I}_{{F}_{\mathrm{c}}}$ is the index function as a CV which is the high correlation of IF, and similarly to Eq. (2), . According to the CV-reduction technique (Rubinstein and Kroese, 2016), the new index function is defined by

$\begin{array}{}\text{(4)}& {I}_{{F}_{\mathrm{CV}}}\left(\mathbit{x}\right)={I}_{F}\left(\mathbit{x}\right)-\mathit{\rho }\left[{I}_{{F}_{\mathrm{c}}}\left(\mathbit{x}\right)-E\left[{I}_{{F}_{\mathrm{c}}}\left(\mathbit{x}\right)\right]\right],\end{array}$

where ρ is the correction coefficient.

Therefore, by substituting Eq. (4) into Eq. (2), the failure probability integral by CV is rewritten as

$\begin{array}{}\text{(5)}& \begin{array}{rl}{P}_{{f}_{\mathrm{CV}}}& ={\int }_{X}{I}_{{F}_{\mathrm{CV}}}\left(\mathbit{x}\right){f}_{X}\left(\mathbit{x}\right)\mathrm{d}x={\int }_{X}{I}_{F}\left(\mathbit{x}\right){f}_{X}\left(\mathbit{x}\right)\mathrm{d}x\\ & -\mathit{\rho }{\int }_{X}{I}_{{F}_{\mathrm{c}}}\left(\mathbit{x}\right){f}_{X}\left(\mathbit{x}\right)\mathrm{d}x+\mathit{\rho }E\left[{I}_{{F}_{\mathrm{c}}}\left(\mathbit{x}\right)\right]\\ & \approx \mathit{\rho }\cdot {P}_{{f}_{\mathrm{c}}}+\frac{\mathrm{1}}{N}\sum _{i=\mathrm{1}}^{N}\left({I}_{F}\left({\mathbit{x}}_{i}\right)-\mathit{\rho }\cdot {I}_{{F}_{\mathrm{c}}}\left({\mathbit{x}}_{i}\right)\right),\end{array}\end{array}$

where ${P}_{{f}_{\mathrm{c}}}$ is the PF of the highly correlated LSF gc(x). To simplify the failure probability function, the posterior portion of Eq. (5) will be zero, i.e., $\frac{\mathrm{1}}{N}\sum _{i=\mathrm{1}}^{N}\left({I}_{F}\left({\mathbit{x}}_{i}\right)-\mathit{\rho }\cdot {I}_{{F}_{\mathrm{c}}}\left({\mathbit{x}}_{i}\right)\right)=\mathrm{0}$. In the light of Eq. (3), the correction coefficient ρ is estimated by the ratio of the expectation of the index functions:

Then, Eq. (5) can be simplified to

$\begin{array}{}\text{(7)}& {P}_{{f}_{\mathrm{CV}}}=\mathit{\rho }\cdot {P}_{{f}_{\mathrm{c}}}.\end{array}$

## 2.2 Quadratic polynomial surface approximation of the LSF

In this study, to ensure the accuracy of reliability analysis for the high-nonlinearity reliability problems, the quadratic polynomial is selected as the highly correlated LSF gc(x). By calculating the second derivative, the quadratic polynomial approximation of g(x) can consider the concave direction and curvature of the limit state surface near the MPP. The highly correlated LSF is expanded into a quadratic term of Taylor series at the MPP x as follows:

$\begin{array}{}\text{(8)}& \begin{array}{rl}{g}_{\mathrm{c}}\left(\mathbit{x}\right)& =g\left({\mathbit{x}}^{\ast }\right)+\mathrm{\nabla }g\left({\mathbit{x}}^{\ast }{\right)}^{\mathrm{T}}\left(\mathbit{x}-{\mathbit{x}}^{\ast }\right)\\ & +\frac{\mathrm{1}}{\mathrm{2}}\left(\mathbit{x}-{\mathbit{x}}^{\ast }{\right)}^{\mathrm{T}}{\mathrm{\nabla }}^{\mathrm{2}}g\left({\mathbit{x}}^{\ast }\right)\left(\mathbit{x}-{\mathbit{x}}^{\ast }\right),\end{array}\end{array}$

where g(x) is the vector of partial derivatives at the MPP; 2g(x) is the Hessian matrix of g(x) at the MPP; and x is the MPP, which is solved by iterative optimization using the mean of x as the initial search point:

where ||x|| is the distance from the origin to the x on the limit state surface.

## 2.3 Latin hypercube sampling

The correction coefficient ρ in Eq. (6) is obtained based on MCS generally. Nevertheless, for complicated high-dimensional problems, MCS requires huge amounts of computation for accurate estimation. To further reduce the computational cost, LHS is employed to generate the initial samples based on the original LSF g(x) and the highly correlated LSF gc(x) denoted by the Taylor expanded quadratic polynomial, respectively.

To sample N initial points of x, the independent U(0,l) (uniformly distributed on [0,1]) random variables are represented as U(i,j) ($i=\mathrm{1},\mathrm{\dots },N,j=\mathrm{1},\mathrm{\dots },n\right)$. The corresponding scaling probability for each random number U(i,j) can be given by

$\begin{array}{}\text{(10)}& P\left(i,j\right)=\frac{U\left(i,j\right)+\left(j-\mathrm{1}\right)}{N}\phantom{\rule{0.25em}{0ex}}\left(i=\mathrm{1},\mathrm{\dots },n,j=\mathrm{1},\mathrm{\dots },N\right).\end{array}$

Then the corresponding components ${x}_{i,j}\left(i=\mathrm{1},\mathrm{\dots },N,j=\mathrm{1},\mathrm{\dots },n\right)$ can be generated by inverse distribution functions. LHS considers the total sampling range of random variables with a small number of samples. Under the premise of ensuring accuracy, the efficiency of the sampling is improved. The initial correction coefficient ρ0 is obtained by the N initial samples x, and the result of ρ is optimized by Nadd by adding samples until more than 20 failure points are reached. The process of LHS is shown in Fig. 1. Figure 1The procedure to generate the initial samples and optimize the correction coefficient by LHS.

## 2.4 SOSPA for the PF of the highly correlated LSF

For the computational advantage of SPA on a quadratic polynomial, SOSPA is more accurate than the other common SORM (Huang et al., 2018). Therefore, it is utilized to approximate the PF of the highly correlated LSF in this study. It employs an inverse Fourier transform to obtain the cumulative distribution function (CDF) of an exponential power series expression based on a saddle point by the moment-generating function (MGF) and the cumulant-generating function (CGF).

The MGF (Broda and Paolella, 2011) of the highly correlated LSF gc(x), where $\mathbit{x}\sim N\left(\mathit{\mu },\mathit{\sigma }\right)$, is given by

$\begin{array}{}\text{(11)}& \begin{array}{rl}& {M}_{{g}_{\mathrm{c}}}\left(t\right)=E\left[{e}^{t{g}_{\mathrm{c}}}\right]={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{t{g}_{\mathrm{c}}}{f}_{X}\left(\mathbit{x}\right)\mathrm{d}\mathbit{x}\\ & =\mathrm{exp}\left\{t\left(a+{\mathbit{b}}^{\mathrm{T}}\mathbit{\mu }+{\mathbit{\mu }}^{\mathrm{T}}\mathbit{c}\mathbit{\mu }\right)+\frac{{t}^{\mathrm{2}}}{\mathrm{2}}\sum _{i=\mathrm{1}}^{n}\frac{{d}_{i}^{\mathrm{2}}}{\mathrm{1}-\mathrm{2}t{\mathit{\lambda }}_{i}}\right\}\\ & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\prod _{i=\mathrm{1}}^{n}{\left(\mathrm{1}-\mathrm{2}t{\mathit{\lambda }}_{i}\right)}^{-\mathrm{0.5}},\end{array}\end{array}$

where ${\mathit{\lambda }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_{n}$ and P are the eigenvalues and eigenvectors of c and $\mathbit{P}\mathbit{c}{\mathbit{P}}^{\mathrm{T}}=\text{diag}\left({\mathit{\lambda }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_{n}\right)$ with P orthogonal, $\mathbit{d}={\left[{d}_{\mathrm{1}},{d}_{\mathrm{2}},\mathrm{\dots },{d}_{n}\right]}^{\mathrm{T}}={\mathbit{P}}^{\mathrm{T}}\left({\mathit{\sigma }}^{\mathrm{0.5}}\mathbit{b}+\mathrm{2}{\mathit{\sigma }}^{\mathrm{0.5}}\mathbit{c}\mathbit{\mu }\right)$, $a=\mathrm{0.5}{\mathbit{x}}^{\ast T}{\mathrm{\nabla }}^{\mathrm{2}}{g}_{\mathrm{c}}\left({\mathbit{x}}^{\ast }\right){\mathbit{x}}^{\ast }-\mathrm{\nabla }{g}_{\mathrm{c}}\left({\mathbit{x}}^{\ast }{\right)}^{\mathrm{T}}{\mathbit{x}}^{\ast }$, $\mathbit{b}=\mathrm{\nabla }{g}_{\mathrm{c}}\left({\mathbit{x}}^{\ast }\right)-{\mathrm{\nabla }}^{\mathrm{2}}{g}_{\mathrm{c}}\left({\mathbit{x}}^{\ast }\right){\mathbit{x}}^{\ast }$ and $\mathbit{c}=\mathrm{0.5}{\mathrm{\nabla }}^{\mathrm{2}}{g}_{\mathrm{c}}\left({\mathbit{x}}^{\ast }\right)$.

Moreover, the CGF of gc is given by

$\begin{array}{}\text{(12)}& \begin{array}{rl}{K}_{{g}_{\mathrm{c}}}\left(t\right)& =\mathrm{ln}\left({M}_{{g}_{\mathrm{c}}}\left(t\right)\right)=t\left(a+{\mathbit{b}}^{\mathrm{T}}\mathbit{\mu }+{\mathbit{\mu }}^{\mathrm{T}}\mathbit{c}\mathbit{\mu }\right)\\ & +\frac{{t}^{\mathrm{2}}}{\mathrm{2}}\sum _{i=\mathrm{1}}^{n}{d}_{i}^{\mathrm{2}}{\mathit{\vartheta }}_{i}+\frac{\mathrm{1}}{\mathrm{2}}\sum _{i=\mathrm{1}}^{n}\mathrm{ln}{\mathit{\vartheta }}_{i},\end{array}\end{array}$

where ${\mathit{\vartheta }}_{i}=\left(\mathrm{1}-\mathrm{2}t{\mathit{\lambda }}_{i}{\right)}^{-\mathrm{1}}$.

Then the derivatives of ${K}_{{g}_{\mathrm{c}}}\left(t\right)$ can be determined as follows.

$\begin{array}{}\text{(13)}& \begin{array}{rl}{K}_{{g}_{\mathrm{c}}}^{\left(\mathrm{1}\right)}\left(t\right)& =a+{\mathbit{b}}^{\mathrm{T}}\mathbit{\mu }+{\mathbit{\mu }}^{\mathrm{T}}\mathbit{c}\mathbit{\mu }\\ & +\sum _{i=\mathrm{1}}^{n}\left(t{d}_{i}^{\mathrm{2}}{\mathit{\vartheta }}_{i}+{t}^{\mathrm{2}}{d}_{i}^{\mathrm{2}}{\mathit{\vartheta }}_{i}^{\mathrm{2}}{\mathit{\lambda }}_{i}+{\mathit{\vartheta }}_{i}{\mathit{\lambda }}_{i}\right)\end{array}\end{array}$

$\begin{array}{}\text{(14)}& {K}_{{g}_{\mathrm{c}}}^{\left(\mathrm{2}\right)}\left(t\right)=\sum _{i=\mathrm{1}}^{n}\left({d}_{i}^{\mathrm{2}}{\mathit{\vartheta }}_{i}+\mathrm{4}t{d}_{i}^{\mathrm{2}}{\mathit{\vartheta }}_{i}^{\mathrm{2}}{\mathit{\lambda }}_{i}+\mathrm{4}{t}^{\mathrm{2}}{d}_{i}^{\mathrm{2}}{\mathit{\vartheta }}_{i}^{\mathrm{3}}{\mathit{\lambda }}_{i}^{\mathrm{2}}+\mathrm{2}{\mathit{\vartheta }}_{i}^{\mathrm{2}}{\mathit{\lambda }}_{i}^{\mathrm{2}}\right)\end{array}$

$\begin{array}{}\text{(15)}& \begin{array}{rl}& {K}_{{g}_{\mathrm{c}}}^{\left(\mathrm{3}\right)}\left(t\right)\\ & =\sum _{i=\mathrm{1}}^{n}\left(\mathrm{6}{d}_{i}^{\mathrm{2}}{\mathit{\vartheta }}_{i}^{\mathrm{2}}{\mathit{\lambda }}_{i}+\mathrm{24}t{d}_{i}^{\mathrm{2}}{\mathit{\vartheta }}_{i}^{\mathrm{3}}{\mathit{\lambda }}_{i}^{\mathrm{2}}+\mathrm{24}{t}^{\mathrm{2}}{d}_{i}^{\mathrm{2}}{\mathit{\vartheta }}_{i}^{\mathrm{4}}{\mathit{\lambda }}_{i}^{\mathrm{3}}+\mathrm{8}{\mathit{\vartheta }}_{i}^{\mathrm{3}}{\mathit{\lambda }}_{i}^{\mathrm{3}}\right)\end{array}\end{array}$

$\begin{array}{}\text{(16)}& \begin{array}{rl}{K}_{{g}_{\mathrm{c}}}^{\left(\mathrm{4}\right)}\left(t\right)=\sum _{i=\mathrm{1}}^{n}\left(& \mathrm{48}{d}_{i}^{\mathrm{2}}{\mathit{\vartheta }}_{i}^{\mathrm{3}}{\mathit{\lambda }}_{i}^{\mathrm{2}}+\mathrm{192}t{d}_{i}^{\mathrm{2}}{\mathit{\vartheta }}_{i}^{\mathrm{4}}{\mathit{\lambda }}_{i}^{\mathrm{3}}\\ & +\mathrm{192}{t}^{\mathrm{2}}{d}_{i}^{\mathrm{2}}{\mathit{\vartheta }}_{i}^{\mathrm{5}}{\mathit{\lambda }}_{i}^{\mathrm{4}}+\mathrm{48}{\mathit{\vartheta }}_{i}^{\mathrm{4}}{\mathit{\lambda }}_{i}^{\mathrm{4}}\right)\end{array}\end{array}$

By numerically integrating the PDF of SPA (Lugannani and Rice, 1980), the CDF of gc can be approximated as

$\begin{array}{}\text{(17)}& {F}_{{g}_{\mathrm{c}}}\left({g}_{\mathrm{c}\mathrm{0}}\right)=P\left\{{g}_{\mathrm{c}}\le {g}_{\mathrm{c}\mathrm{0}}\right\}\approx \mathrm{\Phi }\left(w\right)+\mathit{\varphi }\left(w\right)\left(\frac{\mathrm{1}}{w}-\frac{\mathrm{1}}{v}\right),\end{array}$

where Φ(⋅) and φ(⋅) are the CDF and PDF of the standard normal distribution, respectively, and w and v are parameters expressed as

$\begin{array}{}\text{(18)}& w=\text{sgn}\left({t}_{\mathrm{s}}\right){\left\{\mathrm{2}\left[{t}_{\mathrm{s}}{g}_{\mathrm{c}\mathrm{0}}-{K}_{{g}_{\mathrm{c}}}\left({t}_{\mathrm{s}}\right)\right]\right\}}^{\mathrm{0.5}},\end{array}$

$\begin{array}{}\text{(19)}& v={t}_{\mathrm{s}}{\left[{K}_{{g}_{\mathrm{c}}}^{\left(\mathrm{2}\right)}\left({t}_{\mathrm{s}}\right)\right]}^{\mathrm{0.5}},\end{array}$

where ts is the saddle point, ${K}_{{g}_{\mathrm{c}}}^{\left(\mathrm{2}\right)}\left(\cdot \right)$ are the second-order derivatives of ${K}_{{g}_{\mathrm{c}}}\left(\cdot \right)$, and sgn(ts) is a symbolic function and can be +1, 1 or 0, corresponding to ts>0, ts<0 and ts=0.

Furthermore, the second-order CDF of gc is

$\begin{array}{}\text{(20)}& \begin{array}{rl}{\stackrel{\mathrm{̃}}{F}}_{{g}_{\mathrm{c}}}\left({g}_{\mathrm{c}\mathrm{0}}\right)& =P\left\{{g}_{\mathrm{c}}\le {g}_{\mathrm{c}\mathrm{0}}\right\}\\ & \approx {F}_{{g}_{\mathrm{c}}}\left({g}_{\mathrm{c}\mathrm{0}}\right)-\mathit{\varphi }\left(w\right)\\ & \left({v}^{-\mathrm{1}}\left(\frac{{\mathit{\kappa }}_{\mathrm{4}}}{\mathrm{8}}-\frac{\mathrm{5}}{\mathrm{24}}{\mathit{\kappa }}_{\mathrm{3}}^{\mathrm{2}}\right)-{v}^{-\mathrm{3}}-\frac{{\mathit{\kappa }}_{\mathrm{3}}}{\mathrm{2}{u}^{\mathrm{2}}}+{w}^{-\mathrm{3}}\right),\end{array}\end{array}$

where .

Then, making gc0=0, the PF ${P}_{{f}_{\mathrm{c}}}$ can be easily calculated from the second-order CDF ${\stackrel{\mathrm{̃}}{F}}_{{g}_{\mathrm{c}}}\left(\mathrm{0}\right)$ corresponding to the strongly correlated LSF gc(x)≤0, i.e.,

$\begin{array}{}\text{(21)}& {P}_{{f}_{\mathrm{c}}}={\stackrel{\mathrm{̃}}{F}}_{{g}_{\mathrm{c}}}\left(\mathrm{0}\right).\end{array}$
3 Implementation procedure

To ensure the accuracy of the PF, the coefficient of variation of the PF is calculated as the stop condition for loop operation. The new sample vector xnew is sequentially added to the sample set until the coefficient of variation of the PF has to be convergent.

## 3.1 The estimator statistical property

The expectation of PF in Eq. (7) is estimated as

$\begin{array}{}\text{(22)}& E\left({P}_{{f}_{\mathrm{CV}}}\right)=E\left(\mathit{\rho }\cdot {P}_{{f}_{\mathrm{c}}}\right)=\mathit{\rho }\cdot E\left({\stackrel{\mathrm{̃}}{F}}_{{g}_{\mathrm{c}}}\left(\mathrm{0}\right)\right).\end{array}$

Since the estimation PF by LHS is the unbiased estimation (McKay et al., 1979), the expectations of the PF of the original LSF and the highly correlated LSF are

$\begin{array}{}\text{(23)}& E\left({P}_{{f}_{\mathrm{LHS}}}\right)=E\left({I}_{F}\left(\mathbit{x}\right)\right)\approx \frac{\mathrm{1}}{N}\sum _{i=\mathrm{1}}^{N}{I}_{F}\left({\mathbit{x}}_{i}\right),\end{array}$

$\begin{array}{}\text{(24)}& E\left({P}_{{f}_{{c}_{\mathrm{LHS}}}}\right)=E\left({I}_{{F}_{\mathrm{c}}}\left(\mathbit{x}\right)\right)\approx \frac{\mathrm{1}}{N}\sum _{i=\mathrm{1}}^{N}{I}_{{F}_{\mathrm{c}}}\left({\mathbit{x}}_{i}\right).\end{array}$

Therefore, the correction coefficient ρ is estimated by LHS:

The expectation of PF is given by

$\begin{array}{}\text{(26)}& E\left({P}_{{f}_{\mathrm{CV}}}\right)=\frac{E\left({P}_{{f}_{\mathrm{LHS}}}\right)}{E\left({P}_{{f}_{{\mathrm{c}}_{\mathrm{LHS}}}}\right)}\cdot E\left({\stackrel{\mathrm{̃}}{F}}_{{g}_{\mathrm{c}}}\left(\mathrm{0}\right)\right).\end{array}$

According to the theory of control variate estimator property, the variance of PF by the proposed approach may be estimated as

$\begin{array}{}\text{(27)}& \begin{array}{rl}\text{Var}\left({P}_{{f}_{\mathrm{CV}}}\right)& =\text{Var}\left({P}_{{f}_{\mathrm{LHS}}}\right)+{\mathit{\rho }}^{\mathrm{2}}\text{Var}\left({P}_{{f}_{{c}_{\mathrm{LHS}}}}\right)\\ & -\mathrm{2}\mathit{\rho }\text{Cov}\left({P}_{{f}_{\mathrm{LHS}}},{P}_{{f}_{{c}_{\mathrm{LHS}}}}\right)\\ & =\text{Var}\left({P}_{{f}_{\mathrm{LHS}}}\right)+\frac{{E}^{\mathrm{2}}\left({P}_{{f}_{\mathrm{LHS}}}\right)}{{E}^{\mathrm{2}}\left({P}_{{f}_{{c}_{\mathrm{LHS}}}}\right)}\cdot \text{Var}\left({P}_{{f}_{{c}_{\mathrm{LHS}}}}\right)\\ & -\mathrm{2}\cdot \frac{E\left({P}_{{f}_{\mathrm{LHS}}}\right)}{E\left({P}_{{f}_{{c}_{\mathrm{LHS}}}}\right)}\cdot \text{Cov}\left({P}_{{f}_{\mathrm{LHS}}},{P}_{{f}_{{c}_{\mathrm{LHS}}}}\right).\end{array}\end{array}$

Due to the expectation and the variance of ${P}_{{f}_{\mathrm{CV}}}$, the coefficient of variation ${\mathit{\delta }}_{{P}_{{f}_{\mathrm{CV}}}}$ can be estimated as

$\begin{array}{}\text{(28)}& \begin{array}{rl}& {\mathit{\delta }}_{{P}_{{f}_{\mathrm{CV}}}}=\frac{\sqrt{\text{Var}\left({P}_{{f}_{\mathrm{CV}}}\right)}}{E\left({P}_{{f}_{\mathrm{CV}}}\right)}\\ & =\frac{\sqrt{\text{Var}\left({P}_{{f}_{\mathrm{LHS}}}\right)+\frac{{E}^{\mathrm{2}}\left({P}_{{f}_{\mathrm{LHS}}}\right)}{{E}^{\mathrm{2}}\left({P}_{{f}_{{c}_{\mathrm{LHS}}}}\right)}\cdot \text{Var}\left({P}_{{f}_{{c}_{\mathrm{LHS}}}}\right)-\mathrm{2}\cdot \frac{E\left({P}_{{f}_{\mathrm{LHS}}}\right)}{E\left({P}_{{f}_{{c}_{\mathrm{LHS}}}}\right)}\cdot \text{Cov}\left({P}_{{f}_{\mathrm{LHS}}},{P}_{{f}_{{c}_{\mathrm{LHS}}}}\right)}}{\frac{E\left({P}_{{f}_{\mathrm{LHS}}}\right)}{E\left({P}_{{f}_{{c}_{\mathrm{LHS}}}}\right)}\cdot E\left({\stackrel{\mathrm{̃}}{F}}_{{g}_{\mathrm{c}}}\left(\mathrm{0}\right)\right)}\\ & =\sqrt{\frac{{E}^{\mathrm{2}}\left({P}_{{f}_{{c}_{\mathrm{LHS}}}}\right)\cdot \text{Var}\left({P}_{{f}_{\mathrm{LHS}}}\right)}{{E}^{\mathrm{2}}\left({P}_{{f}_{\mathrm{LHS}}}\right)\cdot {\stackrel{\mathrm{̃}}{F}}_{{g}_{\mathrm{c}}}^{\mathrm{2}}\left(\mathrm{0}\right)}+\frac{\text{Var}\left({P}_{{f}_{{c}_{\mathrm{LHS}}}}\right)}{{\stackrel{\mathrm{̃}}{F}}_{{g}_{\mathrm{c}}}^{\mathrm{2}}\left(\mathrm{0}\right)}-\frac{\mathrm{2}\cdot E\left({P}_{{f}_{{c}_{\mathrm{LHS}}}}\right)\cdot \text{Cov}\left({P}_{{f}_{\mathrm{LHS}}},{P}_{{f}_{{c}_{\mathrm{LHS}}}}\right)}{E\left({P}_{{f}_{\mathrm{LHS}}}\right)\cdot {\stackrel{\mathrm{̃}}{F}}_{{g}_{\mathrm{c}}}^{\mathrm{2}}\left(\mathrm{0}\right)}}.\end{array}\end{array}$

Hereinto, the variances of the PF by LHS are expressed as

Stein (1987) proved that $\frac{\mathrm{1}}{N}\left(N-\mathrm{1}\right)\text{cov}\left({I}_{{F}_{\mathrm{c}}}\left({\mathbit{x}}_{\mathrm{1},j}\right),{I}_{{F}_{\mathrm{c}}}\left({\mathbit{x}}_{\mathrm{2},j}\right)\right)$ is a non-positive value. Therefore, LHS has a smaller variance than MCS and is easier to converge.

## 3.2 Procedure of CVLHS with SOSPA

To solve the reliability problem, the procedure to perform the CVLHS with the SOSPA method is shown in Fig. 2. According to the convergence criteria, the threshold value of the coefficient of variation ${\mathit{\delta }}_{{P}_{{f}_{\mathrm{CV}}}}$ is set to 0.15–0.25 for different cases in this paper.

4 Illustrative examples

In this section, four representative examples are provided to demonstrate the accuracy and efficiency of the proposed method in addressing the characteristics of high nonlinearity and complexity in practical engineering problems. These examples include a highly nonlinear numerical problem, a small failure probability problem, and two continuous-system engineering cases. The results obtained from the proposed method, i.e., CVLHS with SOSPA, are compared with those acquired from other methods, including SORM, CVMCS with SOSPA, IS, and MCS. For the IS, the importance sampling density function kX(x) is selected to have the same form as the PDF fX(x), with equal means and n times the standard deviation, i.e., μk=μf and σk=nσf (where n is set to 2). In four cases, the initial sample size for both CVLHS with SOSPA and CVMCS with the SOSPA methods is set to 102, while, for the MCS method, the initial sample size is 106. To eliminate computational errors, the average of five calculations is taken.

## 4.1 Highly nonlinear numerical problem

To examine the applicability of the proposed method for the highly nonlinear reliability problem, the LSF with several independent normal random variables is analyzed as follows:

$\begin{array}{}\text{(31)}& G\left(\mathbit{x}\right)=\sum _{i=\mathrm{1}}^{\mathrm{10}}{x}_{i}^{i}-{\mathrm{10}}^{\mathrm{8}}.\end{array}$

The means and standard deviations of the random variables Xi ($i=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{10}$) are all μi=10, σi=5.

The results of reliability analysis are listed in Table 1. Figure 3 presents the variation curve of the PF by the proposed CVLHS with the SOSPA method with the sample size. Figure 4 compares the convergence procedure of the proposed method and CVMCS with the SOSPA method, and Fig. 5 shows the CDF curves of the LSF by the proposed method and MCS.

Table 1Comparison of the results of Example 4.1. * LSF call = LSF call(optimize MPP) + LSF call(sampling). Figure 4The curve of the variation coefficient with sample size in Example 4.1.

From Table 1, it can be observed that the number of LSF calls is the lowest in SORM among these methods. However, there is a relatively large relative error in the reliability results compared to MCS. IS yields satisfactory results, but it requires a relatively high number of LSF calls. On the other hand, the proposed CVLHS with SOSPA, compared to the CVMCS with SOSPA, further reduces the number of LSF calls while providing more accurate reliability results. In Figs. 3 and 4, both the proposed CVLHS with the SOSPA method and CVMCS with the SOSPA method demonstrate convergence with a small number of samples, with the proposed method converging more quickly and being more efficient. Furthermore, as depicted in Fig. 5, the results obtained from the proposed method closely align with those obtained from MCS. However, the number of LSF calls for the proposed method is merely 0.85 % of that required by MCS. Hence, the proposed CVLHS with the SOSPA method effectively and accurately solves highly nonlinear reliability problems.

## 4.2 Small failure probability problem

This example illustrates the efficiency and accuracy of the proposed method for the small failure reliability analysis with a high-dimensional LSF. In this example, the LSF, consisting of 40 independent and normally distributed random variables xi ($i=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{40}\right)\sim N\left(\mathrm{1.5},\mathrm{1}\right)$, is defined as (Bichon et al., 2008)

$\begin{array}{}\text{(32)}& G\left(\mathbit{x}\right)=\frac{\left({x}_{\mathrm{1}}^{\mathrm{2}}+\mathrm{4}\right)\left({x}_{\mathrm{2}}-\mathrm{1}\right)}{\mathrm{20}}-\mathrm{cos}\left(\frac{\mathrm{5}{x}_{\mathrm{1}}}{\mathrm{2}}\right)+\sum _{i=\mathrm{1}}^{\mathrm{40}}{x}_{i}^{\mathrm{2}}-\mathrm{71.5}.\end{array}$

For the small failure probability problem, the results in Table 2 show that the proposed CVLHS with SOSPA, CVMCS with SOSPA, and IS are more accurate than the SORM method. However, in order to achieve the corresponding accuracy, the number of calls to the LSF in IS exceeds 105. Figure 6 demonstrates that the proposed CVLHS with the SOSPA method converges when the sample size reaches 10 000. Additionally, Fig. 7 illustrates the convergence curves of the proposed method and the CVMCS with the SOSPA method, indicating that the proposed method converges at a faster rate. Consequently, the proposed method significantly reduces computational consumption for high-dimensional small failure probability reliability problems. Furthermore, Fig. 8 exhibits strong agreement between the results of the proposed CVLHS with the SOSPA method and those of MCS. However, the number of calls to the LSF for the proposed method is only 1.46 % of that required by MCS. Therefore, it can be verified that the proposed method in this paper significantly reduces the number of calls to the LSF while maintaining accuracy, making it suitable for reliability estimation of problems with small failure probabilities.

Table 2Comparison of the results of Example 4.2. * LSF call = LSF call(optimize MPP) + LSF call(sampling). Figure 7The curve of the variation coefficient with a sample size in Example 4.2.

## 4.3 Dynamic reliability analysis of a cantilever beam

This example is to verify the accuracy of the proposed CVLHS with the SOSPA method in an engineering structure. The cantilever beam under a dynamic loading F(t) (Bucher, 2009) is shown in Fig. 9. The dynamic motion equation of the cantilever beam is

$\begin{array}{}\text{(33)}& \mathit{\rho }A\frac{{\partial }^{\mathrm{2}}\mathit{\omega }}{\partial {t}^{\mathrm{2}}}+\text{EI}\frac{{\partial }^{\mathrm{4}}\mathit{\omega }}{\partial {x}^{\mathrm{4}}}=p\left(x,t\right),\end{array}$

where ρ is the constant density, A is the cross-sectional area, and EI is the bending stiffness.

According to the appropriate initial and boundary conditions, the fundamental natural circular frequency ω1 is given by

$\begin{array}{}\text{(34)}& {\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}=\frac{{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{4}}\text{EI}}{\mathit{\rho }A}=\frac{\mathrm{12.362}\text{EI}}{\mathit{\rho }A{L}^{\mathrm{4}}}.\end{array}$

The density and the length are deterministic; i.e., ρ=1 and L=1. The distribution characteristics are shown in Table 3. To avoid resonance failure of the cantilever beam, the ratio is $\mathit{\omega }/{\mathit{\omega }}_{\mathrm{1}}<\mathrm{0.99}$ or $\mathit{\omega }/{\mathit{\omega }}_{\mathrm{1}}>\mathrm{1.01}$.

For an engineering structure, as shown in Table 4 and Figs. 10–12, the reliability results of the proposed CVLHS with the SOSPA method fit well with the MCS results, and the relative error is minimal. When the sample size reaches around 1000, the method starts to converge. Compared to the SORM method, the number of calls to the LSF only increased by less than 1300, but the accuracy significantly improved. Compared to the IS method, in order to achieve the corresponding accuracy, the IS requires more than 21 times the number of LSF calls compared to CVLHS with SOSPA.

Table 4Comparison of the reliability results for the cantilever beam dynamic system. * LSF call = LSF call(optimize MPP) + LSF call(sampling). Figure 10The PF curve by CVLHS with the SOSPA method for the cantilever beam dynamic system. Figure 11The curve of the variation coefficient for the cantilever beam dynamic system.

## 4.4 Reliability of chatter stability in turning

A simplified diagram of a single-degree-of-freedom regenerative chatter dynamics model of the turning process system is shown in Fig. 13.

According to the dynamics analysis of the turning process system (Schmitz and Smith, 2009), the limit cutting width is expressed as

$\begin{array}{}\text{(35)}& \begin{array}{rl}& {a}_{\mathrm{plim}}=\frac{k\left[{\left(\mathrm{1}-{\mathit{\lambda }}^{\mathrm{2}}\right)}^{\mathrm{2}}+{\left(\mathrm{2}\mathit{\xi }\mathit{\lambda }\right)}^{\mathrm{2}}\right]}{{K}_{\mathrm{c}}u\left({\mathit{\lambda }}^{\mathrm{2}}-\mathrm{1}\right)}\\ & \cdot \frac{\mathrm{1}-\mathit{\eta }\mathrm{cos}\left[\mathrm{2}j\mathit{\pi }+\mathrm{arcsin}\frac{\mathrm{2}\mathit{\xi }\mathit{\lambda }}{\mathit{\eta }\sqrt{{\left(\mathrm{2}\mathit{\xi }\mathit{\lambda }\right)}^{\mathrm{2}}+{\left(\mathrm{1}-{\mathit{\lambda }}^{\mathrm{2}}\right)}^{\mathrm{2}}}}-\mathrm{arctan}\left(\frac{\mathrm{2}\mathit{\xi }\mathit{\lambda }}{\mathrm{1}-{\mathit{\lambda }}^{\mathrm{2}}}\right)\right]}{\mathrm{1}+{\mathit{\eta }}^{\mathrm{2}}-\mathrm{2}\mathit{\eta }\mathrm{cos}\left[\mathrm{2}j\mathit{\pi }+\mathrm{arcsin}\frac{\mathrm{2}\mathit{\xi }\mathit{\lambda }}{\mathit{\eta }\sqrt{{\left(\mathrm{2}\mathit{\xi }\mathit{\lambda }\right)}^{\mathrm{2}}+{\left(\mathrm{1}-{\mathit{\lambda }}^{\mathrm{2}}\right)}^{\mathrm{2}}}}-\mathrm{arctan}\left(\frac{\mathrm{2}\mathit{\xi }\mathit{\lambda }}{\mathrm{1}-{\mathit{\lambda }}^{\mathrm{2}}}\right)\right]},\end{array}\end{array}$

and the corresponding spindle speed can be expressed as

$\begin{array}{}\text{(36)}& \mathrm{\Omega }=\frac{\mathrm{60}\mathit{\lambda }{\mathit{\omega }}_{n}}{\mathrm{2}j\mathit{\pi }+\mathrm{arcsin}\frac{\mathrm{2}\mathit{\xi }\mathit{\lambda }}{\mathit{\eta }\sqrt{{\left(\mathrm{2}\mathit{\xi }\mathit{\lambda }\right)}^{\mathrm{2}}+{\left(\mathrm{1}-{\mathit{\lambda }}^{\mathrm{2}}\right)}^{\mathrm{2}}}}-\mathrm{arctan}\left(\frac{\mathrm{2}\mathit{\xi }\mathit{\lambda }}{\mathrm{1}-{\mathit{\lambda }}^{\mathrm{2}}}\right)},\end{array}$

where ωn is the natural frequency of the vibration system in turning, ${\mathit{\omega }}_{n}^{\mathrm{2}}=k/m$, ξ is the damping ratio of the turning system defined as $\mathit{\xi }=c/\mathrm{2}m{\mathit{\omega }}_{n}$, u is the direction coefficient defined as $u=\mathrm{cos}\left(\mathit{\phi }-\mathit{\alpha }\right)\mathrm{cos}\mathit{\alpha }$, η is the turning processing overlap coefficient $\left(\mathrm{0}\le \mathit{\eta }\le \mathrm{1}\right)$, $\mathit{\lambda }=\mathit{\omega }/{\mathit{\omega }}_{n}$ (vibration frequency ω is slightly larger than the natural frequency ωn), i.e., λ>1, and $j=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots }$ is the number of lobes. The cutting force coefficient Kc consists of the tangential cutting force coefficient Kt and the radial cutting force coefficient Kr, and its expression is

$\begin{array}{}\text{(37)}& {K}_{\mathrm{c}}=\sqrt{{K}_{\mathrm{t}}^{\mathrm{2}}+{K}_{\mathrm{r}}^{\mathrm{2}}}.\end{array}$

The minimum aplim appears when chatter is strongest, i.e., η=1, and the first partial derivative of λ in Eq. (35) is taken to obtain the minimum aplim as

$\begin{array}{}\text{(38)}& {a}_{{\mathrm{plim}}_{min}}=\frac{\mathrm{2}k\mathit{\xi }\left(\mathrm{1}+\mathit{\xi }\right)}{{K}_{\mathrm{c}}u}.\end{array}$

The turning process system would be unstable when the limit cutting width aplim is not more than the minimum limit cutting width. The distribution characteristics of the turning process parameters are listed in Table 5.

In Fig. 14, the reliability curve of the turning process at different spindle speeds generated by the proposed CVLHS with the SOSPA method is compared with that by SORM, CVMCS with SOSPA, IS and MCS. The results demonstrate that the proposed method achieves the highest consistency with the results obtained by MCS. Therefore, the method has been proven to be feasible for accurately addressing complex practical mechanical engineering reliability problems.

5 Conclusion

To improve the efficiency and accuracy of reliability analysis for engineering systems with high nonlinearity and complexity in practical applications, a novel method is proposed. The method, CVLHS with SOSPA, utilizes CV to transform the original LSF into a highly correlated LSF that can be approximated using SORM. Sample points are generated using LHS, and the probability of failure for the highly correlated LSF is estimated using SPA. The proposed method aims to reduce the number of function calls required by simulation method that employ variance reduction technique, by integrating approximate analytical method.

Through the reliability analysis of the examples, it is evident that the proposed method achieves accurate approximation of the failure probability with a reasonable number of samples, whereas MCS and IS require a large number of samples to achieve the desired level of accuracy. Additionally, when the relative error of the reliability calculation results obtained by SORM is relatively high, the proposed method still provides accurate estimation with a reasonable sample size. The results demonstrate that, for complex nonlinear problems, the proposed method requires less than 1.5 % of the LSF calls compared to MCS while maintaining a high level of accuracy. Moreover, the comparison between CVLHS with SOSPA and CVMCS with SOSPA results reveals that, even in situations where a large number of samples cannot be utilized in the modeling process, which is often a concern in practical engineering, CVLHS with SOSPA can accurately predict the failure probability using a small number of initial samples generated by LHS, thus further reducing the required number of LSF calls. Hence, the proposed method significantly reduces the number of LSF calls while ensuring high-precision reliability results.

While the proposed method in this study combines the advantages of an approximate analytical method and a simulation method, it is suitable for engineering problems where an explicit limit state function can be established. To further investigate and expand this method, the integration of surrogate model methods can be considered, and the reliability analysis results can be compared to identify the optimal integration approach.

Data availability

The data generated during this study are available from the corresponding author on reasonable request.

Author contributions

XE: Resources, Methodology, Validation, Formal analysis, Writing-Original Draft, Writing-Review & Editing. YZ: Project administration, Funding acquisition, Supervision. XH: Conceptualization, Investigation, Project administration, Funding acquisition, Supervision.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

We would like to express our appreciation to the R&D Program of the Beijing Municipal Education Commission, the National Key R&D Program of China, the National Natural Science Foundation of China, as well as the Beijing Municipal Education Commission and the Beijing Natural Science Foundation Co-financing Project for supporting this research.

Financial support

This research has been supported by the R&D Program of the Beijing Municipal Education Commission (grant no. KM202310015004), the National Key R&D Program of China (grant no. 2019YFB2004400), the National Natural Science Foundation of China (grant nos. 51975110 and U22B2087), as well as the Beijing Municipal Education Commission and the Beijing Natural Science Foundation Co-financing Project (grant no. KZ202210015019) for supporting this research.

Review statement

This paper was edited by Daniel Condurache and reviewed by three anonymous referees.

References

Ameryan, A., Ghalehnovi, M., and Rashki, M.: AK-SESC: a novel reliability procedure based on the integration of active learning kriging and sequential space conversion method, Reliab. Eng. Syst. Safe, 217, 108036, https://doi.org/10.1016/j.ress.2021.108036, 2022.

Au, S. K. and Beck, J. L.: Estimation of Small Failure Probabilities in High Dimensions by Subset Simulation, Probabilistic Eng. Mech., 16, 263–277, https://doi.org/10.1016/s0266-8920(01)00019-4, 2001.

Bichon, B. J., Eldred, M. S., Swiler, L. P., Mahadevan, S., and McFarland, J. M.: Efficient Global Reliability Analysis for Nonlinear Implicit Performance Functions, AIAA J., 46, 2459–2468, https://doi.org/10.2514/1.34321, 2008.

Binder, K. (Ed.): Applications of the Monte Carlo Method in Statistical Physics, in: Topics in Current Physics, Springer Berlin Heidelberg, ISBN: 978-3-540-17650-3, 1987.

Breitung, K.: Asymptotic Approximations for Multinormal Integrals, J. Eng. Mech., 110, 357–366, https://doi.org/10.1061/(asce)0733-9399(1984)110:3(357), 1984.

Broda, S. A. and Paolella, M. S.: Saddlepoint Approximations: A Review and Some New Applications, in: Springer Handbooks of Computational Statistics, edited by: Gentle, J. E., Härdle, W. K., and Mori, Y., Springer, Berlin, Heidelberg, 953–983, https://doi.org/10.1007/978-3-642-21551-3_32, 2011.

Bucher, C.: Computational Analysis of Randomness in Structural Mechanics, CRC Press, Leiden, ISBN 978-0-415-40354-2, 2009.

Chiralaksanakul, A. and Mahadevan, S.: First-Order Approximation Methods in Reliability-Based Design Optimization, J. Mech. Design, 127, 851–857, https://doi.org/10.1115/1.1899691, 2005.

Christensen, E. T., Lund, E., and Lindgaard, E.: Experimental Validation of Surrogate Models for Predicting the Draping of Physical Interpolating Surfaces, J. Mech. Design, 140, 011401, https://doi.org/10.1115/1.4038073, 2017.

Dai, H., Zhang, B., and Wang, W.: A Multiwavelet Support Vector Regression Method for Efficient Reliability Assessment, Reliab. Eng. Syst. Safe, 136, 132–139, https://doi.org/10.1016/j.ress.2014.12.002, 2015.

Der Kiureghian, A., Lin, H. Z., and Hwang, S. J.: Second-Order Reliability Approximations, J. Eng. Mech., 113, 1208–1225, https://doi.org/10.1061/(asce)0733-9399(1987)113:8(1208), 1987.

Dimarco, G. and Pareschi, L.: Multiscale Variance Reduction Methods Based on Multiple Control Variates for Kinetic Equations with Uncertainties, Multiscale Model. Sim., 18, 351–382, https://doi.org/10.1137/18m1231985, 2020.

Du, X.: Unified Uncertainty Analysis by the First Order Reliability Method, J. Mech. Design, 130, 091401, https://doi.org/10.1115/1.2943295, 2008a.

Du, X.: Saddlepoint Approximation for Sequential Optimization and Reliability Analysis, J. Mech. Design, 130, 011011, https://doi.org/10.1115/1.2717225, 2008b.

Du, X. and Sudjianto, A.: First Order Saddlepoint Approximation for Reliability Analysis, AIAA J., 42, 1199–1207, https://doi.org/10.2514/1.3877, 2004.

Echard, B., Gayton, N., and Lemaire, M.: AK-MCS: An Active Learning Reliability Method Combining Kriging and Monte Carlo Simulation, Struct. Saf., 33, 145–154, https://doi.org/10.1016/j.strusafe.2011.01.002, 2011.

Guo, S.: An Efficient Third-Moment Saddlepoint Approximation for Probabilistic Uncertainty Analysis and Reliability Evaluation of Structures, Appl. Math. Model., 38, 221–232, https://doi.org/10.1016/j.apm.2013.06.026, 2014.

Hasofer, A. M. and Lind, N. C.: Exact and Invariant Second-Moment Code Format, J. Eng. Mech. Div.-ASCE, 100, 111–121, 1974.

Helton, J. C., Johnson, J. D., Sallaberry, C. J., and Storlie, C. B.: Survey of Sampling-Based Methods for Uncertainty and Sensitivity Analysis, Reliab. Eng. Syst. Safe, 91, 1175–1209, https://doi.org/10.1016/j.ress.2005.11.017, 2006.

Huang, B. and Du, X.: Probabilistic Uncertainty Analysis by Mean-Value First Order Saddlepoint Approximation, Reliab. Eng. Syst. Safe, 93, 325–336, https://doi.org/10.1016/j.ress.2006.10.021, 2008.

Huang, X., Li, Y., Zhang, Y., and Zhang, X.: A New Direct Second-Order Reliability Analysis Method, Appl. Math. Model., 55, 68–80, https://doi.org/10.1016/j.apm.2017.10.026, 2018.

Kawai, R.: Adaptive Importance Sampling and Control Variates, J. Math. Anal. Appl., 483, 123608, https://doi.org/10.1016/j.jmaa.2019.123608, 2020.

Keshtegar, B. and Chakraborty, S.: A Hybrid Self-Adaptive Conjugate First Order Reliability Method for Robust Structural Reliability Analysis, Appl. Math. Model., 53, 319–332, https://doi.org/10.1016/j.apm.2017.09.017, 2018.

Law, A. M., Kelton, W. D., and Kelton, W. D.: Simulation Modeling and Analysis, McGraw-Hill, New York, 2000.

Lemaire, M.: Structural reliability, John Wiley & Sons, ISBN 978-1-84821-082-0, 2013.

Lugannani, R. and Rice, S.: Saddle Point Approximation for The Distribution of the Sum of Independent Random Variables, Adv. Appl. Probab., 12, 475–490, https://doi.org/10.1017/s0001867800050278, 1980.

McKay, M. D., Beckman, R. J., and Conover, W. J.: A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code, Technometrics, 21, 239–245, https://doi.org/10.2307/1268522, 1979.

Mehni, M. B. and Mehni, M. B.: Reliability analysis with cross-entropy based adaptive Markov chain importance sampling and control variates, Reliab. Eng. Syst. Safe, 231, 109014, https://doi.org/10.1016/j.ress.2022.109014, 2023.

Melchers, R. E.: Importance Sampling in Structural Systems, Struct. Saf., 6, 3–10, https://doi.org/10.1016/0167-4730(89)90003-9, 1989.

Meng, Z., Zhou, H., Hu, H., and Keshtegar, B.: Enhanced Sequential Approximate Programming Using Second Order Reliability Method for Accurate and Efficient Structural Reliability-Based Design Optimization, Appl. Math. Model., 62, 562–579, https://doi.org/10.1016/j.apm.2018.06.018, 2018a.

Meng, Z., Zhang, D., Liu, Z., and Li, G.: An Adaptive Directional Boundary Sampling Method for Efficient Reliability-Based Design Optimization, J. Mech. Design., 140, 121406, https://doi.org/10.1115/1.4040883, 2018b.

Metropolis, N. and Ulam, S.: The Monte Carlo Method, J. Am. Stat. Assoc., 44, 335–341, https://doi.org/10.1080/01621459.1949.10483310, 1949.

Oates, C. J., Girolami, M., and Chopin, N.: Control Functionals for Monte Carlo Integration, J. R. Stat. Soc. B, 79, 695–718, https://doi.org/10.1111/rssb.12185, 2016.

Papaioannou, I., Geyer, S., and Straub, D.: Improved cross entropy-based importance sampling with a flexible mixture model, Reliab. Eng. Syst. Safe., 191, 106564, https://doi.org/10.1016/j.ress.2019.106564, 2019.

Portier, F. and Segers, J.: Monte Carlo Integration with a Growing Number of Control Variates, J. Appl. Probab., 56, 1168–1186, https://doi.org/10.1017/jpr.2019.78, 2019.

Rashki, M.: Hybrid control variates-based simulation method for structural reliability analysis of some problems with low failure probability, Appl. Math. Model., 60, 220–234, 2018.

Rubinstein, R. Y. and Kroese, D. P.: Simulation and the Monte Carlo method, John Wiley & Sons, New York, ISBN 9781118632161, 2016.

Schmitz, L. and Smith, K.: Machining Dynamics: Frequency Response to Improved Productivity, Springer, New York, ISBN 978-0-387-09644-5, 2009.

Shayanfar, M. A., Barkhordari, M. A., and Roudak, M. A.: An Efficient Reliability Algorithm for Locating Design Point Using the Combination of Importance Sampling Concepts and Response Surface Method, Commun. Nonlinear Sci., 47, 223–237, https://doi.org/10.1016/j.cnsns.2016.11.021, 2017.

Stein, M.: Large Sample Properties of Simulations Using Latin Hypercube Sampling, Technometrics, 29, 143–151, https://doi.org/10.1080/00401706.1987.10488205, 1987.

Tvedt, L.: Distribution of Quadratic Forms in Normal Space–Application to Structural Reliability, J. Eng. Mech., 116, 1183–1197, https://doi.org/10.1061/(asce)0733-9399(1990)116:6(1183), 1990.

Xiao, N. C., Huang, H. Z., Wang, Z., Liu, Y., and Zhang, X. L.: Unified Uncertainty Analysis by The Mean Value First Order Saddlepoint Approximation, Struct. Multidiscip. O., 46, 803–812, https://doi.org/10.1007/s00158-012-0794-4, 2012.

Yang, R. J., Wang, N., Tho, C. H., Bobineau, J. P., and Wang, B. P.: Metamodeling Development for Vehicle Frontal Impact Simulation, J. Mech. Design, 127, 1014–1020, https://doi.org/10.1115/1.1906264, 2005.

Zhang, D., Han, X., Jiang, C., Liu, J., and Li, Q.: Time-Dependent Reliability Analysis Through Response Surface Method, J. Mech. Design, 139, 041404, https://doi.org/10.1115/1.4035860, 2017.

Zhang, J. and Du, X.: A Second-Order Reliability Method With First-Order Efficiency, J. Mech. Design, 132, 101006, https://doi.org/10.1115/1.4002459, 2010.

Zhang, J., Xiao, M., Gao, L., and Chu, S.: A Combined Projection-Outline-Based Active Learning Kriging and Adaptive Importance Sampling Method for Hybrid Reliability Analysis with Small Failure Probabilities, Comput. Method. Appl. M., 344, 13–33, https://doi.org/10.1016/j.cma.2018.10.003, 2018.

Zhang, L., Lu, Z., and Wang, P.: Efficient Structural Reliability Analysis Method Based on Advanced Kriging Model, Appl. Math. Model., 39, 781–793, https://doi.org/10.1016/j.apm.2014.07.008, 2015.

Zhao, H., Li, S., and Ru, Z.: Adaptive Reliability Analysis Based on a Support Vector Machine and Its Application to Rock Engineering, Appl. Math. Model., 44, 508–522, https://doi.org/10.1016/j.apm.2017.02.020, 2017.

Zhao, Y. G. and Ono, T.: A General Procedure for First/Second-Order Reliability Method (FORM/SORM), Struct. Saf., 21, 95–112, https://doi.org/10.1016/S0167-4730(99)00008-9, 1999.

Zio, E.: The Monte Carlo Simulation Method for System Reliability and Risk Analysis, Springer, London, ISBN 978-1-4471-4587-5, 2013.

Zou, T., Mourelatos, Z. P., Mahadevan, S., and Tu, J.: An Indicator Response Surface Method for Simulation-Based Reliability Analysis, J. Mech. Design, 130, 071401, https://doi.org/10.1115/1.2918901, 2008.