A quantile-based method to efficiently determine model response at boundary probabilities

Quantiles are values of a function associated with specific probabilities. Two very specific quantiles, that is, ±3, are usually used for estimating probabilistic distribution of the function. Aiming at the engineering complex implicit function, the article provides an algorithm for calculating these values efficiently and rapidly. First, build kriging model with initial high-efficiency samples and the model is taken as the approximate initial analytical expression of the original function. Then, seek for the approximate most probable point with kriging model and performance measure approach. The most probable points found every time and their corresponding function responses are taken as new samples to reconstruct the built kriging model. After several times of reconstructing kriging model, the convergence value for the corresponding quantile can be obtained. Finally, the illustrative cases in this article demonstrate that the quantile-based method can provide accurate results with high efficiency and in this way, the quantile convergence value can be obtained by a few samples.


Introduction
In engineering applications, uncertainties exist in many problems. For example, machining errors in mechanism links lead to uncertainties; internal defects of material lead to uncertainties of physical parameters; and external environment leads to load uncertainties. Under the influence of these uncertain factors, many problems' output also becomes uncertain. Considering its significance to engineering safety design, the uncertainty research has drawn wider attention and given rise to branches such as reliability research and robustness research.
The methods of uncertainty research include probabilistic approaches and non-probabilistic approaches. Non-probabilistic approaches include interval analysis, 1,2 fuzzy theory, 3 possibility theory, 4-6 convex theory, 7 and epistemic parameter uncertainty. 8 Relatively speaking, probabilistic approaches are more widely used. For example, Hamon et al. 9 conducted a sampling-based probabilistic security evaluation (PSE) of a high-power wind system. Sarsria and Azrar 10 studied the dynamic response of large structures under the influence of uncertainty factors with the perturbation 1 method. Sepahvand et al. 11 proposed to use polynomial chaos expansion (PCE) to quantitatively analyze the uncertainties of stochastic systems. Zhao and Ono 12 put forward a universal program to analyze the reliability by the first-order, first-moment method and secondorder, second-moment method. When Lee et al. 13 studied the reliability optimal design, they proposed an uncertainty calculation method of nonlinear and multidimensional systems, which is based on the most probable points (MPPs). Based on the dimension reduction method, Chowdhury and Rao 14 proposed a reliability calculation method of hybrid high-dimensional model under the influence of random factors. Wei et al. 15 proposed an active learning procedure combined with the multiple-response Gaussian process model and Monte Carlo simulation to efficiently and adaptively produce surrogate models for the failure surfaces of systems. Xiao et al. 16 proposed a new adaptive surrogate model-based efficient reliability method to compute structural reliability.
Generally speaking, probabilistic approaches are mainly composed of the following five categories: (1) simulation-based methods: Monte Carlo simulations and improved variance deduction method; [17][18][19] (2) local expansion-based methods: Taylor series, 20 perturbation method, 21 and so on; (3) functional expansionbased methods: Neumann expansion, 22 PCE, 23 and so on; (4) MPP-based methods: first-order reliability method (FORM) 12 and second-order reliability method (SORM); 24 and (5) numerical integration-based methods: dimension reduction. 25 However, despite the big development of uncertainty researches, Methods (2)-(5) are only valid for certain problems under certain conditions. The above four methods are therefore not universal due to their limitations. Fortunately, Method (1) is universal and easy to use. Presently, there are two difficulties in calculating uncertainty for complex largescale engineering problems: (1) the function is a black box and implicit; (2) the calculation requires business software such as finite element, resulting in a large amount of calculation. Under these two circumstances, only Method (1) is suitable. But the sample selection of Method (1) is random, resulting in slow convergence. In consequence, a considerable number of samples are needed to ensure that the results are accurate. Besides, it takes a long time to solve the large-scale and complicated problem every time. In case of conducting response calculation of a large number of samples, the complicated engineering problems will need to be solved repeatedly. Usually, the calculation cost is too high to carry out.
In view of the above, this article presents an efficient analytical method based on quantiles to determine the solutions at critical boundary probabilities. The method contains three ideas: first, evaluating the probability response range of the problem by means of the quantile; second, new samples are chosen close to the boundary probability; and third, reconstructing a surrogate model with new samples in order to make full use of the information provided by new samples and accelerate the surrogate model to approximate the real function and therefore make sure that the convergence of the response is accelerated. The article is organized as follows. Section ''Concept of determining the uncertainty of function probability with quantile'' introduces the concept of determining the probabilistic uncertainty of the function with quantile. Section ''Fast method of calculating function quantile response value'' presents a fast method of calculating the quantile response value of the function. Section ''Case analysis'' verifies the method proposed in this article through case analysis. Finally, section ''Conclusion'' summarizes the method proposed in this article.

Concept of determining the uncertainty of function probability with quantile
For any complex model, suppose its function to be Y = f (X), and X = ½X 1 , X 2 , . . . , X n is a random vector. Y is the output of this complex model. As Figure 1 shows, the input parameter X is uncertain, resulting in the fact that the output Y is also uncertain. Figure 2 presents the probability distribution of function f . p 1 is the quantile of the left tail, and p 2 is the quantile of the right tail. In this article, assume that reliability indexes are equal to 6 3. In statistics, these are z-values for a standard normal distribution. Then,   Figure 2, the probability of the shaded part is 0.9973. Therefore, as long as the function response values f p1 and f p2 corresponding, respectively, to p 1 and p 2 of function f are determined, the probability that function f falls in the range ½f p1 , f p2 is 99.73%.
The quantile response value of the function can be calculated by the performance measure approach (PMA method) as shown in equation (1) min f (u) where u is the mapping value of the random vector X in the standard normal space. When the sub-components of the random vector X are correlated, they can be transformed into independent normal random variables by Rosenblatt 26 transform. In order to better focus on the method proposed in this article, it is assumed that the random vector X is an independent normal random vector. Through the standard normalization of equation (2), the random vector X will be transformed to the standard normal vector u where m X i and s X i are, respectively, the mean and standard deviation of the ith component X i of the random vector X.
which is the quantile f p1 corresponding to p 1 . The quantile f p2 corresponding to p 2 can be obtained by solving equation (3) max f (u) s:t: Obviously, equation (3) There are several advantages of using quantile to determine the function uncertainty. (1) Quantiles describe values of the function associated with defined probabilities. (2) Using the PMA method, the function response value at the quantile can be obtained with a small number of samples. Therefore, the quantile-based method to evaluate function probability distribution needs less calculation and achieves higher efficiency.

Construction of high-efficient initial sample set
The high-dimensional model representation is a generic and quantitative tool for model assessment and analysis. Using this method, the high-dimensional relationship between the model input variable and the output variable can be obtained. 14 And the initial sample used in this method greatly improves its modeling accuracy. Therefore, by referring to the way of selecting initial samples in the above-mentioned high-dimensional model representation method, the initial samples selected of the article in the standard normal space are as follows where m i and s i are, respectively, the mean value and standard deviation of the ith sub-component u i of the random vector U; j ranges from 1 to M, where M is generally set to be 3, 5, 7, or 9; and floor() indicates that the values in parentheses are rounded to the nearest positive integer. For the following cases shown in the article, M = 3. In the standard space, m i = 0 and s i = 1.
Samples of the random vector X in its original space can be transformed by equation (2). Finally, the function response value of the above sample is calculated by performing the actual model simulation and then the samples and their function response values will be combined into initial sample set B 0 = fX, Yg.

Kriging model
After the initial sample set B 0 is established, the author uses kriging method to construct the kriging model for the function. In general, the kriging model consists of two parts: the parameter model and the nonparametric random model, as shown in equation (5). Assume that the input of the model is x = ½x 1 , x 2 , . . . , x n T and the corresponding output is y(x), then the corresponding kriging model is expressed as where F(x) adopts the first-order regression model and is expressed as follows where r 0 , r 1 , . . . , r n are regression parameters.
In equation (5), e(x) is usually a polynomial expression, with the average 0, variance s 2 , and its function structure as follows In this article, y j (u, x iÃ h À x h ) adopts the Gaussian form, and then R i (u, x iÃ , x) can be expressed as In equation (7), g Ã = R À1 (Y À Kr) is the matrix with dimensions being m 3 p, and its child element is K ij = T j (s i ). Y n 3 1 is the output response corresponding to the input sample point x Ã . R is the matrix with dimension being m 3 m, and the child element R ij takes the Gaussian form as shown in equation (8). For any two sample points x iÃ and x jÃ , their R ij form is Fast calculation of function output for the quantile based on kriging modeling method Build the initial kriging model with the initial set of samples B 0 based on the kriging method and assume that the initial model is f kri 0 . The subscript ''0'' indicates the initial modeling model. Equation (1) can be rewritten as follows where k = 0, 1, 2, . indicates the times of model construction. By solving equation (11), the corresponding optimal solution u Ã k and its corresponding function response value f kri k (u Ã k ) (hereinafter abbreviated as f Ãkri k ) can be obtained.
Add the optimal solution u Ã k and its corresponding function response value f Ãkri k to sample set B kÀ1 and mark the new sample set as B k . Then, use the kriging method to build the new kriging model f kri k with the new sample set.
Then compare, respectively, the convergence of the optimal solution u Ã k and u Ã kÀ1 , as well as f Ãkri k and f Ãkri kÀ1 , according to equations (12) and (13) to judge whether the calculation is convergent. e is a relatively small value. In the later instance of this article, e is taken as 0.05 f Ãkri When equation (12) or (13) is satisfied, the calculation is stopped; otherwise, the above calculation is repeated. Finally, the resulting convergence value is the function response value of the corresponding quantile. Similarly, the convergence value of f p2 in equation (14) can be calculated according to the calculated method of f p1 Case analysis Case 1 Example 1. Figure 3 shows a roof truss structure with a uniform load q and a node load P. 22,27 The limit state function of the truss structure is g(x) = 0:2 À D C (Ac, Ec, As, Es, l) The above formula is called the limit state function which indicates that the downward deflection D C of the calculated point C does not exceed 0.2. If the downward deflection D C is bigger than 0.2, the structure is considered to be unsafe. In D C , the Ac, Ec, As, Es, and l, respectively, represent the cross-sectional area, elastic modulus, and length of the concrete and steel pole. The above random variables are normally distributed, and the specific distribution parameters are shown in Table 1. The above structure is relatively complex, and its limit state function is not analytic. In the calculation of the structure, the finite element method is used to calculate its limit state function.
By substituting the average value in equation (15), the approximate mean value of the function response is 0.1766. Figures 4 and 5, respectively, represent the iteration calculation of the function response values of the left and right quantiles in Case 1. As for the left quantile g 1 , the iteration number is 3, in which one time is for initial construction and two times for reconstructions. For initial construction, the number of initial samples is 13 (= 2 3 6 + 1). While for reconstructions, the two extra new samples are derived. Hence, the number of calculated samples is 15 (= 13 + 2). The same is true for the right quantile g 2 . Since the initial samples are also used for constructing the right quantile g 2 , the repetition for computing their response values can be saved. Hence, the total number of calculated samples is 17 (= 15 + 2). Hence, the proposed method is quite efficient for obtaining accurately converged results.
In order to verify the correctness of this method, the Monte Carlo method (hereinafter referred to as the MC method) is used to compute the structure, and then its mean value and standard deviation are calculated. Finally, the left and right quantiles are obtained according to equations (16) and (17), with the details shown in Tables 2 and 3. The MC calculation of the left and right quantiles in Cases 2 and 3 is the same as Case 1 In Table 3, the calculated results of the MC method are used as the standard results and the relative errors are obtained by calculating the relative errors between the calculated results and the standard results. It can be found that the calculation results between the MC method and the presented method in this article are close to each other. In addition, the number of samples required for this method is much less than those of the MC method. Therefore, the calculation method used in this article is very efficient.     Figure 6 shows a double-lap joint design of rubbermodified epoxy-based adhesive. 28 The design consists of aluminum outer adherent and an inner steel adherent. The completed bond is subjected to an axial load p. S a is the lap-shear strength of adhesive. The parameters of probability distribution of the independent random variables are shown in Table 4. The limit state function of the structure is the strength of the joint, that is, the maximum shear stress is less than the shear strength of the binder. When x = 0:5, the shear stress of the structure is the largest, and the analytic expression of its limit state function is as follows.
The limit state function is the safety margin for strength requirement of the joint, which is defined by the difference between the lap-shear strength of adhesive and the maximum shear stress. The equation is obtained at x = 0:5 where the maximum shear stress occurs. The limit state function is given by where t max = t(x) By taking the mean value of equation (18), the approximate mean value of the function response is 524.24. Figures 7 and 8, respectively, represent the iterative computation of the left and right quantile response of the function in Case 2. On the basis of 19 (2 3 9 + 1) initial samples, the left quantile g 1 and the right quantile g 2 converge to 253.6 and 782.64, respectively, after three iterations. As for the left quantile g 1 , one time is for initial construction and two times for reconstructions. For initial construction, the number of initial samples is 19 (= 2 3 9 + 1). While for reconstructions, the two extra new samples are derived. Hence, the number of calculated samples is 21 (= 19 + 2). The same is true for the right quantile g 2 . Since the initial samples are also used for constructing the right quantile g 2 , the repetition for computing their  Table 4. Distribution of the independent random variable parameters in Case 2.

Random variables Probability distribution type Average
Standard deviation Normal distribution 2e5 lb/in 2 2000 lb/in 2 x 6 (b) Normal distribution 1 in. 0.01 in. response values can be saved. Hence, the total number of calculated samples is 23 (= 21 + 2). Hence, the proposed method is quite efficient for obtaining accurately converged results ( Table 5). The MC method is adopted for computing the structure and its computed results are used as the standard for comparing the results of the method used in this article. More detailed results are shown in Table 6. It can be found that the calculated results of this method are similar to those of the MC method, but the computational efficiency is much higher than that of the MC method. Figure 9 shows a structure consisting of three beams with two degrees of freedom. The details of the structural parameters are listed in Table 7.

Case 3
The structure requires that the maximum deformation under loads cannot exceed 0.2 m, as follows By substituting the average value in equation (19), the approximate mean value of the function response is 0.1548. Figures 10 and 11, respectively, represent the iterative computation of the left and right percentile response of the function in Case 3. On the basis of 15 (= 2 3 7 + 1) initial samples, they respectively converge to 0.1188 and 0.1718 after three and four times of iterations. The details of the results are shown in Table  8. A total number of 20 (= 15 + 3 + 2) calculated samples are needed to make sure the calculated results accurately converge to the real function response value for the corresponding quantile. Here, 15 represents the calculated number for initial samples; 3 and 2 represent the new generated samples when reconstructing the left quantile g 1 and right quantile g 2 . All in all, the case shows that the proposed method is quite efficient.
The structure is also computed by the MC method, and the corresponding results are used as the standard to compare the calculated results of the presented method. More detailed results are shown in Table 9. It can be found that the calculated results of this method are similar to those of the MC method, but the computational efficiency is much higher than those of the MC method.

Conclusion
In this article, a framework of calculating the function based on the quantile method is used to describe the stochastic uncertainty of the function. The main features of this method are as follows: (1) use the quantile conception to estimate the probability distribution of the function, namely, to determine the probability distribution range of the function based on the idea of the quantile. (2) Purposely select new samples, namely, according to certain strategy, that is, the PMA method, the new samples are chosen, thus overcoming the randomness in sampling and conducive to accelerate the calculation convergence. (3) Reconstruct the function, namely, by reconstructing the kriging model; it can fully absorb the information provided by the new sample and accelerate to approximate the real response value of the function at the quantile. The results of this article show that the method is accurate and reliable, the actual calculation time of the function is much lower than the MC method, and the calculation efficiency is extremely high.   The above-mentioned three features of this method have a strong reference. If this method is used, the computational efficiency of reliability and robustness optimization design can be greatly improved because these designs are involved with the uncertainty assessment of functions. Currently, the reliability and robustness optimization design of complex engineering still has     problems, such as black boxes and excessive computation. In the future, related researches will be carried out to further expand the application of this method.