Second-order reliability methods: a review and comparative study

Second-order reliability methods are commonly used for the computation of reliability, defined as the probability of satisfying an intended function in the presence of uncertainties. These methods can achieve highly accurate reliability predictions owing to a second-order approximation of the limit-state function around the Most Probable Point of failure. Although numerous formulations have been developed, the lack of full-scale comparative studies has led to a dubiety regarding the selection of a suitable method for a specific reliability analysis problem. In this study, the performance of commonly used second-order reliability methods is assessed based on the problem scale, curvatures at the Most Probable Point of failure, first-order reliability index, and limit-state contour. The assessment is based on three performance metrics: capability, accuracy, and robustness. The capability is a measure of the ability of a method to compute feasible probabilities, i.e., probabilities between 0 and 1. The accuracy and robustness are quantified based on the mean and standard deviation of relative errors with respect to exact reliabilities, respectively. This study not only provides a review of classical and novel second-order reliability methods, but also gives an insight on the selection of an appropriate reliability method for a given engineering application.


Introduction
Reliability analysis is widely used in industrial applications, ranging from structural (Chojaczyk et al. 2015;Mansour et al. 2019a;Hultgren et al. 2021a) and mechanical engineering (Papadimitriou et al. 2020;Sandberg et al. 2017) to materials design (Liu et al. 2018;Mansour et al. 2019b;Alzweighi et al. 2020). Regardless of the application, the reliability is defined as the probability that an intended function is satisfied in the presence of uncertainties. It is typically predicted during the design stage in order to avoid the catastrophic consequences of failure with substantial certainty. The fundamental task of reliability analysis is to compute the reliability defined by the probability integral (Choi et al. 2007).
or the probability of failure given by X, and g(X) is a limit-state function. The limit state is defined by g(X) = 0 and the failure event is indicated by g(X) < 0. Generally, the probability of failure is computed for an implicitly given g(X), since the explicit analytical expression is unknown. This is typical for numerous engineering applications where each evaluation of g(X) involves a black-box simulation such as Finite-Element Analysis (FEA) (Gu et al. 2001;Hultgren et al. 2021b). Direct assessment of the multidimensional integral is generally difficult or even impossible to perform. One method to evaluate the integration is Monte-Carlo (MC) simulation (Wang et al. 2019). However, when the reliability is very high (approaching 1.0), the computational effort of MC simulation is prohibitively expensive (Du and Chen 2000). Therefore, Hasofer and Lind (Hasofer and Lind 1974) proposed the concept of the Most Probable Point of failure (MPP) to approximate the integration. The most widely used approximation methods based on this concept are the First-Order Reliability Method (FORM) (Hasofer and Lind 1974;Hohenbichler and Rackwitz 1982;Der Kiureghian and Dakessian 1998) and second-order reliability methods (SORM) (Breitung 1984;Tvedt 1990;Tvedt 1984;Zhao and Ono 1999a;Hohenbichler et al. 1987;Köylüoǧlu and Nielsen 1994). In these methods, a transformation is first performed from the original distributed variables to independent standard normal variables. Thereafter, the MPP, which has the highest probability density at the limit state in the standard normal space, is found by an optimization procedure. This optimization necessitates the iterative evaluations of the generally computationally expensive limit-state function g(X), and specialized algorithms aiming at maintaining the convergence efficiency have therefore been developed (Wu and Wirsching 1987;Liping and Grandhi 1994). In FORM, the limit-state function is linearized at the MPP and the probability of failure is directly computed from the distance from the origin to the MPP (Mansour and Olsson 2018), denoted as the first-order reliability index. However, approximating the probability of failure using a linear function can result in high inaccuracies, which limits the applicability of FORM for engineering applications (Eamon et al. 2005;Schuëller et al. 2004). In SORM, the limit-state is approximated by a quadratic function using a second-order Taylor expansion at the MPP, in order to achieve higher accuracy in the probability estimates. In the traditional SORMs proposed by Breitung (Breitung 1984) and Tvedt (Tvedt 1990;Tvedt 1983), a rotational transformation is performed after the second-order Taylor expansion (Madsen et al. 2006). The general quadratic function is thereafter approximated by a paraboloid, by ignoring the last row and last column in the transformed Hessian matrix (Choi et al. 2007;Mansour and Olsson 2014). Finally, the probability of failure is analytically computed from the curvatures of the parabola and the first-order reliability index based on asymptotic formulas (Breitung 1984;Zhao and Ono 1999a;Tvedt 1983).
Various modifications have been made to the traditional SORMs proposed by Breitung and Tvedt. Hohenbichler and Rackwitz (Hohenbichler and Rackwitz 1988) improved the accuracy of reliability analysis when all curvatures at the MPP have the same sign. Cai and Elishakoff (Cai and Elishakoff 1994) proposed a refined second-order reliability method used for curvatures with different signs. Koyluoglu and Nielsen (Köylüoǧlu and Nielsen 1994) developed three approximations for the SORM integrals, applied for positive, negative, and mixed curvature signs at the MPP. Zhao and Ono (Zhao and Ono 1999b) provided an empirical closedform formula which does not necessitate the rotational transformation involved in the aforementioned SORMs.
Recently, SORMs based on the general quadratic approximation of the limit-state have been proposed. These methods take into account all elements of the Hessian matrix as opposed to the parabolic SORMs. In this category, Mansour and Olsson (Mansour and Olsson 2014) proposed closed-form probability expressions based on asymptotic expansions of the sum of non-central chi-squared distributions. Hu and Du (Hu and Du 2018a) and Huang et al. (Huang et al. 2018) integrated the saddlepoint approximation method with SORM to further improve the accuracy of the reliability analysis. Lee et al. (Lee et al. 2012) proposed a general quadratic SORM using non-central chi-squared distributions which was later improved by Park and Lee (Park and Lee 2018).
In numerous engineering applications, reliability may vary with respect to time if the associated limit-state function involves dynamic parameters, such as time-dependent load processes and deterioration of material performance. The timedependent reliability analysis methods can roughly be classified into three groups. The first group of methods convert time-dependent reliability analysis into time-independent reliability analysis using the extreme value of the time-dependent limit-state function. Thereafter, static reliability analysis methods such as FORM and SORM can be applied. For instance, the time-variant reliability analysis method based on stochastic process discretization method (TRPD) (Jiang et al. 2014) discretizes the random process into random variables, and transforms a time-variant reliability problem into a conventional time-invariant system reliability problem. The extreme value response method (Wang and Wang 2012) converts the time-dependent reliability analysis to a timeindependent one by building a Kriging model between the design variables and the time when the system response approaches its extreme value. Hu and Du (Hu and Du 2015) proposed a mixed efficient global optimization method employing adaptive Kriging-Monte Carlo simulation to build the surrogate model of the extreme value function as well as a second-order approximation to the extreme value using sequential efficient global optimization . Singh and Mourelatos (Singh and Mourelatos 2010) introduced the composite limit-state function and identified the significant MPP to calculate the time-dependent reliability. Du (Du 2014) uses the envelope of instantaneous limit-state functions, which is explicitly related to time, to deal with timedependent mechanism reliability. The second group includes methods based on the Rice's formula, in which the key idea is to calculate the upcrossing rate or outcrossing rate using a Poisson distribution assumption. For example, Lutes and Sarkani (Lutes and Sarkani 2003) derived an analytical outcrossing rate for Gaussian stationary processes as well as for general Gaussian stochastic processes (Lutes and Sarkani 2009). Singh and Mourelatos ) developed a lifecycle cost optimization approach using the outcrossing rate as the failure rate, and proposed an importance sampling method to approximate the outcrossing rate (Singh et al. 2011).  derived analytical upcrossing and downcrossing rates followed by a numerical procedure to integrate the two rates to obtain the kinematic reliability. Hu and Du (Hu and Du 2012) employed the upcrossing rate method to calculate the time-dependent reliability of the hydrokinetic turbine blade over its design life period. Moreover, numerous empirical modifications (Hagen and Tvedt 1991;Engelund et al. 1995;Streicher and Rackwitz 2004; have been made to the upcrossing rate methods. The third group contains the simulation-based surrogate modeling methods. These methods generally build a surrogate model to replace the limit-state function by evaluating the limit-states at a number of points predefined through Design of Experiment (DoE) (Condra 2018;Holman 2012;Lundstedt et al. 1998). Thereafter, MC simulation can be performed on the surrogate model to compute the reliability. The typical methods include polynomial chaos expansions (PCE) (Hu et al. 2013;Li and Ghanem 1998), Kriging or Gaussian process-based methods (Hu and Mahadevan 2016a;Martin and Simpson 2005;Zhang and Du 2015;Zhu and Du 2016), and machine learning methods (Papadrakakis and Lagaros 2002;Deng et al. 2005).
Besides time, the limit-state function may also involve interval variables that cannot be quantified by the probabilistic model. This is often the case when the information of variables is limited or is expensive to obtain in the early design stage. Therefore, the limit-state function may include mixed variables with both random and interval variables. Several techniques have been developed for reliability analysis with mixed variables. For example, Kang et al. (Kang and Luo 2010;Luo et al. 2009) developed a probability and a convex set mixed model to search for the worst-case point (WCP) and the MPP. A sequential single-loop (SSL) procedure was proposed by Du et al. (Du et al. 2005;Guo and Du 2010) to perform reliability-based design optimization (RBDO) when interval variables are considered. Li et al. (Li et al. 2019) combined the sequential sampling strategy with convex set mixed models to improve the efficiency of RBDO when random variables and interval uncertain variables coexist. Furthermore, Li et al. (Li et al. 2020) proposed a sequential iterative strategy to calculate the worst-case reliability after transferring time-dependent reliability problem into a static reliability one using TRPD. When the worst-case scenario is considered, the reliability problem with mixed variables can be converted to time-dependent problem (Hu and Du 2016).
Although methodologies have been developed for timedependent and mixed variables reliability problems, the applicability of SORM for these problems is not mature. However, since mixed variables problems can be converted to timedependent reliability, which in turn can be converted to a static reliability problem, the accuracy of SORM for static reliability evaluation is crucial. Therefore, this comparative study is limited to static problems considering random variables with given probability distributions. As introduced earlier, although a number of SORMs are available for static reliability analysis, their applicability and relative merits have not yet been investigated in a full-scale comparative study. In this paper, the performance of commonly used SORMs is assessed using a large number of quadratic mathematical problems categorized based on the problem scale, curvatures at the MPP, first-order reliability index, and limit-state contour. Each category is further divided into the following classes: small and large for the problem scale; small, large, positive, negative, and mixed for the curvatures; small and large for the first-order reliability index and parabolic or full quadratic for the limit-state contour. Using this detailed testing scheme, the applicability of the different SORMs is evaluated based on the characteristics of the second-order Taylor approximation of the limit-state function. The assessment of the methods is based on three performance metrics: Capability, Accuracy and Robustness. The Capability is a measure of the ability of a method to compute feasible probabilities, i.e., probabilities between 0 and 1. A high Capability is crucial in engineering design, especially in the context of RBDO (Choi et al. 2007;Mansour and Olsson 2016), in which a non-real or negative probability in any iteration may result in non-convergence. The Accuracy is quantified based on the mean relative error with respect to exact reliabilities. The Robustness metric is a measure of the variability in the performance of the method and is therefore quantified based on standard deviation of relative errors.
The remainder of this paper is organized as follows. Section 2 gives a review of the most representative secondorder reliability methods, see Fig. 1. This includes three relatively new methods, Mansour and Olsson's method, the second-order saddlepoint approximation (SOSPA), and Park and Lee's method, as well as a number of traditional methods, i.e., Breitung's method, Tvedt's method, Hohenbichler and Rackwitz's method, Koyluoglu and Nielsen's method, Cai and Elishakoff's method, and Zhao and Ono's method. The performance of these SORMs is evaluated in Section 3 based on a wide range of quadratic limit-states. In Section 4, three engineering examples are solved in order to study the performance of SORMs for engineering examples and quantify the error from approximating the limit-state with a quadratic function at the MPP.

Second-order reliability methods
In SORM, the original random variables X are first transformed into standard normal uncorrelated variables U. This can be performed using the Rosenblatt transformation (Rosenblatt 1952;Hohenbichler and Rackwitz 1981), the Nataf transformation (Nataf 1962;Liu and Der Kiureghian 1986) which is used in this work, or alternatively using the recently developed analytical transformation proposed by Lu et al. (Lu et al. 2020). The limit-state function g u (U) = g(X) is approximated with a second-order Taylor expansion q u (U) at the Most Probable Point of failure (MPP), denoted by u * , according to where ∇g u (u * ) is the gradient vector at the MPP and ∇ 2 g u (u * ) is the Hessian matrix. The MPP is defined as the point with the highest probability density on limit-state boundary g u (U) = 0 and therefore g u (u * ) = 0. It corresponds to the shortest distance from the origin to the limit-state and is therefore obtained by solving a constrained non-linear optimization problem according to All SORMs are derived based on the second-order Taylor approximation of the limit-state according to (3). Next, we will briefly review the commonly used methods, i.e., Breitung's method, Tvedt's method, Hohenbichler's method, Koyluoglu's method, Cai's method, Zhao's method, Mansour's method, SOSPA, and Park's method. These methods fall into two categories, parabolic SORMs and general quadratic SORMs.

Parabolic second-order reliability methods
An exact solution to the probability of failure cannot be directly derived from (3) due to the second-order terms. Therefore, two variable transformations are performed in traditional SORMs. The first transformation Z = TU, where T is obtained by a Gram-Schmidt orthogonalization (Choi et al. 2007), rotates the coordinate system so that last coordinate Z n coincides with the vector u * from the origin to the MPP. The limit-state function is rewritten as where β = ∥u * ∥ is the first-order reliability index, z * = [0, 0, ⋯, β] T is the MPP in the Z-space and is the transformed Hessian matrix. The second equality in (6) and the columns of T ′ ∈ ℝ (n − 1) × (n − 1) are the eigenvectors of e H with corresponding eigenvalues k i (Mansour and Olsson 2014). The limit-state function in (5) can therefore be written as It is noted from (7) that γ i ¼ n define the components of the last column of the symmetric Hessian matrix ∇ 2 g y (y * ) evaluated at the MPP. The traditional parabolic SORMs further approximates q y (Y) by a paraboloid p(Y) by setting γ i = 0, ∀i and λ = 0, i.e., by neglecting the last row and column of the Hessian matrix, resulting in The expression in (8) is a function of the distance from the origin to the MPP in U-space defined by the first-order reliability index β = ∥u * ∥, the main curvatures k i of the limit-state function g u (u) at the MPP, and the number of random variables n. An illustration of p y = 0 in the y 1 y 2 y n -space for the case k 1 > 0 and k 2 < 0 is shown in Fig. 2 with marked MPP and first-order reliability index. As can be seen, the projected curve on the y i y n -planes is concave when k i > 0 and convex when k i < 0. Although the expression is substantially simpler Fig. 1 Timeline with studied SORMs in the present review paper Z. Hu et al. than the general quadratic approximation according to (3), there is no exact solution for the reliability integral defined in (1). In Sections 2.1.1-2.1.6, the most common approximations of the multidimensional reliability integral based on the above parabolic approximation are reviewed.
2.1.1 Breitung's method (Breitung 1984) The probability of failure is computed based on (2) with the integration domain defined using the parabolic approximation p(Y) in (8), i.e., A variable transformation Y¼β −1 Y results in The integral in (10) is of the form ∫ D exp β 2 F yÞ ð ÞdyÀ for which an asymptotic expansion for β → ∞ and two times continuously differentiable functions F yÞ ð has been derived (Bleistein 1986). The domain D is defined by the parabolic equation −y n þ 1 þ 1 2 ∑ n−1 i¼1 βk i y 2 i which has curvatures βk i and a Most Probable Point of failure defined by y * ¼ 1. Direct application of these properties in the asymptotic expansion (Bleistein 1986) results in Breitung further approximates the above expression using the asymptotic relation 2π where Φ(·) is the standard normal Cumulative Distribution Function (CDF), yielding Based on the above derivation, it is clear that the Breitung formula is applicable for large β. This follows from the approximation of the integral in (10) by the asymptotic expansion in (11) as well as the approximation of the latter by (12). Furthermore, the Breitung formula is only applicable for k i > −1/β, ∀i, as is seen by (12).
2.1.2 Tvedt's method (Tvedt 1990;Tvedt 1983) In Tvedt's method, an asymptotic expansion based on the reliability integral is derived. Since the safe set is given by where ϕ(·) is the standard normal PDF. Tvedt derived a threeterm approximation for this equation by a power series expansion in terms of 1 2 ∑ n−1 i¼1 k i Y 2 i , ignoring terms of orders higher than two (Madsen et al. 2006). The resulting closed-form expression for p f is given by the three-term approximation where in which Re(·) denotes the real part of a complex number. As can be seen, the first term in (14) is equal to the probability formula proposed by Breitung, i.e., T 1 = p f,Breitung . The second and third terms, T 2 and T 3 , can be interpreted as correctors to Breitung's formula aiming at increasing the accuracy for moderate values of β.
2.1.3 Hohenbichler and Rackwitz's method (Hohenbichler and Rackwitz 1988) Similar to Tvedt, Hohenbichler and Rackwitz derived a closed form expression based on a Taylor expansion in term of 1 2 ∑ n−1 i¼1 k i Y 2 i . The probability integral p f is first written as in (13) where the inner integral is replaced by Second-order reliability methods: a review and comparative study with the integration bounds corresponding to the failure set p(Y) < 0, i.e., Y n < −β− 1 2 ∑ n−1 i¼1 k i Y 2 i . A Taylor expansion of the right-hand side of (16) in terms of 1 2 ∑ n−1 i¼1 k i Y 2 i is given by which follows from the first-order Taylor expansion Replacing the inner integral according to (16) by its firstorder Taylor approximation according to (17) and inserting into the probability integral yields This probability formula aims at improving the reliability estimate for moderate β compared to Breitung's formula. It should be noted that the above expression is asymptotically 2.1.4 Koyluoglu and Nielsen's method (Köylüoǧlu and Nielsen 1994) In Koyluoglu and Nielsen's method, the integration domain Thereafter McLaurin series in such a way that (19) can be integrated analytically. Three formulas are proposed for different signs of curvatures as follows.
1. For the case when all curvatures are positive, i.e., k i > 0, ∀ i = 1, 2, …, n − 1 the probability of failure is given by 2. For the case when all curvatures are negative, i.e., k i < 0, ∀ i = 1, 2, …, n − 1, the probability of failure is given by 3. For the case of both positive and negative curvatures k i > 0, ∀ i = 1, 2, …, m − 1 and k i < 0, ∀ i = m, m + 1, …, n − 1, a one-term approximation of the probability of failure is given by 2.1.5 Cai and Elishakoff's method (Cai and Elishakoff 1994) Similar to Koyluoglu and Nielsen's method, in Cai and Elishakoff's method, the inner integral in (13) is first divided into two parts according to where v = y n − β. The exponential term exp −βv− v 2 2 is thereafter expanded into a Taylor series at v = 0. The integral in the right-hand side of (23) can therefore be integrated analytically, which results in a three-term approximation of the probability of failure according to Z. Hu et al.
2.1.6 Zhao and Ono's method (Zhao and Ono 1999a) Rather than using the main curvatures shown in (8), in Zhao and Ono's method, the limit-state function is approximated by a rotational parabolic function with curvatures equal to the average of the main curvatures, i.e., where K s ¼ ∑ n−1 i¼1 k i is the sum of the main curvatures. The latter can be computed as where tr(∇ 2 g u (u * )) is the trace of the Hessian matrix evaluated at the MPP and α ¼ − ∇g u u * ð Þ ∇g u u * ð Þ is the unit normal vector to the MPP. The limit-state function p(Y) = 0 in (26) therefore becomes a combination of a normal random variable −Y n + β and a random variable 1 2 K s n−1 ∑ n−1 i¼1 Y 2 i with a central chi-square distribution with n − 1 degrees of freedom. The probability of failure is given by where f χ 2 n−1 t ð Þ is the PDF of the central chi-square distribution and R ¼ n−1 K s is the average main curvature radius. Since a closed form for (28) cannot be obtained, Zhao proposed two empirical closed-form expressions depending on the sign of K s according to As can be seen from (27), Zhao's method directly uses the gradient and Hessian in u-space and bypasses the transformations and eigenvalue analysis in (5)-(8).

General quadratic second-order reliability methods
Recent research efforts have been focusing on developing closed-form reliability expressions based on the general quadratic approximation q u (U) according to (3), which can be written in standard form as Three methods are studied and these are presented in Sections 2.2.1-2.2.3. These methods use all components of the Hessian matrix as opposed to the parabolic SORMs presented in Section 2.1. It is noted that, inserting u * = βα where α = − ∇ g u (u * )/‖∇g u (u * )‖ in (32), the SORM limit-state according to (31) is commonly written as

2014)
In Mansour and Olsson's method, a variable transformation e U ¼ D −1 U removing the second-order cross terms in (31) is first performed, where D is an orthogonal matrix whose column vectors are the eigenvectors of C, and e U ¼ e U 1 ; e U 2 ; …; e U n is a n-dimensional vector with independent standard normal random variables. Thus, the limit-state function becomes where e b¼D T b and e C¼D T CD . Since e C is diagonal, (34) can be written as.
Second-order reliability methods: a review and comparative study The probability of failure is thus given by The function h(·) in (36) is a linear combination of noncentral chi-squared variables for which Konishi (Konishi et al. 1988) proposed asymptotic expansions of its probability distribution. Based on Konishi's work, Mansour and Olsson proposed two probabilities of failure formulas depending on the limit-state contour.
For non-elliptic limit-states, i.e., when e C ii have different signs as well as for the paraboloid case according to (8) regardless of the signs of k i , the probability of failure is approximated by the closed-form expression where H i (·) is the ith probabilists' Hermite polynomials and For the elliptic case, i.e., when all e C ii have the same sign, the probability is approximated by p f = P when sgn e C ii Á s < 0 and p f = 1 − P when sgn e C ii Á s > 0, where sgn(·) is the signum function, and 2.2.2 Second-order saddlepoint approximation (Hu and Du 2018a) In the Second-Order Saddlepoint Approximation (SOSPA) proposed by Hu and Du, the Cumulant Generating Function (CGF) of the limit-state in standard normal space is first analytically expressed. Thereafter, the saddlepoint approximation is applied to estimate the probability of failure.
To analytically derive the CGF, the cross terms in the secondorder Taylor expansion in (31) are first eliminated with the same variable transformation e U¼D −1 U as in Mansour and Olsson's method. Equation (34) is then rewritten as where a ei =a/n, which is further rewritten as The CGF of q u e e U is given by (42) according to (Du and Zhang 2010).
Z. Hu et al. where The CGF can equivalently be written in a more compact form as (Huang et al. 2018).
Once K q (t) is available, the saddlepoint t s is found by solv- respect to t. According to the Lugannani and Rice's (Lugannani and Rice 1980), p f is computed by In (46), sgn(t s ) = + 1, − 1 or 0, depending on whether t s is positive, negative, or zero, respectively, and K ″ q t ð Þ is the second derivative of K q (t) with respect to t.

Park and Lee-convolution integral (Park and Lee 2018)
Similar to Mansour and Olsson's method as well as SOSPA, Park and Lee's method starts with the variable transformation e U¼D −1 U in order to remove the second-order cross terms in (31). Using (33), the limit-state is given by where and λ i is the ith eigenvalue of the Hessian matrix ∇ 2 g u (u * ). Equation (47) can therefore be written as where the first and second summation extend over all positive integers i satisfying λ i ≠ 0 and λ i = 0, respectively. This is further rewritten as where (51) is computed by convolution according to Fig. 3 Flowchart of first-and second-order reliability methods Second-order reliability methods: a review and comparative study where the PDF of Q i is expressed using the PDF of W i as The convolution integral in (52) is computed numerically. Both the efficiency and accuracy of the method is therefore directly affected by the chosen integration bounds and the discretization step used in the numerical convolution. A smaller discretization step and larger integration bounds may be required for the computation of small probabilities with sufficient accuracy. In the remainder of this paper, the integration bounds [−50,50] and [0,50] are chosen for the normal and chi-squared variables, respectively, and the discretization step is 0.001.

Summary of second-order reliability methods
The different steps involved in the most widely used SORMs are summarized in the flowchart in Fig. 3. Given an arbitrary and implicitly given limit-state function g(X) where X is a vector of arbitrarily distributed and correlated random variables, the first steps in all SORMs are a transformation to standard normalized independent space U. The MPP, u * , is thereafter located based on a constrained optimization problem, which is followed by a quadratic approximation q u (U) of the limit-state function g u (U) at the MPP. For all methods, except Zhao's, Mansour's, Park's and SOSPA, see Fig. 3, a variable transformation U → Z aiming at rotating the coordinate system so that the last coordinate coincides with the vector u * from the origin to the MPP is performed. The leading submatrix e H of the Hessian matrix in Zspace is diagonalized and the limit-state function is approximated by a parabola by neglecting γ i and λ in the quadratic function according to (7). The probability of failure is finally computed from the curvatures of the parabola and the first-order reliability index. For Zhao's method, a rotational transformation to Z-space is not necessary since the probability of failure is a function of the average curvatures at the MPP computed in U-space, see (27). Similarly to Zhao's method, in Mansour's method Park's method and SOSPA, a rotational transformation to Z-space is not necessary. These methods are based on a full quadratic approximation of the limit-state function and not on parabolic expressions. The full Hessian matrix 2C in U-space, see (31), is diagonalized, followed by approximate expressions of the probability of failure.

Features of the testing problems
The quadratic testing problems are classified based on the following four representative features, which are summarized in Table 1.
Problem scale: It is characterized by the number of random variables n. The scale is considered to be small when n < 10 and large when n ≥ 10. Curvature: The curvature is classified as positive, negative, or mixed, if the main curvatures are positive, negative, or have different signs, respectively. It is further Large β>2

Contour
Parabolic Hu et al. classified as small or large if the magnitude of all curvatures is smaller or larger than 1/5, respectively. First-order reliability index: It is characterized as small if β ≤ 2 and large if β >2 . Contour type: The problems are classified into parabolic contour and non-parabolic contour where the latter includes ellipse, hyperbola, and mixed contours. The quadratic limit-state function is derived in Section 2, see (7), where λ = 0 and γ i = 0, ∀i, result in a parabolic contour.

Performance metrics
Three performance metrics are used to evaluate the performance of each SORM. These are derived from the relative error in the logarithm of the probability of failure given by in which p fi is the SORM result for the ith problem, and p fi,exact is the result from MC simulation or an analytical solution if available. This expression, which may seem unorthodox as a definition of relative error, is fair with regard to the more demanding cases of extremely small probabilities of failure. Next, three performance metrics are described: Capability, Accuracy, and Robustness. The Capability is the ability of a SORM to achieve a probability of failure p f within the interval [0, 1]. It is defined as the ratio of the number of problems N c fulfilling this requirement and the total number of problems N, i.e.
where I (p fi ∈ [0, 1]) is an indicator function that takes the value 1 when p fi ∈ [0, 1] and 0 when p fi ∉ [0, 1]. The Accuracy is the complementary mean of the relative error defined in (53), evaluated among the N c problems for which the method has been deemed capable, i.e., It should be noted that if the mean relative error mean(ε) of a method is larger than 1, the accuracy is assigned to be 0, which is the lowest possible level.
The Robustness is a quality measure of the method describing its ability to achieve good accuracy for different types of problems. It is measured by the complementary standard deviation of the relative error, evaluated among the N c problems for which the method has been deemed capable, i.e., It is noted from (57) that if the standard deviation of the relative error std(ε) of a method is larger than 1, the robustness is assigned to be 0, which is the lowest possible level. An illustration of (55) and (57) is shown in Fig. 4. In the computation of the mean and standard deviation of the relative error used in the accuracy and robustness metrics, statistical outliers corresponding to values outside the 5th and 95th and percentile are disregarded.

Results and comparison
In this section, mathematical quadratic testing problems are selected to compare the performance of various second-order reliability Second-order reliability methods: a review and comparative study methods based on (7). Seven problem sets are defined based on the four features in Table 1. In Table 2, a summary of the testing problems used in this comparative study is presented. The total number of testing problems is 5520.

Problem Set 1
In Problem Set 1, all curvatures k i ¼ 1 R > 0 are equal and positive, and the scale of the problem is small with n = 8. The probability of failure p f = Pr[q y = 0] is computed for a wide range of combinations of curvature radii R = 1,1.1, ..10 and first-order reliability indices β = 0.1,1, 2, 3, 4. Three contour types are considered: parabolic, elliptic, and hyperbolic, see Table 2.
For the parabolic case (γ i = 0, ∀ i, and λ = 0, see (7)), the limit-state function is a function of the first-order reliability index β and the curvature radius R according to The contour for q y = 0 is shown in Fig. 5a for β = 1 and curvatures R = 10 and R = 1, where hatched region represents the failure domain q y = 0. It should be observed that the problem is rotationally parabolic, i.e., the contour in Fig. 5a is identical in  all directions y i , i = 1, …, n − 1. A subset of the results is shown in Fig. 6a computed using the methods presented in Section 2 as well as MC simulation, where the latter is considered as the exact result. As can be seen, for a small reliability index β = 1, all methods have good accuracy when the curvature radius is large, while the performance strongly deteriorates for all methods except SOSPA, Zhao's, Mansour's, and Park's as the curvature radius decreases. For a large β = 3, Cai's method yields large inaccuracies in the whole range of studied R and Mansour's method results in large deviations from the exact results for small curvature radii. The elliptic case is studied by adding the term 1 2 λ y n −β ð Þ 2 to the rotationally parabolic equation according to (59) with R > 0, i.e., The elliptic case corresponds to (7) with γ i ¼ 0; Ák i ¼ 1 R > 0, and λ > 0 . As can be seen from Fig.  5b, increasing λ > 0 yields a decrease of the area of the elliptic failure region.
The problem is solved for λ = 0.5 and λ = 1, where the results for the latter case are presented in Fig. 6b when β = 1.
As can be seen, the parabolic SORMs (ref. to Section 2.1) mostly result in substantially higher probability of failure values compared to exact results, since they approximate the elliptic limit-state by a parabola by neglecting the influence of λ. It is also observed that although Mansour's method captures the effect of λ on the probability of failure, its accuracy is lower than SOSPA and Park's method.
The hyperbolic contour is studied by adding the term y n −β ð Þ∑ n−1 i¼l γ i y i to the rotationally parabolic equation according to (59) with R > 0, where it is assumed that γ i ≡ γ, ∀ i, i.e., This corresponds to (7) with λ = 0 and k i ¼ 1 R > 0. The failure domain is shown in Fig. 5c for γ = 0.2 and γ = 0.5. It can be seen that, as γ increases, the second branch of the failure region approaches the origin from the fourth quadrant. However, for γ ≤ 0.5, the contribution of this second branch to the failure probability is small. The computed probability of failures for γ = 0.5 and β = 1 is shown in Fig. 6c. As can be seen, the parabolic SORMs (which sets γ = 0) yield lower p f estimates compared to both exact probabilities (MC) and the full quadratic SORMs (Mansour, SOSPA, and Park).

Problem Set 2
In Problem Set 2, the effect of problem scale (3 ≤ n ≤ 40) on the performance of the different SORMs is investigated, see Table 2. Parabolic, elliptic, and hyperbolic contours are considered for β = 0.1,1, 2, 3, 4 where all curvatures k i = 1/R, R = 10, are equal and positive.
For the parabolic case, the limit-state function according to (59) is applied. The results for β = 2 and β = 4 are presented in Fig. 7a. It is noted that Cai's method yields non-real results for large-scale problems. These results, which adversely affect the capability metric (see Section 3.2), are not marked in Fig. 7a. The performance for Tvedt's method also strongly deteriorates as n increases. The highest accuracy is achieved by SOSPA, Mansour's, and Zhao's methods followed by Park's and Hohenbichler's method. These methods show good performance also for the very challenging small probability estimates (≈10 −8 ).

Problem Set 3
In Problem Set 3, all curvatures k i = 1/R < 0 are equal and negative, and the problems are small scaled with n = 8. The probability of failure is computed for all combinations of R = − 1, − 1.1, …, − 10 and β = 0.1,1, 2, 3, 4. Two contour types are considered: parabolic and hyperbolic, see Table 2. The parabolic case is studied using (59) with R < 0. The contours for R = 1 and R = − 10 are shown in Fig. 8a when β = 1. A subset of the results (β = 2and β = 4) are shown in Fig. 9a. Most notably is that Cai's method has an overall superior performance for negative curvatures compared to positive ones, especially in terms of capability, as can be seen by comparing its results to Fig. 6a. It is also noted that Breitung's, Tvedt's, Hohenbichler's, and Koyluoglu's methods deviate from the exact probabilities for small |R|. The best accuracy is observed for Zhao's, Mansour's, SOSPA, and Park's method which perform well for both large and small curvatures as well as large and small β.
The hyperbolic case is studied using (60) with R < 0 for λ = 0.5 and λ = 1 It should be observed that, while R > 0 results in an elliptic contour as is seen in Fig. 5b, the case R < 0 corresponds to a hyperbolic contour, see Fig. 8b. It is noted that the parabolic SORMs approximate the failure region using the parabolic contour by assuming λ = 0 which corresponds to the dashed blue line in Fig. 8b. A subset of the results (β = 2 and λ = 0.5) is presented in Fig. 9b.

Problem Set 4
In Problem Set 4, the same types of contour as in Problem Set 3 are studied but the effect of problem scale is investigated instead of the curvature radii. The probability of failure is computed for n = 3, 4, …, 40 and for a fixed negative curvature radius R = − 10.
Results for the parabolic case are shown in Fig. 10a. As can be seen, the highest accuracy is achieved by Zhao's, Mansour's, and SOSPA followed by Park's method. It is also seen that Hohenbichler's and Tvedt's method yield unfeasible probabilities (p f > 1) for large-scale problems. Furthermore, the overall accuracy of Cai's method is superior to Breitung's, Hohenbichler's, and Tvedt's methods.
In Fig. 10b, results for the hyperbolic case corresponding to (60) are shown for β = 2 and λ =1 . All methods except Mansour's, SOSPA, and Park's assume λ = 0. For this particular case, it is observed that the error in Zhao's method is solely due to its approximation of the hyperbolic contour by a parabola (λ = 0). This can be concluded by comparing the results to Second-order reliability methods: a review and comparative study the parabolic case in Fig. 10a for which Zhao's method has a high accuracy in the whole range of n. However, for Breitung's, Hohenbichler's, Cai's, and Tvedt's methods, the highly erroneous results for large-scale problems are predominantly due to inaccuracies in the probability formula as n increases.

Problem Set 5, Set 6, and Set 7
In Problem Sets 5, 6, and 7, the accuracy of the studied SORMs is investigated for parabolic limit-state functions according to (59). In Problem Set 5 and Set 6, the curvatures are equal in all directions but with opposite signs. In Problem Set 7, the case of curvatures with different magnitudes and signs is investigated.
In Problem Set 5, a small-scale problem with n = 8 and k i = (−1) i /R is studied for a wide range of R and β = 0.1,1, 2, 3, 4, see Table 2. The limit-state is given by The corresponding contours are shown in Fig. 11a for R = 10 and R = 0.5 when β = 1. A subset of the results, β = 2 and β = 4, are presented in Fig. 12a for R = 1,1.1, …, 10. Most notably is that the performance of Zhao's method, although satisfactory in many of the previous studied cases, yields low accuracy in the present case. The highest accuracy is observed for SOSPA, Mansour's, and Park's methods followed by Cai's method. For small curvature radii, Breitungs's and Hohenbichler's methods yield probabilities of failure larger than unity (when β = 2) and Tvedt's method results in non-real probabilities. In Problem Set 6, the influence of problem scale n is studied for a fixed curvature radius R = 10, see Table 2. The results for β = 2 are shown in Fig. 12b. It is noted that p f is not monotonically increasing with n. When the problem scale is even, adding another variable corresponds to the addition of a negative curvature −1/R to the problem. However, when the problem scale is uneven, a positive curvature 1/R is added when the problem scale is increased by one. This results in the non-monotonically increasing pattern seen in Fig. 12b. It is observed that Koyluoglu's and Zhao's methods yield results that are close or identical to FORM, respectively. These methods therefore have low accuracy for the case of equal curvatures of opposite signs, regardless of the problem scale. Good accuracy is however achieved with Breitung's, Cai's, Mansour's, SOSPA, and Park's methods, while Tvedt's and Hohenbichler's methods overestimate the probability of failure.
In Problem Set 7, problems with curvatures that have different magnitudes and signs in different directions are investigated. The limit-state function is given by where −1 < a < 0, see Table 2. In Fig. 11b, contours of the limit-state in (63) are shown for a = − 0.7 and a = − 0.3. As can be seen, a = − 0.3 results in a high nonlinearity in the y 1 y 8 -space but an almost linear contour in the y 5 y 8 -space. In Fig. 12c, computed probabilities of failure are shown as a function of a for β = 2 and β = 4. Overall, the best performance is achieved by SOSPA, Mansour's, and Park's methods followed by Cai's method. It is also noted that, although the accuracy of Koyluoglu's method is low, its capability is high. The lowest capability is observed for Breitung's, Hohenbichler's, and Tvedt's methods with negative or nonreal probabilities, which are not marked in the figure.

Performance evaluation
The performance metrics defined in Section 3.2, i.e., Capability, Accuracy, and Robustness, are evaluated from the results of the test problems in Table 2. In Sections 3.4.1-3.4.3, the influence of curvature sign and magnitude, problem scale, and first-order reliability index as defined in Table 1 is presented based on the results for the parabolic limit-states. In Section 3.4.4, the overall performance is presented based on the results for both parabolic and general quadratic limit-states. The reader is here reminded that the Capability is the proportion of problems for which the method yields a feasible probability and that the Accuracy and Robustness are computed for problems with feasible probabilities. It should also be noted that the Robustness metric measures the consistency in the performance of the method. Therefore, a high Robustness simply signifies that the accuracy metric is representative for all studied problems.

Influence of curvature sign and magnitude
In the following, the Capability, Accuracy, and Robustness of the SORMs are obtained based on curvature signs and curvature magnitude. For a better visualization, the initial letter of each method is used in the metrics plot. In Fig. 13a, the performance metrics for the positive curvature case are presented. As can be seen, the best performance is achieved by SOSPA followed by Zhao's and Park's methods. Mansour's, Hohenbichler's, and Koyluoglu's methods also have a high performance where the two latter methods use the same expression for the probability computation, see Section 2.1.3 and Section 2.1.4. The lower Capability of Cai's method for positive curvatures is in accordance with the results in Fig. 6a and Fig. 7a.
It is observed that the performance of Breitung's, Tvedt's, Hohenbichler's, and Koyluoglu's methods deteriorate for negative curvatures, see Fig. 13b. This is in contrast to Cai's method which has a substantially higher capability compared to the positive curvature case. Similarly to Cai's method, Mansour's method has a slightly higher performance for negative curvatures compared to the positive curvature case, especially in terms of robustness.
For the mixed curvature case, see Fig. 13c, the best performance is observed for Mansour's, SOSPA, and Park's methods followed by Cai's method. It is also seen that the performance of Zhao's method is lower than for both the positive (Fig. 13b) and negative curvature case (Fig. 13c), since the method only considers the average of curvatures and thus the effect of sign is ignored (see Section 2.1.6). For Breitung's, Tvedt's, Hohenbichler's, and Koyluoglu's methods, the overall performance is lower than for the positive curvature case but substantially superior to problems where all curvatures are negative.  Figure 14 further illustrates the performance of SORMs evaluated for the curvature magnitude instead of sign, where small and large curvatures are defined by |k i | ≤ 1/5 and |k i | > 1/5, ∀ i, respectively. It is noted that all methods have a satisfactory performance for small curvatures (Fig. 14a). However, all methods except SOSPA and Park's show a lower performance for larger curvature size (Fig. 14b). It is also observed that a very high capability for large curvatures is achieved by Koyluoglu's and Zhao's methods.

Influence of problem scale
The three performance metrics evaluated among all parabolic problems are shown for small (n < 10) and large (n ≥ 10) scale problems, in Fig. 15a and b, respectively. The results show that the accuracy of all SORMs except Zhao's, Mansour's, and SOSPA methods is slightly higher for small-scale problems. It is instead noted that the overall performance of Zhao's and Mansour's methods improves with increased problem scale. The slightly lower performance of Park's method for largescale problem is attributed to errors in the 1D convolution integrals that accumulates as n increases. Overall, the influence of problem scale on the performance of studied SORMs is relatively small in comparison to the effect of curvature size and sign presented in Section 3.4.1.

Influence of first-order reliability index
In Fig. 16, the performance is presented for problems with small (β ≤ 2) and large (β > 2) first-order reliability indexes. It is noted that the accuracy and robustness of Breitung's method improve for large β, as is seen by comparing Fig. 16a and Fig. 16b. A small improvement is also observed for Tvedt's method. It is however noted that Zhao's, Cai's, and Mansour's methods have a lower performance for larger β, where the largest decrease in terms of capability is observed for Cai's method. The best performance for both small and large β is observed for SOSPA followed by Mansour's and Park's methods.

Overall performance for parabolic and general quadratic limit-states
In Fig. 17a, the overall performance for all studied parabolic problems is shown. In Fig. 17b, the performance for all studied problems, including parabolic and non-parabolic limit states, is presented. As can be seen by comparing the performance of Breitung's method to Tvedt's and Hohenbichler's methods, the corrections proposed by the two latter methods result in a higher accuracy but to the price of a lower capability. It is also seen that the best capability is achieved by Koyluoglu's, Zhao's, SOSPA, and Park's methods. The best accuracy and robustness is seen for SOSPA, Mansour's, and Park's methods. Overall, these three most recent methods have the best performance for both parabolic and general quadratic limit-states.

Computational efficiency
The computational cost of reliability assessment generally has two main contributors: 1) The computation of the full-quadratic approximate limit-state function ((3)), where the MPP, as well as the gradient and Hessian at the MPP, are determined from the true limit-state function (see Section 4.6.1) and 2) the computation of the probability of failure based on the quadratic limit-state approximation. In this section, the second contribution is discussed.
In order to compare the computational efficiency of different SORMs, the CPU time corresponding to the probability computations using the limit-state function is measured. This corresponds to (7) with β = 1, k i = 1/10, γ i = 0, ∀i, and λ = 0.5, i.e., Problem Set 2 (see Table 2). The CPU [Intel(R) Core(TM) i5-7400 (4Core, 3.00GHz)/8GB RAM] time is measured in MATLAB for n = 10, n = 20, and n = 30. The results are presented in Table 3. The MC computed probabilities with confidence interval and corresponding execution times are also added in the table. The absolute error in the MC computed probabilities is given by ε MC ≈2 Â ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi p f 1−p f ð Þ=N MC p (Du and Chen 2000;Mansour and Olsson 2014) with 95% confidence, where N MC is the number of function evaluations on the quadratic limit-state. The allowable error is set to ε MC = p f /10 2 and the required number of function evaluations is given by solving for N MC , yielding N MC ≈ 4.10 4 (1 − p f )/p f.
As can be seen, for all parabolic SORMs (Section 2.1) as well as the full quadratic SORM proposed by Mansour and Olsson (Section 2.2.1), the computational cost to compute the probability of failure given the quadratic approximation of the limit-state is negligible, even for large number of variables.
The computational effort needed in SOSPA is due to the numerical solution of K 0 q t ð Þ ¼ 0 to find the saddle point t s , where K q t ð Þ ¼ ∑ n i¼l K qi t ð Þ and K qi (t) are given by (43). As can be seen from Table 3, this effort increases as the number of variable increases, but the efficiency remains high.
In the SORM proposed by Park and Lee, see Section 2.2.3, the probability of failure is computed using a numerical convolution integral. In this work, the overall performance of the method has been evaluated using the integration bounds [−50,50] and [0,50] for normal and chi- Table 3 Comparison of execution time necessary to compute probability of failure based on the quadratic limit-state squared variables, respectively, and the discretization step 10 −3 . As has been shown, the performance of the method using these parameters is high. In Table 3, it is also seen that the computational efficiency is higher than SOSPA.
However, for small probabilities and large number of random variables, the accuracy of SOSPA is higher than all other SORMs, see Table 4. It should also be noted that the accuracy in Park and Lee's method can increase until convergence to the exact solution if the discretization step is decreased. However, the computational burden also increases non-linearly with the decreased discretization step. This is demonstrated by noting that the execution time increases by a factor 100 when the discretization step is decreased by a factor 10 (Table 3).

Problem characteristic and error metrics
Engineering applications are generally non-quadratic and explicit expressions for the limit-states are typically not available. The error in p f using SORM has therefore three main contributions: 1) the approximation of the true Hessian matrix at the MPP by quasi-Newton methods for increased efficiency; 2) the approximation of the true limit-state by a quadratic function at the MPP; and (3) the inability of the p f expressions to exactly compute the probability content of a quadratic approximation. The last point has been extensively studied in Section 3. In this section, the average relative importance of the three contributions is quantified based on three engineering problems with characteristics summarized in Table 5. The available analytical expression for the limitstate functions in these engineering problems is not explicitly used. Instead, the gradient and Hessian are computed using the finite difference method and the quasi-Newton approximation, respectively.

Error from Hessian approximation using Symmetric Rank-1 update
For the Hessian computation in (3), the Finite Difference Method (FDM) is generally time consuming for engineering applications, especially for high-dimensional problems. Therefore, a commonly used quasi-Newton approach Arora 2004;Lim et al. 2014) is applied in the following approximate Hessian evaluation. The approximate Hessian matrix is updated using the Symmetric Rank-1 (SR1) update at each iteration k during the MPP search (see Section 4.6.1) according to where s (k) and r (k) are the differences of the MPP vectors and gradient vectors between two successive MPP iterations, i.e., To preserve numerical stability, an additional condition is required as (Nocedal and Wright 2006) If (67) is not satisfied, the SR1 update in (65) is skipped, . When k = 0, a zero matrix is used as the initial Hessian matrix, i.e., ∇ 2 g u u * 0 ð Þ ¼ 0.
This Hessian approximation does not include any additional function evaluations on the true limit-state after the MPP Table 4 Computed probabilities of failure given the quadratic limit-state q y ¼ −y n þ 3 þ 1 20 ∑ n−1 i¼l y 2 i þ 1 4 y n −3 ð Þ 2 . MC computed probabilities with confidence interval are presented in Table 3 n  Second-order reliability methods: a review and comparative study search since it simply uses the MPP search history. However, for small errors. These relative errors can therefore be approximately interpreted in an additive sense. In the following examples, it should be emphasized that, even if the true limit-state functions are available, they are only used implicitly in SORM. Therefore, the gradient ∇g u (u * ) and Hessian ∇ 2 g u (u * ) evaluated at the MPP, see (3), are computed using the forward FDM and the SR1 quasi-Newton approximation, respectively.

A slider-crank mechanism
A slider-crank mechanism is shown in Fig. 18a. The mechanism fails if the acceleration a of the slider at the initial time t = 0 exceeds an allowed value a all = 2 cm/s 2 . Thus, the limitstate function of the slider crank mechanism is given by X 1 and X 2 are lengths of beam AB and beam BC, respectively, X 3 is the installation error and ω = π/3 rad/s is the angular velocity of beam AB. All random variables are assumed to be independent with mean μ x = [1, 2, 0.5] T cm and standard deviation σ x = [0.05,0.1,0.1] T cm. The random variable X 1 follows a truncated normal distribution with lower limit 0 cm and upper limit 1.35 cm, X 2 follows a truncated normal (a) A roof truss structure (b) Contour plot of limit-state and failure domain Fig. 19 A roof truss structure with corresponding limit-state function in rotated normalized space Second-order reliability methods: a review and comparative study distribution with lower limit 1.25 cm and X 3 follows a normal distribution.
The contours of the limit-state function in the transformed Y -space are plotted in Fig. 18b. As can be seen, the slider crank system non-linearity is well captured in the vicinity of the MPP by both the parabolic ((8)) and the full quadratic approximation ( (7)). From Table 6, it is further observed that the parabolic approximation is slightly more accurate (ε g = −0.0613% compared to ε g = 0.207%). This example demonstrates that it is not always beneficial to use all elements of the Hessian matrix when locally approximating the limit-state function at the MPP, which can also be observed from Fig. 18b. It is also seen that all SORMs have good accuracy and that Tvedt's method achieves the best performance in terms of relative error ε tot with respect to MC simulations performed on the exact limit-state. The results using FORM are also added for comparison, where it is noted that ε p f ¼ 0 is based on the linear approximation at the MPP.

Roof truss
A small-scale roof truss structure (modified from (Song et al. 2009;Zhao et al. 2015)) with non-normal input variables is solved as a second example, see Fig. 19a. The top chords and compression bars of the truss are reinforced by concrete, and the bottom chords and tension bars are made of steel. A uniformly distributed load q is assumed to be applied on the roof truss, which is equivalent to the nodal load P = ql/4. A failure occurs if the perpendicular deflection of the truss peak node is larger than 3.5 cm. The limit-state function of the truss structure is therefore defined by g X ð Þ ¼ 0:035− ql 2 2 3:81 where A s = 9.82 · 10 −4 m 2 and A C = 0.04 m 2 are deterministic variables, X = [q, l, E S , E C ] is the vector of random variables with means μ x1 = 20 · 10 3 N/m, μ x2 = 12.5m, μ x2 = 1 · 10 11 Pa, and μ x4 = 2 · 10 10 Pa and standard deviations σ x1 = 2.2 · 10 3 N/m, σ x2 = 0.125m, σ x3 = 1.5 · 10 10 Pa, and σ x4 = 3 · 10 9 Pa. The random variables X 1 , X 3 , and X 4 follows a lognormal distribution and X 2 is normally distributed. The contours of the limitstate function in the rotated standard normal Y-space are plotted in Fig. 19b. As can be seen, the roof truss structure is approximately linear in y 2 y 4 -plane. Furthermore, the failure region is overestimated in the y 1 y 4 -plane and underestimated in the y 3 y 4 -planes. From Table 7, it is seen that the relative errors from approximating the limit-state by a parabola ε g = 0.190 and a general quadratic function ε g = 0.145 are small. Furthermore, it is observed that Zhao's method shows the highest accuracy ε tot among all SORMs. However, Cai's and Tvedt's methods have the highest accuracy in terms of the proposed probability of failure formula, as is seen from the low relative error jε p f j ¼ 0:00380.

Cantilever beam with high dimensions
A cantilever beam (Du 2010) is used to investigate the SORMs for engineering problem with high dimensions, see Fig. 20. It is subjected to two external forces F 1 and F 2 , two external moments M 1 and M 2 , and two external distributed loads denoted by (q Ll , q Rl ) and (q L2 , q R2 ). These forces, moments, and distributed loads as well  Table 8. A failure occurs if the deflection of the right tip of the beam is larger than the allowable deflection Δ all = 1.3 cm. The limit-state function of the cantilever beam is defined by g = Δ all − Δ where Δ is the tip deflection given by E = 200 GPa is the Young's modulus, I = wh 3 /12 is the moment of inertia, and M and R are the fixed end bending moment and reaction force given by and The results are presented in Table 9. It is observed that the relative error ε g = 5.20% for the general quadratic SORMs is larger than the relative error |ε g | = 2.48% for the parabolic SORMs. It is further noted that the relative error ε g constitutes the largest contribution to total relative error ε tot in the SOSPA method. The method also has the lowest relative error ε pf = − 0.403%, i.e., the accuracy of the method in computing the probability of failure from the general quadratic limit-state is high. In contrast to SOSPA, Hohenbichler's method achieves the best accuracy ε tot = 0.234% among all the SORMs. The reason is that the error resulting from the parabolic approximation of the limit-state, ε g , is almost corrected by the error resulting from the inaccuracy in the probability formula, ε p f . Finally, it can also be observed that the widely used FORM expectedly yields a large error as a result of the non-linearity of the beam system but also due to the large number of random variables, which is in accordance with the observations in Section 3.

Accuracy and error contributions
The average performance of all SORMs in terms of absolute relative errors from the studied engineering problems is shown in Fig. 21. It is observed that the relative errors |ε g | resulting from approximating the limit-state by a parabola or a general quadratic function at the MPP are low. However, this error is on average larger for the general quadratic case. This demonstrates that the increased local accuracy at the MPP achieved by the general quadratic approximation may be detrimental in terms of accuracy compared to the simpler parabolic approximation. This is for instance clearly observed in Fig. 18b. Overall, the best performance is achieved by Tvedt's, Hohenbichler's, and Cai's methods. It should however be noted that the level of non-linearity of the three studied engineering problems is less severe compared to many of the mathematical problems solved in Section 3. It is also noted that the relative error from the SR1 Hessian approximation, In industrial applications, each function evaluation of the true limit-state may involve expensive FE-simulations. For these cases, a computationally intensive part of the reliability assessment is related to the number of function evaluations needed in the MPP search.
The MPP search is performed using the Hasofer-Lind-Rackwitz-Fiessler (HL-RF) algorithm (Hasofer and Lind 1974;Rackwitz and Flessler 1978). The HL-RL algorithm is an iterative MPP search process, which can be summarized as where u * k ð Þ represents the MPP at the current iteration (k) and ∇g u u * k ð Þ is the corresponding gradient vector at u * k ð Þ computed using the forward FDM as   with step size Δu i = 0.0005 . The MPP search therefore necessitates n + 1 function evaluations at each iteration k. It should be noted that several optimization algorithms can be used for the MPP-search, such as the gradient-projection method (Rosen 1961), sequential quadratic programming (SQP) (Schittkowski 1983), the improved HL-RF method (Liu and Der Kiureghian 1991), and adaptive HL-RF (Yang et al. 2020). It is also noted that, since the forward FDM in (80) may be sensitive to the step size (Kim 2010), it may be beneficial to use the central FDM for the purpose of increased accuracy and the convenient choice of the step size. This would however result in an increased computational cost. For all three studied problems in Sections 4.2-4.4, the termination criteria of the MPP search are chosen as |β (k + 1) − β (k) | < 0.001 and the initial values are  Table 10 for the three studied benchmark problems. It is noted that in SORM, no additional function evaluations are necessary if the quasi-Newton SR1 Hessian approximation is used ( (65)). If the forward FDM is used ( (69)), an additional n(n + 1) is however needed. In Table 11, the simultaneous convergence of the MPP and the Hessian using the SR1 method is shown for the slider crank mechanism example (Section 4.2).

Probability of failure computation
After construction of the second-order approximate limit-state function at the MPP using (3), the probability of failure can be computed using different SORM formulas. This step does not involve any additional

2.264¯¯0
Second-order reliability methods: a review and comparative study function evaluations on the true limit-state function. The computational cost related to the computation of the probability of failure given the explicit quadratic approximation at the MPP for the different SORMs has been studied in Section 3.4.5.

Concluding remarks
This comparative study has addressed nearly 50 years of research in the field of SORM. These methods are generally more efficient than MC simulations, even in the presence of a metamodel, especially in an optimization framework where small probabilities and gradients need to be iteratively computed for a large number of constraints. While early efforts focused on the development of probability of failure expressions based on parabolic approximations of the limit-state at the MPP, modern approaches are capable of addressing higher non-linearities based on full quadratic approximations of the limit-state.
The progress achieved in the field has been quantified using the three metrics: Accuracy, Capability, and Robustness. Overall SOSPA method appears to be the most promising method for reliability assessment, regardless of the size or sign of the curvatures at the MPP, the problem scale, or the magnitude of the first-order reliability index. It has however been noted that the use of a general quadratic approximation of the limit-state, based on the full Hessian at the MPP, may at instance yield lower accuracy in terms of computed probability compared to the simpler parabolic approximation. In terms of simplicity, Breitung's, Tvedt's, Hohenbichler's, and Cai's methods are superior since they only include one probability of failure formula regardless of the sign of the curvatures, as opposed to Koyluoglu's, Mansour's, and Zhao's methods. Furthermore, these methods are also simpler than SOSPA since the latter requires the numerical solution of the saddle point.
In this study, we investigated the applicability of the existing methods for large-scale problems with less than 40 variables. For engineering systems which might involve a larger number of random variables, the studied SORMs can be used to estimate the component level reliability which in general has fewer random variables. Thereafter, system reliability analysis (Hao et al. 2020;Song and Der Kiureghian 2006;Radhika et al. 2008;Yu et al. 2018;Hu and Du 2018b;Cheng and Du 2016) for series and parallel systems can be used to compute the system level reliability. Furthermore, sensitivity analysis (Song et al. 2009;Hu and Mahadevan 2016b;Xiao et al. 2011;Jacques et al. 2006;Mara and Tarantola 2012;McRae et al. 1982;Sudret 2008;Hu and Xiong 2015) and dimension reduction techniques (Li et al. 2019;Hawchar et al. 2017;Hu and Youn 2011;Sargsyan et al. 2014;Guo et al. 2019;Yadav and Rahman 2012;Huang and Du 2006;Yadav and Rahman 2014;Bin et al. 2017) can also be adopted to reduce the complexity of large-scale engineering system before using SORM. To the best of the author's knowledge, the use of SORM for very large number of random variables, i.e., thousands, has been limited. Such high dimension typically occurs when spatial field uncertainties are included, in which each spatial location is represented by a random variable. For these problems, finding the MPP becomes inefficient since (1) the computational effort of computing finite-difference gradients during the MPP search increases linearly with the number of random variables, (2) a large number of iterations may be necessary to converge to the MPP in such high dimension, and (3) convergence to a global solution using conventional MPP-optimization algorithms becomes increasingly unlikely as the dimension increases, due to the potential presence of local minima.
Finally, to extend our study, the different SORMs should be applied to RBDO of large-scale engineering problems and their performance compared to advanced meta-modeling techniques combined with MC simulations. This should be performed for both time-independent and time-dependent reliabilities as well as system reliabilities to truly assess the applicability of the existing methods for engineering optimization under uncertainties.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.