A High-Dimensional Reliability Analysis Method for Simulation-Based Design Under Uncertainty

where the system safety indicator $IΩS$ equals 1 (0) given a safety (failure) state. Although the direct MCS can yield an accurate reliability estimate, achieving satisfactory accuracy often requires a large number of G function evaluations [2]. A popular alternative to the direct MCS is to use a surrogate model, $Ĝ(x)$⁠, which is built with a much smaller number of sample points, to approximate the original performance function, $G(x)$ [3]. A large number of surrogate modeling methods have been proposed in the literature, among which our review is mainly focused on two categories of these methods: high-dimensional model representation (HDMR)-based methods [4–19] and kriging-based methods [20–24].

The HDMR-based methods simplify the construction of the surrogate model by decomposing a high-dimensional performance function into several low-dimensional functions (or low-order component functions) [4,5]. In most practical applications, the response of an engineered system can be largely determined by the low-order interactions among the input variables, which makes HDMR an attractive option in these applications. HDMR was initially applied for efficient multivariate representations in chemical system modeling [6–8]. Numerical experiments showed that HDMR-based methods are promising for simulation-based design of high-dimensional systems [9]. Two of the primary and well-known HDMR-based methods are univariate dimension reduction (UDR) and bivariate dimension reduction (BDR), which get their names on the basis of whether the performance function is decomposed into first- or second-order component functions [10,11]. Each component function can be approximated with an interpolation technique, and a global surrogate model of the performance function can then be built based on the summation of the approximated component functions. The HDMR-based methods thereafter have been further developed for the application in an expanded range of areas (e.g., design optimization [12], metamodeling [13], and reliability analysis [14]). In design optimization, Li et al. [12] integrated HDMR with an expected improvement sampling strategy that compensates for the underestimation of the response uncertainty by the kriging predictor. The resulting HDMR-based method achieved satisfactory search performance in deterministic optimization. In metamodeling, Foo and Karniadakis [13] proposed the combined use of HDMR and stochastic collocation, where HDMR is first used to decompose a high-dimensional function into several low-dimensional functions, and multi-element probabilistic collocation is then used to approximate these low-dimensional functions. In reliability analysis, the random input variables could be statistically dependent and in the form of intervals. Xie et al. [14] proposed an efficient interval analysis algorithm based on the HDMR expansion to handle-dependent interval input variables in reliability analysis. Hu and Youn [15] proposed the adaptive dimension decomposition and reselection method that can automatically identify potentially important low-order component functions and adaptively reselect the truly important ones. Li et al. [16] coupled HDMR with an intelligent sampling strategy for building global surrogate models on nonlinear problems. To integrate HDMR with the sampling strategy, a projection-based intelligent method was introduced to locate the sample points along the corresponding cuts (lines, planes, or hyperplanes of high orders). This method treats all dimensions equally and is able to track the nonlinear regions in the input space. The component functions in HDMR, resulting from the decomposition of a performance function, can also be approximated by a family of linearly independent basis functions and then metamodeling or reliability analysis can be performed over the approximated component functions [17–19]. The basis functions, such as orthogonal polynomials and cubic splines, are defined as uni- or multivariate functions in a way that a proper combination of these basis functions can reflect the true behavior of the performance function. Hajikolaei and Wang [17] proposed the integration of principal component analysis with HDMR for finding the coefficients of orthogonal basis functions that provide the best linear combination of the bases with minimum error. Liu et al. [18] proposed the generalized radial basis function-based HDMR method, which employs the virtual regular points, projected from the random points, to update the predictions by the basis functions. The radial basis function-based HDMR method was further developed by using a recently proposed proportional sampling strategy that allows adding first-order sample points to efficiently construct each component basis function [19].

Given the negligibly weak high-order interactions, HDMR-based methods with low-order component functions would provide an efficient and accurate approximation of the performance function. However, several limitations of these HDMR-based methods do exist: (i) they usually consider only the low-order interactive effects among input variables; (ii) increasing the order of component functions may require a dramatic increase in the number of function evaluations (e.g., building the first-order component functions for a d-dimensional performance function only needs nd + 1 function evaluations with n + 1 equal to the number of sample points along each dimension; however, building the mth-order component functions with m ≥ 2 requires $dmnm$ function evaluations¹); and (iii) most of these methods build a surrogate of the performance function based on the function values determined at a set of structured and predefined grid points, thus lacking the flexibility in choosing the locations of the sample points.

Kriging-based methods first construct an initial surrogate model of the performance function with an initial set of observations (e.g., at a set of sample points generated by Latin Hypercube sampling), and then continuously refine the surrogate model by identifying and adding new sample points until adequate accuracy is achieved [20–29]. Two well-known kriging-based reliability analysis methods are efficient global reliability analysis and maximum confidence enhancement, both of which consider the prediction uncertainty and LSF proximity in selecting a new sample point [20–22]. A recently developed method, namely, sequential exploration–exploitation with dynamic trade-off (SEEDT), introduces an exploration–exploitation trade-off coefficient that adaptively weighs exploration and exploitation throughout sampling process to achieve a more optimal balance between these two criteria [23,24]. Two unique features that kriging-based methods possess are: (i) they are flexible in choosing the location of a new sample point, which enables these methods to sample in an optimized sequence; and (ii) they provide a probabilistic surrogate model that quantifies both the prediction of a performance function and its uncertainty. With these features, kriging-based methods have been successful in reliability analysis involving system performance functions that are computationally expensive to evaluate, and implicit and unknown to the user. Additionally, kriging-based surrogate modeling has also been successfully applied in real-time processing of multisensor signals for structural health monitoring [25]. However, these success stories have been on problems of low (typically < 10) dimension. When it comes to high-dimensional problems, the efficiency of these methods diminishes dramatically. It has been reported in the machine learning area that kriging does not scale to high-dimensional problems [30,31]. This severely limits the application of kriging to reliability analysis of systems that involve the use of complex simulation models with large dimensions (e.g., > 30 parameters) and expensive simulations (e.g., hours of run-time for numerically solving differential equations in one simulation).

In this paper, we tackle the above challenges by addressing efficient reliability analysis involving high-dimensional, computationally expensive, highly nonlinear performance functions. Due to the above mentioned limitations, neither HDMR- nor kriging-based methods are capable of providing an accurate and computationally efficient solution for these challenging problems. However, we find that appropriately modifying these two methods and optimally combining their strengths can lead to the development of a feasible solution that is capable of (i) capturing the strong interactions among variables; (ii) alleviating the curse of dimensionality; and (iii) flexibly choosing the locations of new sample points. Built upon this finding, this study proposes a new method for reliability analysis, named as high-dimensional reliability analysis (HDRA). The proposed HDRA method decomposes the task of surrogate model construction for reliability analysis into two sequential steps: (i) adaptive univariate sampling (with UDR and SEEDT) for building a global surrogate model; and (ii) kriging-based multivariate sampling (with SEEDT) to refine the global surrogate model in highly uncertain regions close to the LSF. At the first step, a newly proposed method, named as adaptive UDR (AUDR), adaptively locates significant univariate sample points by decomposing the important multivariate points identified by SEEDT. A UDR-based global surrogate model is then built based on these univariate sample points. At the second step, the global surrogate model is used as the trend function in the kriging model, and SEEDT further refines the surrogate model by sequentially locating critical multivariate sample points in highly uncertain regions close to the LSF.

The reminder of this paper is organized as follows: Section 2 gives a brief review on two related methods, UDR and SEEDT. Section 3 introduces the proposed HDRA method, a hybrid of AUDR and SEEDT, for reliability analysis. The effectiveness of the proposed method is demonstrated using four mathematical and engineering examples in Sec. 4. Section 5 provides several concluding remarks.

where $ϵ(x)$ is the truncation error of the univariate approximation, caused by dropping the interactive effects among variables on the response. Table 1 describes explicitly the UDR method in a pseudo-code. The algorithm starts by generating 3 or 5 sample points along each dimension (line 2). After evaluating the performance function at the sample points, each component function $Gi(xi)$ is approximated based on the cubic spline interpolation (line 5). At the end, the global surrogate model is built based on Eq. (4) (line 7).

The SEEDT method first builds an initial surrogate model with kriging and then sequentially refines the model. For the updating of the surrogate model, expected utility is used as the acquisition function to locate the next sample point. The sample points are sequentially chosen according to a new exploration–exploitation scheme that adaptively weighs exploring the region with high prediction uncertainty and exploiting the probable location of the LSF throughout the sampling process [23]. The exploration criterion ensures that the sequential sample points would not cluster in one region because the algorithm in SEEDT explores the whole input space for selecting the new sample points. In other words, if the algorithm places a number of sample points in one nonlinear area, the acquisition function will tend to guide the exploration to other nonlinear areas with higher prediction uncertainty. The sequential sampling procedure of this method is listed in Table 2. The algorithm starts by generating several initial Latin hypercube points using the Latin hypercube sampling technique. The initial surrogate model is built with kriging upon the initial sample points (line 1). At each iteration, a new sample point is located by maximizing the expected utility function EU (lines 3–6). The surrogate model is then updated with the augmented data set (lines 7, 8).

In both kernel functions, $θ$ is the hyper-parameter vector that determines the smoothness of the prediction. The maximum likelihood estimation technique is performed to estimate these hyper-parameters. Figure 1 schematically illustrates the correlation between two arbitrary points $xandx′,$using both the squared exponential kernel (k_se) and additive kernel (k_add) functions with $θ=1,1T$ in a two-dimensional (2D) space. It can be observed that k_se is related to the effect of the direct bivariate distance, l, between $xandx′$⁠, while k_add corresponds to the sum of the effects of the decomposed univariate distances, l₁ and l₂, between these two points.

where $y=y1,…,ytT$ is a vector of t responses, $M$ is the trend function value at sample points {$x1,…,xt$}, and $r(x)=ψ(x,x1),…,ψ(x,xt)T$ is a correlation vector [26].

The overall flowchart of the proposed HDRA method is shown in Fig. 2. The method is composed of two sequential steps. At the first step, AUDR adaptively identifies significant univariate sample points by first identifying important multivariate sample points with SEEDT and then decomposing the multivariate points into univariate points along all active dimensions. This step builds a global surrogate model with Eq. (4) while ensuing satisfactory approximation accuracy of the univariate component function along each dimension. Upon the completion of adaptive univariate sampling in AUDR, the global surrogate model is used in step 2 as the trend function in the kriging model, and SEEDT further refines the surrogate model by sequentially locating critical multivariate sample points in regions near the LSF and with high prediction uncertainty. Table 3 describes explicitly the procedure of HDRA. Steps 1 and 2 are described in lines 1–14 and 15–24 of the pseudo-code, respectively. In what follows, the two main steps are explained in further detail.

It should be mentioned that the first part of the algorithm (lines 1–14) does not consider the interactive effects among variables, while the second part (lines 15–24) efficiently captures the interactions among variables by identifying multivariate sample points in regions near the LSF and with high uncertainty in prediction. We use the additive kernel function, Eq. (9), in the correlation matrix $Ψ$ for determining the correlation between sample points.

It can be observed from the LSF that the performance function is more highly nonlinear along the $x1$ dimension. UDR treats both dimensions to be of equal importance, while AUDR identifies the higher degree of nonlinearity along the $x1$ dimension and adaptively assigns more univariate points to the $x1$ dimension and fewer points to the $x2$ dimension. Furthermore, AUDR allocates more points to regions that potentially contribute more to improving the accuracy in the LSF and reliability predictions.

AUDR does not incorporate any multivariate sample points during the sampling process, and thus the resulting global surrogate of the performance function lacks the ability to capture the interactive effects among variables. To address this limitation, after step 1 is completed, the SEEDT method is utilized again to adaptively refine the surrogate model by sequentially locating multivariate sample points in regions close to the LSF and with high prediction uncertainty. For the refinement of the global surrogate, multivariate sample points are chosen according to a sequential exploration–exploitation scheme in SEEDT that dynamically weighs exploring the regions with high prediction uncertainty and exploiting the ones close to the LSF. The procedure of SEEDT (see Algorithm 2 in Sec. 2.2) is used here as step 2 of the HDRA method, and the surrogate model built in step 1, $ĜAUDRx,$ is employed as the trend function of the kriging model in Eq. (6). To allow for the consideration of the interactions among variables, this step uses the squared exponential kernel in Eq. (8) for building the correlation matrix. After building the initial kriging model, a new sample point $xk$ is chosen by maximizing the expected utility function (line 18). This new point is then added to the data set, and the kriging model is updated with the augmented data set (lines 19–21). This procedure is repeated until the target accuracy (see the convergence estimator in line 22) is achieved or the maximum number of iterations is reached (line 17). Interested readers are referred to Ref. [23] for more information about the convergence estimator in SEEDT.

As described earlier, HDRA consists of two sequential steps: (i) adaptive univariate sampling (with AUDR) for building a global surrogate that captures the global trend of a performance function (see Sec. 3.1) and (ii) kriging-based multivariate sampling (with SEEDT) to refine the global surrogate for capturing the variate interactions in the performance function (see Sec. 3.2). AUDR in step 1 does not add any multivariate samples to the sample data set and thus the resulting global surrogate of the performance function lacks the ability to capture the interactive effects among variables. To address this limitation of the global surrogate, SEEDT in step 2 sequentially identifies multivariate samples in regions close to the LSF and with high prediction uncertainty. The addition of these multivariate samples produces a refined surrogate that is capable of capturing the variate interactions in the performance function.

It is important to note that the initial samples (i.e., the univariate samples added in step 1) greatly influence the accuracy and efficiency of the proposed HDRA method. In what follows, we will elaborate this influence based on two extreme scenarios:

• Scenario 1—considering too many initial samples: The addition of more initial samples in step 1 (AUDR) would result in higher accuracy in the approximation of the first-order component functions in the UDR model of the performance function (see Eq. (4)). This UDR model, regardless of how many univariate samples are incorporated, cannot capture the interactive effects among variables. If the performance function contains strong variate interactions, the kriging-based multivariate sampling in step 2 (SEEDT) is always essential, as it overcomes the limitation of the UDR approximation in step 1. In this regard, considering too many initial univariate samples in step 1 could exhaust the computing budget, thereby allowing for the addition of only few multivariate samples in step 2 and producing an inaccurate surrogate model that does not adequately represent the variate interactions.

• Scenario 2—considering too few initial samples: The smallest number of initial samples in step 1 is 2d + 1, where d is the number of random variables in the performance function. Assuming only 2d + 1 univariate samples are considered under this scenario, step 1 produces an initial global surrogate model built with UDR without adaptive enrichment of univariate samples along any dimension. This initial UDR model may fail to capture the global trend of the performance function, especially along dimensions with high degrees of nonlinearity. In such cases, no matter how many multivariate samples are added in step 2 to refine the global surrogate, the surrogate may always fail to capture both the high dimensionality and strong variate interactions. Under this scenario, the final HDRA model would likely be a slightly improved version of ordinary kriging with the UDR model as the trend function.

To conclude, adaptive univariate sampling (AUDR) is used in step 1 to efficiently build a global trend function for dealing with high dimensionality, and kriging-based multivariate sampling (SEEDT) is adopted in step 2 to efficiently refine the global trend function for capturing the interactive effect among variables. It is important to determine a proper number of initial samples in step 1 and this sample size determination is automatically performed in AUDR by stopping the adaptive univariate sampling when the information gain in all the dimensions by adding a new sample is lower than a predefined threshold (see Eqs. (12) and (13)).

Three mathematical examples and one real-world problem are used to evaluate the performance of the proposed HDRA method. The purpose of each example is summarized in Table 4. Example 1 compares the reliability errors by four methods (UDR, AUDR, SEEDT, and HDRA) for different degrees of response nonlinearity. Example 2 employs a highly nonlinear performance function and investigates the effect of the reliability level on the efficiency of UDR, AUDR, SEEDT, and HDRA. Example 3 is also a mathematical example chosen to evaluate the effect of dimensionality on the performance of the three methods. Then, Example 4 estimates the practicality of the proposed HDRA method with a real-world application. In this example, the performance function is defined as the power generation of a piezoelectric energy harvester and can be represented by a function of six intrinsic input variables (i.e., three geometric design variables and three material property variables), embedded in a space of extrinsic dimensionality d = 15. That is, we add nine extrinsic input variables, each of which contains uncertainty but does not have any effect on the performance function.

In these four examples, the reliability estimate by each method is compared with that by the direct MCS to evaluate the accuracy of the method. If one wants to directly estimate the accuracy of the surrogate model produced by one method, the LSF error estimator defined in Ref. [23] can be used that may address the potential issue that two methods give very similar reliability estimates but produce very different surrogate models and thus estimates of the LSF. In the examples of this study, the reliability estimation results by the proposed method have consistently shown better accuracy than those by the benchmark methods. As a result, only the reliability estimation accuracy is used to compare the performance of these methods.

Due to the randomness in the initial sample selection and MCS (for reliability analysis), the reliability estimation errors by UDR, AUDR, SEEDT, and HDRA contain uncertainty. To capture this uncertainty, we repeatedly run each method with the same parameter setting for ten times for all the examples. The performance metric of each method is presented by the mean and uncertainty (⁠$±σ$⁠) of the metric over the ten repeated runs.

where the coefficient b is used to adjust the nonlinearity of the performance function along the first direction. All 40 input variables are independent and normally distributed (mean μ = 1.5, and standard deviation σ = 1). The corresponding reliability is defined as R = P(G(x) < 0), which is kept in the same level for all chosen values of b. MCS is performed to set the benchmark results of reliability analysis for each level of nonlinearity (sample size: 1 million). The relative reliability error, which is the ratio of the absolute difference between the reliability estimate by each method and that by the direct MCS to the reliability estimate by the direct MCS, serves as the criterion for evaluating the performance of each method. The relative reliability errors by UDR, AUDR, SEEDT, and HDRA for different nonlinearity levels are graphically compared in Fig. 4. The numbers of sample points (or function evaluations) required by UDR, SEEDT, and HDRA are 161, 95, and 92, respectively. The number of sample points for AUDR varies from 83 to 86. From Fig. 4, HDRA, in general, is more accurate than UDR, with considerably fewer sample points, and more accurate than AUDR and SEEDT, with approximately the same number of sample points. SEEDT, as a kriging-based reliability analysis method, fails to produce satisfactory accuracy (i.e., <1% reliability error) in this high-dimensional problem, and this suggests that the limitation of kriging in high-dimensional problems reported in the machine learning area also applies to reliability analysis.

where $ωp=kp/mp$⁠, $ωs=ks/ms$⁠, $ωa=(ωp+ωs)/2$⁠, $ξa=(ξp+ξs)/2$⁠, $γ=ms/mp$ and $θ=(ωp−ωs)/ωa$⁠, x=[$ωp,ωs,ωa,ξa,γ,θ$]. Table 5 summarizes the statistical information of the input variables.

With the reliability estimate by MCS as the benchmark, changing the mean value of $Fs$ from 15.0 N to 27.5 N results in a reliability interval between 98.22% and 99.77%. Figure 5 graphically compares the relative reliability errors of UDR, AUDR, SEEDT, and HDRA at seven reliability levels within the interval. Table 6 summarizes the numbers of function evaluations, reliability estimates, and relative errors by the three methods. As can be seen in the figure and table, the proposed HDRA method yields the minimum error among the four methods with the fewest function evaluations. It is important to emphasize that although UDR and AUDR succeed in accurately approximating all first-order component functions by considering only a small number of sample points at each dimension, they fail to produce an accurate prediction of the reliability. This indicates that the performance function in Eq. (18) may not show high nonlinearity along each dimension, but it is most affected by high-order interactions among the random variables. HDRA captures these interactive effects by adding to the sample data set multivariate sample points in regions close to the LSF and with high prediction uncertainty.

The above definition allows us to create performance functions of varying dimensions by simply changing the value of input dimension d. Figure 6 compares the reliability errors of UDR, SEEDT, and HDRA under four settings of input dimension d. To make the comparison fair, we keep the numbers of function evaluations the same (i.e., 20d + 1) for the three methods. The number of initial sample points for HDRA is set as 2d + 1, and in most cases, AUDR automatically considers three more sample points along each direction to refine the univariate component functions. As can be seen in Fig. 6, HDRA yields the lowest error value among the three methods when d = 6, 10, and 14. SEEDT produces slightly better accuracy than HDRA when d = 2, and somewhat comparable accuracy when d = 6. These results suggest that kriging-based methods have a strong capability to handle low-dimensional reliability analysis problems. However, SEEDT fails to accurately predict the reliability for the cases of d = 10 and d = 14. Furthermore, the relative performance of UDR, as compared to SEEDT, improves with the input dimension, d, which shows the unique strength of UDR in handling high-dimensional problems. By taking the advantages of the unique strengths of SEEDT and UDR for low- and high-dimensional problems, HDRA is able to achieve satisfactory performance in all cases.

Vibration energy harvesting is widely used to transform commonly wasted energy of vibration into accessible energy, which can then be applied to charge supercapacitors, batteries, or enable self-powering sensors. The typical configuration of an energy harvester is shown in Fig. 7. It consists of a shim laminated by piezoelectric materials at one end and a tip mass attached at the other end. The tip and shim mass are constructed from tungsten/nickel alloy and Blue steel, respectively [36].

Under the piezoelectric effect, mechanical strain is transferred into electric voltage or current. This study considers 31 modes, which allows for higher longitudinal strain when energy harvester is subject to smaller input forces. Under longitudinal stress/strain, voltage is produced along the thickness direction. The piezoelectric harvester plays similar role as a transformer circuit. Per the Kirchhoff's voltage Law, the circuit can be expressed by the coupled differential equations that describe the conversion from mechanical stress/strain to electrical voltage. The conversion process can be simulated by matlabsimulink. The harvester output is determined by the geometries of input and material properties. A detailed description of the energy harvester model can be found in the authors' previous publication [36]. The study in Ref. [36] aims to optimize the design of the energy harvester with reliability-based design optimization (RBDO). Three geometric terms, $x=lb,lm,hm$⁠, and three material properties, $υ=υ1,υ2,υ3$⁠, were considered as the random design variables and the random parameters, respectively. For all these random terms, the means and standard deviations were considered fixed. In this RBDO problem, the size of the energy harvester was minimized while fulfilling the reliability requirement with respect to power generation.

Besides $x$ and $υ$⁠, we add nine random input variables (⁠$z1,…,z9$⁠) that do not affect the power output of the energy harvester, making the performance function of high extrinsic dimension (i.e., d = 15). This treatment is to demonstrate the capability of HDRA in dealing with high-dimensional reliability analysis problems. The distributional information of the 15 random variables is summarized in Table 7.

The reliability level is estimated with the UDR, SEEDT, and HDRA method with a total of 121, 120 and 120 sample points, respectively. In step 1 of HDRA, 37 initial sample points are selected to build the global trend function. In HDRA, it is observed that in step 1 of the algorithm, all nine ineffective dimensions are considered as inactive dimensions. Therefore, no more sample points are selected along these directions. This feature helps the algorithm focus more on the effective variables (i.e., the first six variables) by locating more sample points along them. The direct MCS serves as the benchmark method for the reliability estimation. Since the uncertainty of mean reliability estimates of the direct MCS with 10,000 sample points is much smaller than that of the HDRA method and the process of function evaluation is computationally expensive, it is concluded that 10,000 sample points are enough for the benchmark estimation of reliability. Figure 8 graphically summarizes the reliability estimates by UDRA, SEEDT, HDRA, and MCS for different levels of P₀. It can be observed that the reliability estimates by HDRA are closer to the benchmark reliability levels acquired by the direct MCS in comparison with those by UDR and SEEDT. The reliability analysis results of this real-world engineering application further demonstrate the effectiveness and accuracy of the proposed method in high-dimensional reliability analysis.

Table 8 lists the main features of four reliability analysis methods, UDR, AUDR, SEEDT, and HDRA. These features are identified by analyzing the intrinsic structures of the algorithms in the methods and interpreting the results from the above examples. UDR is capable of alleviating the curse of dimensionality, and AUDR further enhances this capability by incorporating an adaptive sampling strategy. Both methods are particularly useful for high-dimensional problems with weak variate interactions. However, neither possesses the capability to handle performance functions with strong variate interactions. Although this limitation can be alleviated by considering higher order component functions (e.g., in BDR and trivariate DR), the computational efficiency quickly diminishes with the increase of the input dimension. Kriging-based methods such as SEEDT can capture strong variate interactions and are favored in low-dimensional problems (typically d < 6) with strong variate interactions. But similar to the limitation of BDR and trivariate DR, these methods lose efficiency in high-dimensional problems. HDRA optimally combines the strengths of UDR and SEEDT, and can achieve satisfactory accuracy and efficiency for problems of both ranges of dimension and with various degrees of variate interactions. The HDRA method provides a better alternative to these popular existing methods by alleviating the fundamental limitation of UDR in tackling strong variate interactions and high computational cost of SEEDT in high-dimensional problems.

The high-dimensional reliability analysis (HDRA) method has been proposed to solve reliability analysis problems involving high dimensionality, computationally expensive simulations, high nonlinearity, and strong variate interactions. The basic idea of HDRA is to first build a global surrogate model with AUDR and then locally refine the surrogate model using SEEDT. AUDR adaptively identifies important univariate sample points considering both response nonlinearity and criticality (with respect to accurate reliability prediction), while SEEDT captures the interactive effects among variables on the performance function by adding multivariate sample points in highly uncertain regions close to the LSF. Results from four mathematical and real-world examples with up to 40 dimensions suggest that the HDRA method can achieve significantly higher efficiency in reliability analysis than the existing DR- and kriging-based methods, and is especially useful for high-dimensional, computationally expensive problems with strong interactions among random variables. Future research will further develop the proposed method for a wider range of applications, such as sequential experimental design for high-dimensional uncertainty quantification and simulation-based design under high-dimensional uncertainty.

This research was in part supported by the U.S. National Science Foundation (NSF) and the NSF I/UCRC Center for e-Design. Any opinions, findings or conclusions in this paper are those of the authors and do not necessarily reflect the views of the sponsoring agency.

• National Science Foundation (Grant Nos. CNS-1566579 and ECCS-1611333).