An efficient surrogate-aided importance sampling framework for reliability analysis
AAn efficient surrogate-aided importance sampling framework for reliability analysis
Wang-Sheng LIU , Wen-Jun CAO , and Sai Hung CHEUNG School of Civil & Environmental Engineering, Nanyang Technological University, Block N1, 50 Nanyang Avenue, Singapore 639798 Department of Civil and Environmental Engineering, National University of Singapore, Block E1A, No.1 Engineering Drive 2, Singapore 117576 ETH Zurich, Future Cities Laboratory, Singapore-ETH Centre , Singapore 138602
Abstract:
Surrogates in lieu of expensive-to-evaluate performance functions can accelerate the reliability analysis greatly. This paper proposes a new two-stage framework for surrogate-aided reliability analysis named Surrogates for Importance Sampling (S4IS). In the first stage, a coarse surrogate is built to gain the information about failure regions; the second stage zooms into the important regions and improves the accuracy of the failure probability estimator by adaptively selecting support points therein. The learning functions are proposed to guide the selection of support points such that the exploration and exploitation can be dynamically balanced. As a generic framework, S4IS has the potential to incorporate different types of surrogates (Gaussian Processes, Support Vector Machines, Neural Network, etc.). The effectiveness and efficiency of S4IS is validated by five illustrative examples, which involve system reliability, highly nonlinear limit-state function, small failure probability and moderately high dimensionality. The implementation of S4IS is made available to download at (the link will be added once the paper is accepted for publication).
Keywords: reliability analysis, stochastic sampling, metamodel, active learning, design of experiment
Introduction (Background: reliability analysis) Reliability analysis is very essential for ensuring a system to perform as designed when subjected to different sources of uncertainties. In the context of engineering, the measure of the reliability is the failure probability. Given a model whose uncertain parameters are represented by a vector of ๐ random variables ๐ฏ = [ฮ & , โฏ , ฮ ) ] + , the failure probability ๐ - for a specific failure event ๐น can be written as: ๐ - = / ๐ - (๐)๐ (๐)d๐ = ๐ผ [๐ - (๐)] (1) where ๐ = [๐ & , โฏ , ๐ ) ] + โ ๐บ โ โ ) denotes the parameter value of ๐ฏ in the parameter space ๐บ ; ๐ - (โ) is the binary indicator function of the failure event ๐น which equals to one if ๐ belongs to the failure region ๐บ - and zero otherwise; ๐ (โ) is the joint probability density function (PDF) of ๐ฏ ; and ๐ผ [โ] denotes the expectation for ๐ฏ~๐ (๐) . The failure region ๐บ - is defined by the performance function ๐(๐) , i.e., ๐บ - = {๐: ๐(๐) โค 0} . The corresponding limit-state function ๐(๐) = 0 divides the parameter space ๐บ into the failure and safety regions. (Reliability methods) A considerable amount of literature has been focusing on developing efficient and accurate ways to calculate the failure probability. One widely used set of reliability methods is based on the linear or quadratic approximation of the limit-state function around the most probable point (MPP) using Taylor series approximation, usually known as the first or second order reliability methods (FORM or SORM), respectively [1, 2]. But the error brought by the approximation tends to be large when the limit-state function is highly nonlinear and/or there exists multiple MPPs in the failure region; furthermore, the computational cost in the search of MPP increases drastically with the number of random variables. Another set of reliability methods is based on stochastic sampling. Among them, Monte Carlo Simulation (MCS) has been acknowledged as a robust technique where the failure probability is approximated statistically through stochastic sampling. In its crude form, MCS first simulates ๐ samples from the PDF ๐ (๐) ; then evaluates the corresponding values of the performance function ๐(๐) and counts the number of failure samples ๐ - ; finally, an estimator ๐I - can be given by ๐ - /๐ . Though ๐I - is asymptotically unbiased, the convergence rate of the crude MCS is low. For example, to estimate a small failure probability of LM , it takes a very large ize of MNO samples for the estimator to achieve a Coefficient of Variation (
CoV ) close to 10%. Advanced stochastic simulation techniques have been developed to reduce the required sample size, including Subset Simulation (SS) [3], Importance Sampling (IS) [4], and Spherical Subset Simulation (S ) [5]. Despite the remarkable improvement in efficiency, they still require a moderately large number of samples to be evaluated (e.g., in SS, the sample size in each level is typically set to 1000). (Motivation for surrogates) In practice, the performance function ๐(๐) of a complex system often involves implicit and time-demanding numerical models such as finite element models which are widely used by engineers to predict the responses of a system. Whatโs worse, in the attempt to better mimic the behavior of a system, numerical models tend to be more and more sophisticated, making the model evaluation increasingly computationally expensive. Consequently, the aforementioned reliability methods are likely to be inapplicable for complex systems when the computation resource is limited. (Review of surrogates and learning strategies) To alleviate this problem, surrogates (also referred to as metamodels or response surfaces) can be built in lieu of the original expensive performance functions. In the literature, different surrogates have been employed to perform reliability analysis, such as polynomials [6, 7], Gaussian process (GP or Kriging) [8-12], polynomial chaos expansion (PCE) [13, 14], artificially neural networks (ANN) [15, 16] and support vector machines (SVM) [17]. The most critical part of efficiently building a surrogate is probably the adaptive infill strategy, also referred to as adaptive design of experiment or refinement scheme in some literature. Instead of densely filling the parameter space with support points, an adaptive infill strategy often starts with a coarse surrogate built upon a small number of support points, and then refines it progressively using judiciously selected support points. Many GP-based adaptive infill strategies have been reported in the last decade in conjunction with stochastic-sampling-based reliability methods. For example, Bichon et al. [9] presented an Efficient Global Reliability Analysis (EGRA) method where an Expected Feasibility Function (EFF) was maximized by using the global optimization algorithm to select the best next infill support points in each iteration. More recently, Echard et al. [18] proposed a new method AK-MCS which combines GP with crude Monte Carlo Simulation. In this method, an active learning function called U function was developed to guide the selection of equential support points. Later, a modified version AK-IS [12] was proposed using the same learning function, which resorts to importance sampling for solving small-failure-probability cases. Alternatively, Dubourg et al. [11] adaptively built a GP surrogate as a substitute of the performance function and used it to construct the quasi optimal importance sampling density. Surrogate-aided reliability analysis is very parsimonious with the help of progressive refinement, but most adaptive infill strategies in the current literature are customized for GP and specific reliability analysis problems only. The paper presents a new framework for surrogate-aided reliability analysis named Surrogates for Importance Sampling (S4IS). This framework is robust to different types of reliability analysis problems (e.g., high dimensionality and small failure probability) and has the potential to incorporate different types of surrogates. The estimation of the failure probability and the building of the surrogate are performed in two stages. In Stage 1, support points are generated to build a coarse surrogate which can identify multiple failure regions; the next stage focuses more on the important regions to improve the failure probability estimator. The adaptive infill of support points in two stages is under the guidance of two different learning functions such that the dynamic balance between the exploration and exploitation can be achieved. The proposed learning functions are the convex combination of several objectives. The remainder of this paper is organized as follows. Section 2 introduces the importance sampling technique for reliability analysis as well as the existing learning functions used to guide the infill of support points. In Section 3, the proposed two-stage framework S4IS is discussed in detail. To validate the proposed framework, five illustrative examples are investigated in Section 4, which involves system reliability problem, highly nonlinear limit-state function, different reliability levels and moderately high dimensionality. The performance of S4IS is compared with the crude MCS, FORM and AK-IS. Surrogate-aided importance sampling
Importance sampling
It is common that the system failure ๐น in Equation (1) is a rare event, meaning the failure region ๐บ - is in the tail part of ๐ (๐) . In this case, the crude MCS may cannot afford to imulate sufficient failure samples in order to achieve a small variance of the estimated ๐ - . Alternatively, Importance Sampling (IS) [4] is widely used as a variance reduction technique which computes the expectation with respect to a proposed instrumental PDF. Given an instrumental PDF ๐(๐) supported on ๐บ โ โ ) , the failure probability in Equation (1) can now be rewritten as: ๐ -,QR = / ๐ - (๐) ๐ (๐)๐(๐) ๐(๐)d๐ = ๐ผ S T๐ - (๐)๐ (๐)๐(๐) U (2) Instead of sampling from the original PDF ๐ (๐) , we can sample from ๐(๐) and adjust ๐ - (๐) with the likelihood ratio ๐ (๐)/ ๐(๐) . Now the IS estimator ๐I -,QR can be expressed as: ๐I -,QR = 1๐ QR V ๐ - W๐ (X) Y ๐ W๐ (X) Y๐(๐ (X) ) Z [\ X]& , ๐ (X) ~๐(๐) (3) where ^๐ (X) ; ๐: 1 โ ๐ QR b denotes a set of ๐ QR i.i.d (independent and identical distributed) samples from ๐(๐) . The IS estimator is unbiased and its variance is estimated by: ๐๐๐f๐I -,QR g = 1๐ QR (๐ QR โ 1) V T๐ - W๐ (X) Y๐ W๐ (X) Y๐(๐ (X) ) โ ๐I -,QR U OZ [\ X]& = 1๐ QR โ 1 i 1๐ QR V ๐ - W๐ (X) Y๐ W๐ (X) Y O ๐(๐ (X) ) OZ [\ X]& โ ๐I -,QRO j (4) The selection of the instrumental PDF is of crucial importance in reducing the variance. It is not difficult to derive that minimizing the variance gives the optimal instrumental PDF ๐ โ (๐) = ๐ - (๐)๐ (๐)โซ ๐ - (๐)๐ (๐)d๐ = ๐ - (๐)๐ (๐)๐ - (5) Although this zero-variance optimal density in Equation (5) is not usable without knowing ๐ - , it can guide the selection procedure of ๐(๐) and provide insights to the design of the importance sampling scheme. In practice, quasi-optimal instrumental PDFs are commonly used. Building the surrogate adaptively (Surrogate for what) Assume the performance function of interest ๐(๐) is a time-demanding black-box function. The main goal is to build a surrogate ๐m(๐) to the original function from a dataset of support points
๐ฎ = ^W๐ (X) , ๐ฆ (X)
Y; ๐: 1 โ ๐ p b so that all failure regions ๐บ - over the parameter space ๐บ can be identified accurately. In other words, in contrast with general-purpose surrogates that are optimized to be accurate over the whole parameter space, surrogates or reliability analysis focus on the accuracy of failure region boundaries โ๐บ - defined by the limit-state function ๐(๐) = 0 . Given the surrogate ๐m(๐) , the binary indicator function in Equation (2) and (3) can be approximated by: ๐r - (๐) = s1, ๐m(๐) โค 00, ๐m(๐) > 0 (6) (Active and passive infill strategy) An infill strategy should be employed to select a finite number of support points aiming at gaining the maximal information for building the surrogate. Instead of using dense support points for fitting a general-purpose surrogate, an adaptive infill strategy starts with a small size of support points to fit an initial surrogate and then new support points are sequentially selected to refine the previous surrogate. This way is proved to be far more efficient than the non-adaptive one since it is customized to perform the analysis of interest (e.g., reliability analysis in this paper). (How to build adaptively) Finding the best next support point(s) in each iteration is guided by some criteria or objectives in terms of optimization. In general, on one hand, we expect the new support point to explore the under-explored regions where existing support points are sparse; one the other hand, the previously explored regions should be exploited to improve the accuracy over the regions of interest. These two criteria are referred to as exploration and exploitation , respectively. In the context of reliability analysis, exploration involves adding more support points to sparsely populated regions in order to identify the failure region(s) and exploration involves zooming into the regions that contribute most to the integral for estimating the failure probability. As exploration and exploitation are conflicting goals, the problem of finding the next support point(s) is essentially a multi-objective optimization problem. Various Learning functions have been proposed as objectives to solve the multi-objective optimization problem. An early learning function called Expected Feasibility Function (EFF) was introduced by Bichon et al. [9] in the Efficient Global Reliability Analysis (EGRA). EFF converts a multi-objective problem to a single objective problem, acting as an indicator of how well the true value is expected to satisfy the limit-state function ๐(๐) = 0 in its vicinity of ยฑ๐ ๐ธ๐น๐น(๐) = / (๐(๐) โ |๐m(๐)|)๐ zI d๐m(๐) N{L{ (7) here ๐m(๐) is a realization of a Gaussian Process (GP) distributed as ๐ zI ; and ๐ is set to be proportional to the standard deviation of ๐m(๐) (e.g., ๐ = 2๐ ~m in [9]). The EFF in Equation (7) can be expressed in the analytical form when random variables ๐ฏ are transformed into the uncorrelated standard normal space via isoprobabilistic transformations [19-21]. The support point with the maximum expected feasibility is found by a global optimization method. It can be observed from Equation (7) that the optimization of EFF favors points that are close to the limit-state boundary (small |๐(๐)| ) and have large variance in the prediction (large ๐ ~m ) Similarly, another learning function in terms of the prediction mean and variance of a GP surrogate, called U-function, was introduced in AK-MCS [18] and AK-IS [12] ๐(๐) = |๐ zI (๐)|๐ zI (๐) (8) where ๐ zI (๐) and ๐ zI (๐) denote the prediction mean and standard deviation function of GP, respectively. The next point(s) can be obtained by minimizing ๐(๐) . A sampling-based approach is used to replace the global optimization step in EGRA, where the values of
๐(๐) are calculated for a large number of candidate support points drawn from the original distribution of random variables. The above two learning functions interpreted exploitation as the closeness to limit-state boundary but neglected the fact that the high-density regions in the failure domain contribute more to the failure probability than the low-density ones. To address this issue, Sun et al. [22] combined the statistical information provided by the GP surrogate with the joint PDF of random variables by proposing a new learning function named Least Expected Improvement Function (LIF). The proposed framework
In the existing framework described in [12] and [23], the first stage consists in the search of the most probable point (MPP) and reliability analysis using FORM; in the second stage, IS is performed by proposing a multivariate Gaussian distribution as the instrumental PDF, which is centered on the MPP and with covariance being the identity matrix. However, this framework can only incorporate Gaussian Process (GP) as the surrogate; furthermore, it is not able to identify multiple MPPs or multiple failure regions. n this paper, a new two-stage framework for surrogate-aided importance sampling called Surrogates for Importance Sampling (S4IS) is presented as illustrated in Figure 1. In the first stage, a coarse surrogate is built adaptively to identify all failure regions in the parameter space. The second stage starts with the construction of the Gaussian Mixture (GM) instrument PDF for IS based on the identified failure regions, then proceeds to the refinement of the surrogate mainly over the important regions to reduce the error of the IS estimator. Two different learning functions proposed in this framework are formulated by the convex combination of different exploration and exploitation objectives and are applicable to different types of surrogates. Also, the learning functions in the two stages balance the exploration and exploitation dynamically, with the one in the second stage focusing more on exploitation to ensure the fast convergence. As shown in Figure 1, infill support points in Stage 1 scatter more uniformly over the whole parameter space than those in Stage 2.
Figure 1: Schematic diagram of the proposed two-stage framework Surrogate for Importance sampling (S4IS)
Stage 1 3.1.1
Initialization
Throughout this paper, an isoprobabilistic transformation ๐ [19-21] is applied to the original random variables ๐ฏ to transform them into independent standard normal random variables ๐ , i.e., ๐ = ๐(๐ฏ) . The IS estimator in Equation (3) can be rewritten as: ๐I -,QR = 1๐ QR V ๐ - W๐ฎ (X) Y ๐ (cid:132) W๐ฎ (X) Y๐(๐ฎ (X) ) Z [\ X]& , ๐ฎ (X) ~๐(๐ฎ) (9) . ๐ " ๐ . Stage 1 ๐บ % Initial support pointsSequential infill support points
Stage 2 ๐๐บ % ๐บ % here ๐ (cid:132) (โ) is the multivariate standard normal PDF and the failure region in ๐ -space can be given by ๐บ - = ^๐ฎ: ๐W๐ L& (๐ฎ)Y โค 0b . To explore the parameter space, a large number ๐ (cid:133)& (in this paper, ๐ (cid:133)& = minW10 (cid:137) , 10 ) Y ) of samples are randomly generated from a distribution ๐ & (๐ฎ) in the first stage. They are treated as the inputs of candidate support points ๐ฎ (cid:133)& = ^๐ฎ (X) ~๐ & (๐ฎ); ๐: 1 โ ๐ (cid:133)& b . After selecting a small subset ๐ฎ p&,5 from ๐ฎ (cid:133)& , the outputs of the subset are evaluated ๐ฆ p&,5 = (cid:138)๐ฆ (X) = ๐ (cid:139)๐ L& W๐ฎ (X) Y(cid:140) ; ๐ฎ (X) โ ๐ฎ p&,5 (cid:141) (10) where ๐ฆ p&,5 denote the outputs of the time-demanding performance function for ๐ฎ p&,5 . Given the dataset of initial support points ๐ฎ &,5 = ^W๐ p&,5 , ๐ฆ p&,5 Yb = ^W๐ L& W๐ฎ p&,5 Y, ๐ฆ p&,5 Yb , the initial surrogate ๐m &,5 can be built to predict the output ๐ฆm = ๐m &,5 (๐) = ๐m &,5 W๐ L& (๐ฎ)Y (11) In this paper, we suggest the size of ๐ฎ p&,5 to be ๐ p&,5 = max[12, (๐ + 1)(๐ + 2)/2] where ๐ is the dimension of ๐ฎ . The inputs of candidate support points ๐ฎ (cid:133)& in the first stage can be simulated from the standard normal distribution, i.e., ๐ & (๐ฎ) = ๐ (cid:132) (๐ฎ) . However, for cases of small failure probabilities, only a few or even not a single ๐ฎ (cid:133)& will scatter in the failure regions, leading to at least two problems: (i) since the inputs of support points is selected from the candidates, the surrogate fitted using the corresponding support points will hardly be accurate near the limit-state boundary; (ii) multiple failure regions cannot be effectively identified. In this paper, ๐ฎ (cid:133)& are generated uniformly in the hypercubic space ranging from -5 to 5 in each dimension so as to explore the whole parameter space, i.e., ๐ & (๐ฎ) = &&5 (cid:145) , ๐ฎ โ [โ5,5] ) . Adaptive refinement of the surrogate
At the beginning of iteration ๐ (๐ = 1,2 โฏ ) , the current surrogate ๐m &,(cid:148) has been fitted using the dataset of support points ๐ฎ &,(cid:148) = ^W๐ p&,(cid:148) , ๐ฆ p&,(cid:148) Yb = ^W๐ L& W๐ฎ p&,(cid:148) Y, ๐ฆ p&,(cid:148) Yb (for ๐ = 1 , ๐m &,(cid:148) = ๐m &,5 and ๐ฎ &,(cid:148) = ๐ฎ &,5 ). Our aim is to find the best next input of support point from the candidates ๐ฎ (cid:133)& \๐ฎ p&,(cid:148) . The existing inputs of support points ๐ฎ p&,(cid:148) is excluded from the candidate set ๐ฎ (cid:133)& to avoid duplication. he learning function is essential to the selection of the location of new support points. For the first stage, the learning function achieves the balance between exploration and exploitation by minimizing the distance to the limit-state boundary (exploration) and meanwhile maximizing the distance to the existing support points (exploitation) ๐ฟ๐น & (๐ฎ) = (cid:151)๐m &,(cid:148) W๐ L& (๐ฎ)Y(cid:151) โ (cid:152)๐ฎ, ๐ฎ p&,(cid:148) (cid:152) (cid:153)X(cid:132) (12) where (cid:152)๐ฎ, ๐ฎ p&,(cid:148) (cid:152) (cid:153)X(cid:132) is the minimum distance between ๐ฎ and ๐ฎ p&,(cid:148) in iteration ๐ . Therefore, the new input of support point selected from ๐ฎ (cid:133)& \๐ฎ p&,(cid:148) can be obtained as: ๐ฎ (cid:148)โ = argmin ๐ฎ ((cid:156)) โ๐ฎ (cid:157)(cid:158) \๐ฎ (cid:159)(cid:158),(cid:160) (cid:139)๐ฟ๐น & W๐ฎ (X) Y(cid:140) (13) Equation (12) differs from the learning functions in Equation (7) and (8) in a way that the exploration is conducted by maximizing the minimal distance between the new and existing inputs of support points rather than maximizing the prediction variance. Since the learning function in Equation (12) does not rely on the prediction statistics provided by GP, it can be applied to other types of surrogates (SVM, ANN, PCE, etc.). Also, it converts the multi-objective optimization problem into a single objective one by equally weighted sum of objectives. The convex combination of objectives ensures that the optimal solution is a Pareto optimal. Various distance metrics can be used to measure the distance between two points ๐ฎ (X) and ๐ฎ (ยก) in real-valued vector spaces. The Euclidean distance is suggested in this paper, which can be calculated by: (cid:152)๐ฎ (X) , ๐ฎ (ยก) (cid:152) = ยขV (cid:139)u M(X) โ u
M(ยก) (cid:140)
O)M]& (14) Therefore, the minimal distance (cid:152)๐ฎ, ๐ฎ p&,(cid:148) (cid:152) (cid:153)X(cid:132) where ๐ฎ p&,(cid:148) is a set of inputs of support points can be given by: (cid:152)๐ฎ, ๐ฎ p&,(cid:148) (cid:152) (cid:153)X(cid:132) = argmin ๐ฎ ((cid:156)) โ๐ฎ (cid:159)(cid:158),(cid:160) W(cid:152)๐ฎ, ๐ฎ (X) (cid:152)Y (15) The selected input ๐ฎ (cid:148)โ in Equation (13) is evaluated by the implicit performance function to obtain the output ๐ฆ (cid:148)โ = ๐W๐ L& (๐ฎ (cid:148)โ )Y . For the next iteration, the dataset of support points is pdated as ๐ฎ &,(cid:148)N& = ๐ฎ &,(cid:148) โช {(๐ L& (๐ฎ (cid:148)โ ), ๐ฆ (cid:148)โ )} and the corresponding updated surrogated is built as ๐m &,(cid:148)N& accordingly. Stopping criteria
In each iteration ๐ (๐ = 1,2 โฏ ) , the failure probability is estimated by: ๐I -&,(cid:148) = 1๐ (cid:133)& V ๐r - (cid:139)๐ L& W๐ฎ (X) Y(cid:140) ๐ (cid:132) W๐ฎ (X) Y๐ & (๐ฎ) Z (cid:157)(cid:158) X]& , ๐ฎ (X) โ ๐ฎ (cid:133)& (16) where the indicator function based on the surrogate ๐m (cid:148) is written as: ๐r - W๐ L& (๐ฎ)Y = ยฅ1, ๐m &,(cid:148) (๐ L& (๐ฎ)) โค 00, ๐m &,(cid:148) (๐ L& (๐ฎ)) > 0 (17) Note that ๐I -&,(cid:148) may not be a good estimator because ๐ (cid:132) (โ) and ๐ & (โ) can have different supports, but it can be treated as a measure of how well all candidates ๐ฎ (cid:133)& for the first stage have been classified into the safety and failure regions. To decide in which iteration to terminate the Stage 1, early stopping criterion is used in this paper. We stop in iteration ๐ XM& , if ยง๐I -&,(cid:132) (cid:156)ยค(cid:158) โ 1๐ & โ ๐I -&,(cid:148)(cid:132) (cid:156)ยค(cid:158) (cid:148)](cid:132) (cid:156)ยค(cid:158) N&Lโ (cid:158) ยง1๐ & โ ๐I -&,(cid:148)(cid:132) (cid:156)ยค(cid:158) (cid:148)](cid:132) (cid:156)ยค(cid:158) N&Lโ (cid:158) โค ๐ & (18) where ๐ & is a positive integer and ๐ & is the tolerance to the relative error of the estimated failure probability with respect to the mean failure probabilities from previous ๐ & iterations. The setting ๐ & = 5 works well for the proposed framework. And ๐ & controls the accuracy of the built surrogate and is set to 0.01. Another stopping criterion is ๐ XM& exceeding the maximal allowable number of iterations [๐ XM& ] for Stage 1, i.e., ๐ XM& > [๐
XM& ] . At the end of the first stage, we can have (i) a coarse estimator of the failure probability ๐I -&,(cid:132) (cid:156)ยค(cid:158) ; (ii) a coarse surrogate ๐m &,(cid:132) (cid:156)ยค(cid:158) ; (iii) failure samples which are generated by evaluating all candidates ๐ฎ (cid:133)& using ๐m &,(cid:132) (cid:156)ยค(cid:158) . The results can be further utilized in the second stage to improve the estimation of the failure probability as we will discuss in the next section. It is also interesting to find that Stage 1 degrades to AK-MCS [18] if ๐ & (๐ฎ) = ๐ (cid:132) (๐ฎ) in Equation (16) but with a different learning function. .1.4 High dimensionality
Sufficient inputs of candidate support points should be generated to ensure an effective exploration. Considerable attention must be paid to the input size ๐ (cid:133)& when the dimension ๐ is large ( ๐ โฅ 10 ). Specifically, for the space defined by [โ5,5] ) , an input size ) is required to achieve a density of one point per unit volume. Even though these inputs of candidates are evaluated only by the surrogate which requires relatively small computing efforts, the computational burden is likely to be prohibitive when ๐ becomes large for non-high-performance-computer users. An alternative solution for exploration is to combine the search of multiple most probable points (MPP) with FORM as discussed in [1]. Unlike the sampling-based exploration, [1] searches for multiple MPPs via optimization by taking advantage of the first order gradient information of the performance function, making it very efficient even in high dimensional space. As a result, MPPs as well as a FORM estimator of the failure probability can be obtained. The support points generated during the optimization procedure are used to fit the initial surrogate for Stage 2. This exploration solution is suggested when ๐ is large and computers with large memory and computing power are not available. Stage 2 3.2.1
Instrumental probability density function
Similar to [10], in this paper we use Gaussian Mixture (GM) as the instrumental PDF ๐ O (๐ฎ) in Stage 2 in order to take multiple MPPs into account. The GM is made up of ๐พ equally weighted Gaussian distributions centered on the MPPs ๐ฎ ๏ฌ๏ฌ๏ฌ ๐ O (๐ฎ) = 1๐พ V ๐ฉW๐ฎ(cid:151)๐ฎ ๏ฌ๏ฌ๏ฌ,M , ๐Y โ M]& (19) where the mean ๐ฎ ๏ฌ๏ฌ๏ฌ,M is the MPP in the cluster ๐ก and represents the identity matrix. The MPPs ๐ฎ ๏ฌ๏ฌ๏ฌ are determined approximately based on the failure samples from Stage 1 which are denoted as ๐ฎ (cid:133)&,- = ^๐ฎ (X) โ ๐ฎ (cid:133)& ; ๐m &,(cid:132) (cid:156)ยค(cid:158) (๐ L& W๐ฎ (X) Y) โค 0b (or based on [1]). Specifically, first ๐พ -means algorithm [24] is used to partition ๐ฎ (cid:133)&,- into ๐พ clusters. In each cluster ๐ก , the failure sample with the largest ๐ (cid:132) (๐ฎ) is selected as the MPP for this cluster, i.e., ๏ฌ๏ฌ๏ฌ,M = argmax ๐ฎ ((cid:156)) โ๐ฎ (cid:157)(cid:158),ยท(ยค) (cid:139)๐ (cid:132) W๐ฎ (X) Y(cid:140) (20) where ๐ฎ (cid:133)&,-(M) is the set of failure samples in cluster ๐ก . It should be noted that the number of clusters ๐พ may not equal to the number of failure regions. But it is not a problem as shown in the examples if we set a reasonably large ๐พ . The reason is that ๐พ being larger than the number of failure regions simply leads to partitioning the connected failure regions and samples simulated from the GM ๐ O (๐ฎ) will still scatter over all the failure regions. For initialization in Stage 2, the initial surrogate is ๐m O,5 = ๐m &,(cid:132) (cid:156)ยค(cid:158) and the initial dataset of support points ๐ฎ O,5 = ^W๐ p&,(cid:132) (cid:156)ยค(cid:158) , ๐ฆ p&,(cid:132) (cid:156)ยค(cid:158)
Yb = ^W๐ L& W๐ฎ p&,(cid:132) (cid:156)ยค(cid:158) Y, ๐ฆ p&,(cid:132) (cid:156)ยค(cid:158) Yb . To calculate the IS estimator and further refine the surrogate, a large number of samples ๐ฎ (cid:133)O =^๐ฎ (X) ~๐ O (๐ฎ); ๐: 1 โ ๐ (cid:133)O b are generated from ๐ O (๐ฎ) and evaluated by ๐m O,5 . Then, the initial IS estimator in the second stage can be given by: ๐I -O,5 = 1๐ (cid:133)O V ๐r - (cid:139)๐ L& W๐ฎ (X) Y(cid:140) ๐ (cid:132) W๐ฎ (X) Y๐ O (๐ฎ) Z (cid:157)(cid:181) X]& , ๐ฎ (X) โ ๐ฎ (cid:133)O (21) Also, the samples ๐ฎ (cid:133)O serves as the inputs of candidate support points in Stage 2 whose outputs will be evaluated by the original performance function if selected. The sample size of ๐ฎ (cid:133)O is large in general (e.g., ๐ (cid:133)O = 10000 ). Adaptive refinement of the surrogate
Differing from the learning function ๐ฟ๐น & (โ) for the first stage, the learning function in the second stage puts more emphasis on the exploitation task by zooming into the important regions. Observing Equation (21), we can find that in the failure regions the larger the likelihood ratio ๐ (cid:132) / ๐ O is, the more it contributes to the IS estimator. By taking this into account, the learning function in Stage 2 adds another exploitation term to ๐ฟ๐น & (โ) , which can be written as: ๐ฟ๐น O (๐ฎ) = (cid:151)๐m O,(cid:148) W๐ L& (๐ฎ)Y(cid:151) โ (cid:152)๐ฎ, ๐ฎ p&,(cid:148) (cid:152) (cid:153)X(cid:132) โ log โ๐ (cid:132) (๐ฎ) ๐ O (๐ฎ)โ (22) Like ๐ฟ๐น & (โ) , the multi-objective optimization is converted to a single objective optimization problem by weighted convex combination of objectives. In iteration ๐ (๐ = 1,2 โฏ ) , the current surrogate ๐m O,(cid:148) has been fitted using the dataset of support points ๐ฎ O,(cid:148) = ^W๐ pO,(cid:148) , ๐ฆ pO,(cid:148)
Yb = ^W๐ L& W๐ฎ pO,(cid:148) Y, ๐ฆ pO,(cid:148) Yb (for ๐ = 1 , ๐m O,(cid:148) = ๐m
O,5 and
O,(cid:148) = ๐ฎ
O,5 ). The new input of support points can be selected from ๐ฎ (cid:133)O \๐ฎ pO,(cid:148) by minimizing the learning function ๐ฎ (cid:148)โ = argmin ๐ฎ ((cid:156)) โ๐ฎ (cid:157)(cid:181) \๐ฎ (cid:159)(cid:181),(cid:160) (cid:139)๐ฟ๐น O W๐ฎ (X) Y(cid:140) (23) Again, ๐ฎ (cid:148)โ is evaluated by the original performance function. New support point (๐ L& (๐ฎ (cid:148)โ ), ๐ฆ (cid:148)โ ) is appended to ๐ฎ O,(cid:148) to update the surrogate.
Stopping criteria
Similar to the first stage, the iterative procedure terminates when either the IS estimator converges or the number of iterations ๐ XMO exceeds the predefined maximal allowable number of iterations [๐ XMO ] . The convergence criterion is expressed as: ยง๐I -O,(cid:132) (cid:156)ยค(cid:158) โ 1๐ O โ ๐I -O,(cid:148)(cid:132) (cid:156)ยค(cid:181) (cid:148)](cid:132) (cid:156)ยค(cid:181) N&Lโ (cid:181) ยง1๐ O โ ๐I -O,(cid:148)(cid:132) (cid:156)ยค(cid:181) (cid:148)](cid:132) (cid:156)ยค(cid:181) N&Lโ (cid:181) โค ๐ O (24) where ๐ O = 5 and ๐ O = 0.001 . Flowchart
The implementation of S4IS is summarized in Figure 2. To begin the process, a large number of inputs of candidate support points (or candidates) are generated from ๐ & (โ) for Stage 1, which will later be selected adaptively to refine the surrogate in each iteration. To explore the whole parameter space effectively, heavy-tailed distributions ๐ & (โ) such as uniform distributions are suggested. Based on the information about the failure regions obtained in Stage 1, Stage 2 starts with the construction of the instrumental probability density function ๐ O (โ) for IS. Again, after selecting inputs of support points from the candidates generated from ๐ O (โ) , they are evaluated by the performance function to further refine the surrogate. Balancing dynamically between exploration and exploitation is achieved in the two stages by (i) different candidates to be evaluated by the performance function; (ii) different learning functions to select the best inputs of support points. Figure 2: Flowchart of the proposed framework S4IS. Filled rectangles represent where the performance function is evaluated. Illustrative examples
In this section, five examples are discussed to illustrate the effectiveness and efficiency of S4IS for different types of reliability analysis problems, which are featured by the system reliability, the involvement of nonlinear limit-state functions, small failure probability, and moderately high dimensionality. The efficiency is measured in terms of ๐ ยปโฆโโฐ , i.e., the number of performance function evaluations required to achieve a small Coefficient of Variation ( CoV ). For all examples below, we use Gaussian Process (GP) as the surrogate and ten replicates are run to compute the statistics of the final estimator of the failure probability.
Example 1: a series system with four failure regions [18, 25]
This example involves a series system, in which the performance function returns the smallest value of four component performance function. The purpose of this example is to validate the proposed framework for system reliability analysis and multiple failure regions. The performance function is: ๐(๐) = min โฃโขโขโขโขโขโก3 + 0.1(๐ & โ ๐ O ) O โ โ2(๐ & + ๐ O )23 + 0.1(๐ & โ ๐ O ) O + โ2(๐ & + ๐ O )2(๐ & โ ๐ O ) + 3โ2โ(๐ & โ ๐ O ) + 3โ2 โฆโฅโฅโฅโฅโฅโค (25) Initialization
1. Generate a large number of inputs of candidate support points ๐ฎ ๐1 = ๐ฎ (๐) ~๐ ๐ฎ ;๐:1 โ ๐ ๐1 from ๐ ๐ฎ
2. Select the inputs of initial support points ๐ฎ ๐ 1,0 โ ๐ฎ ๐1
3. Evaluate the performance function for ๐ฎ to obtain the initial support points ๐ฎ =๐ โ1 ๐ฎ ๐ 1,0 ,๐ ๐ โ1 ๐ฎ
4. Build the initial surrogate ๐9 based on ๐ฎ
1. Evaluate the current surrogate ๐9 for candidates ๐ฎ ;3
2. Estimate the failure probability ๐= ๐น1,๐
3. Find the input of the next support point by minimizing the learning function ๐ฎ ๐โ = argmin ๐ฎ (๐) โ๐ฎ ๐1 \๐ฎ ๐ 1,๐ ๐ฟ๐น ๐ฎ (๐)
4. Check the stopping criteria. If continue, a) Update the support points by adding the new one ๐ฎ = ๐ฎ โช ๐ โ1 ๐ฎ ๐โ ,๐ ๐ K3 ๐ฎ Lโ b) Update the surrogate ๐9 c) Return to Step 1 Adaptive refinement of the surrogate Initialization
1. Construct the instrumental probability density function ๐ N ๐ฎ
2. Generate a large number of inputs of candidate support points for Stage 2 ๐ฎ ๐2 = P๐ฎ (๐) ~๐ ๐ฎ ;๐:1 โ๐ ๐2 Q
3. The initial surrogate of Stage 2 ๐9 = ๐9 ๐๐ก1 where ๐9 UVW is the final surrogate of Stage 1
Adaptive refinement of the surrogate
1. Evaluate ๐9 N,L for candidates ๐ฎ ;N
2. Estimate the failure probability ๐= XN,L
3. Find the input of the next support point by minimizing the learning function ๐ฎ L โ = argmin ๐ฎ (๐) โ๐ฎ ๐2 \๐ฎ ๐ 2,๐ ๐ฟ๐น ๐ฎ (๐)
4. Check the stopping criteria. If continue, a) Update the support points by adding the new one ๐ฎ N,๐+1 = ๐ฎ
N,๐ โช ๐ โ1 ๐ฎ ๐โ ,๐ ๐ K3 ๐ฎ Lโ b) Update the surrogate ๐9 N,LM3 c) Return to Step 1
Stage 1 Stage 2 here ๐ & and ๐ O are parameter values of two independent standard normal random variables ฮ & ~๐(0,1) and ฮ O ~๐(0,1) , respectively. Figure 3 shows how the support points are adaptively selected to construct the surrogate in S4IS. They are scattering over the all four failure regions especially the four failure boundaries. As a result, as can be seen from Figure 3, the original limit-state function can be well approximated by the surrogate. Table 1 compares the results from S4IS with those from the crude MCS, FORM and AK-IS [18]. The failure probability from MCS ๐I -,๏ฌหR is considered as the ground truth of the problem and the relative error ๐ ห is calculated by: ๐ ห = (cid:151) ๐ m ๐น,๐๐ถ๐ โ ๐ m ๐น (cid:151) ๐ m ๐น,๐๐ถ๐ (26) where ๐I - is the failure probability obtained from the reliability analysis by other methods. Due to the nonlinearity of the limit-state function and the multiple failure regions as show in Figure 3, the failure probability from FORM is inaccurate with ๐ ห = 69.8% . No improvement is observed in AK-IS as the importance sampling in the second stage is based on the FORM results in the first stage. For S4IS, an estimator with a large CoV was obtained in Stage 1 based on the coarse surrogate; the estimator was improved significantly in Stage 2 resulting in ๐ ห = 0.5% and ๐ถ๐๐ = 1.06% . Also, S4IS is more efficient compared with AK-IS, requiring only a total number of 60.6 evaluations of the performance function at the end of Stage 2. It should be pointed out that for system reliability problems, we suggest to use the multi-output surrogate or multiple surrogates to approximate all component performance functions in the system then call the min (for series systems) or max (for parallel systems) function, rather than approximating the aggregated performance function only. The reason is that the min or max function of simple component performance functions may become too complex to be approximated accurately, leading to a large error of the final failure probability estimator.
Figure 3: Adaptively selected support points and the limit-state function using S4IS in Example 1 Table 1: Comparison of the results using different methods for Example 1 Methods MCS FORM AK-IS S4IS
Stage 1 Stage 2 ๐I - L(cid:210)
L(cid:210)
L(cid:210)
L(cid:210)
L(cid:210) ๐ ห -- 69.8% 73.6% 17.1% 0.5% ๐ถ๐๐ ๐ ยปโฆโโฐ (cid:212)
12 71.1 35.7 60.6
Example 2: a nonlinear oscillator [18, 25, 26]
Example 2 investigates the effect of the highly nonlinear limit-state function on the performance of S4IS. The oscillator of interest is shown in Figure 4 and the statistical properties f its six random variables are summarized in Table 2. The performance function can be expressed as: ๐(๐) = 3๐ โ (cid:213) 2๐น & ๐๐ ๐ ๐๐ (cid:217)๐ ๐ก & (27) where ๐ = [๐ & , ๐ O , ๐, ๐, ๐ก & , ๐น & ] + and ๐ = (cid:220) (cid:133) (cid:158) N(cid:133) (cid:181) (cid:153) . Figure 4: Schematic diagram of the oscillator in Example 2 Table 2: Statistical characteristics of random variables in Example 2 Random variable Distribution Mean Standard deviation ๐ & Normal 1 0.1 ๐ O Normal 0.1 0.01 ๐ Normal 1 0.05 ๐ Normal 0.5 0.05 ๐ก & Normal 1 0.2 ๐น & Normal 1 0.2
The results of the reliability analysis using different methods are reported in Table 3. While the relative error of the estimated failure probability from FORM is large, both AK-IS and S4IS results in accurate estimators with
๐ถ๐๐ < 5.0% and ๐ ห < 1.0% . Furthermore, S4IS only requires a total number of 53.3 calls to the performance function, being about half of the number required by AK-IS. Table 3: Comparison of the results using different methods for Example 2 Methods MCS FORM AK-IS S4IS ๐ " ๐ ๐น ๐ก
๐น ๐ก ๐ก๐ก " ๐น " ๐ tage 1 Stage 2 ๐I - ๐ ห -- 9.1% 0.2% 17.2% 0.9% ๐ถ๐๐ ๐ ยปโฆโโฐ (cid:212)
39 91.4 34.2 53.3
Example 3: a multimodal nonlinear function [9]
This example involves a highly nonlinear and multimodal performance function defined by: ๐(๐) = โ (๐ &O + 4)(๐ O โ 1)20 + sin (cid:217)5๐ & (28) Two independent random variables follow the normal distribution ฮ & ~๐(1.5,1) and ฮ O ~๐(2.5,1) . As can been seen from Figure 1, as a whole, the original limit-state function can be approximated by the surrogate accurately with the multiple modals identified. What is interesting to observe in this figure is that the surrogate has a very high accuracy over the important region where the density value of the distribution of random variables is large. It is a natural result of the adaptive infill strategy as in S4IS the second stage zooms into the important region and produces support points therein. Table 3 summarizes the results by using different reliability analysis methods. Both AK-IS and S4IS can achieve the small relative error and variance in the failure probability estimation. In this example, S4IS outperforms greatly AK-IS in terms of efficiency as it takes hundreds of iterations for the optimization procedure in the first stage of AK-IS to find the most probable point (MPP). Figure 5: Adaptively selected support points and the limit-state function using S4IS in Example 3. The dotted lines are the contour lines of the joint probability density function. Table 4: Comparison of the results using different methods for Example 3 Methods MCS FORM AK-IS S4IS
Stage 1 Stage 2 ๐I - ๐ ห -- 277.6% 0.2% 34.5% 1.7% ๐ถ๐๐ ๐ ยปโฆโโฐ (cid:212)
695 985.9 46.2 71.4
Example 4: the influence of reliability levels [11, 27]
The purpose of this example is to investigate the influence of different reliability levels on the performance of S4IS, particularly when the reliability level is very high (small failure probability). The performance function reads: ๐(๐) = min โฉโชโจโชโง๐ โ 1 โ ๐ O + exp โโ ๐ &O
10โ + (cid:217)๐ & (cid:137) ๐ O & ๐ O โญโชโฌโชโซ (29) Two independent random variables in the performance function follow the standard normal distribution. By increasing the constant ๐ , the reliability level can be increased accordingly. o investigate the influence of different reliability levels, three cases ( ๐ = 3 , ๐ = 4 and ๐ = 5 ) are discussed. Comparative results are reported from Table 5 to Table 7. S4IS performs consistently well for different reliability levels, in which the coarse estimator ๐I - in Stage 1 can be improved significantly in Stage 2. In this example, increasing the reliability level from ๐I - = 3.470 ร 10 L(cid:210) to ๐I - = 9.485 ร 10 Lล only add a little computational effort (from ๐ ยปโฆโโฐ = 72.8 to ๐ ยปโฆโโฐ = 118.6 ) to achieve the high accuracy. Again, FORM and AK-IS tend to be mistaken by the only MPP they can find, resulting in the inaccurate estimators. Table 5: Comparison of the results using different methods for Example 4 ( ๐ = ๐ ) Methods MCS FORM AK-IS S4IS
Stage 1 Stage 2 ๐I - L(cid:210)
L(cid:210)
L(cid:210)
L(cid:210)
L(cid:210) ๐ ห -- 61.1% 57.9% 13.6% 1.8% ๐ถ๐๐ ๐ ยปโฆโโฐ (cid:212)
7 97.6 44.9 72.8
Table 6: Comparison of the results using different methods for Example 4 ( ๐ = ๐ ) Methods MCS FORM AK-IS S4IS
Stage 1 Stage 2 ๐I - L(cid:236)
L(cid:236)
L(cid:236)
L(cid:236)
L(cid:236) ๐ ห -- 65.5% 50.8% 1.2% 0.6% ๐ถ๐๐ ๐ ยปโฆโโฐ (cid:212)
7 110.3 54.4 83.2
Table 7: Comparison of the results using different methods for Example 4 ( ๐ = ๐ ) Methods MCS FORM AK-IS S4IS
Stage 1 Stage 2 ๐I - Lล Lล Lล Lล Lล ๐ ห -- 69.8% 76.0% 36.1% 4.7% ๐ถ๐๐ ยปโฆโโฐ (cid:238)
7 92.4 63.4 118.6
Example 5: the influence of dimensionality [11, 28]
The last example is characterized by an analytical performance function where the number of random variables ๐ can be changed without altering the reliability level a lot. The performance function reads as follows: ๐(๐) = W๐ + 6โ๐Y โ V ๐ X)X]& (30) This example involves ๐ independent lognormal random variables with the mean values 1 and standard deviations 2. Three cases ( ๐ = 2 , ๐ = 10 and ๐ = 50 ) are investigated by using different methods for this example. The results are summarized from Table 8 to Table 10. In S4IS, for moderately high dimensional cases ( ๐ = 10 and ๐ = 50 )), the combination of the search of multiple MPPs and FORM [1] is adopted in Stage 1 as discussed in Section 3.1.4. It is shown that on this example, the total required number of evaluations ๐ ยปโฆโโฐ in S4RS increases slowly with the number of random variables. For ๐ = 50 , less than two hundred evaluations are required, which is generally affordable in the setting of the simulation of the real-world systems. The results using other reliability analysis methods are compared to S4IS. For FORM, though the performance space in the original space is linear, it becomes nonlinear when the original random variables are transformed into the standard normal ones. As it can be seen, FORM fails when the dimension is high since the nonlinearity of the limit-state function in the U-space. AK-IS can produce the accurate estimators of the failure probability for this example even when the dimension is high, but AK-IS requires much more calls of the performance function than S4IS (for ๐ = 50 , ๐ ยปโฆโโฐ = 1845.2 in AK-IS). It is not surprising as the stopping criterion in AK-IS rely on the convergence of the minimum U-function value over the all IS samples while S4IS checks the convergence of the estimator directly. Table 8: Comparison of the results using different methods for Example 5 ( ๐ = ๐ ) Methods MCS FORM AK-IS S4IS
Stage 1 Stage 2 I - L(cid:210)
L(cid:210)
L(cid:210)
L(cid:210)
L(cid:210) ๐ ห -- 22.0% 0.04% 0.2% 0.1% ๐ถ๐๐ ๐ ยปโฆโโฐ (cid:212)
20 59.0 16.4 23.9
Table 9: Comparison of the results using different methods for Example 5 ( ๐ = 10 ) Methods MCS FORM AK-IS S4IS
Stage 1 Stage 2 ๐I - L(cid:210)
L(cid:210)
L(cid:210)
L(cid:210)
L(cid:210) ๐ ห -- 63.4% 1.2% 44.5% 0.2% ๐ถ๐๐ ๐ ยปโฆโโฐ (cid:212)
35 678.2 38 48.6
Table 10: Comparison of the results using different methods for Example 5 ( ๐ = 50 ) Methods MCS FORM AK-IS S4IS
Stage 1 Stage 2 ๐I - L(cid:210)
L(cid:137)
L(cid:210)
L(cid:137)
L(cid:210) ๐ ห -- 92.0% 1.6% 88.0% 1.0% ๐ถ๐๐ ๐ ยปโฆโโฐ (cid:212)