Quantile-based versus Sobol sensitivity analysis in limit state design

Quantile-based versus Sobol sensitivity analysis in limit state design KALA, Z. Structures 202O, vol. 28, December 2020, pp. 2424-2430 ISSN: 2352-0124 DOI: https://doi.org/10.1016/j.istruc.2020.10.037 Accepted manuscript © 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license http://creativecommons.org/licenses/by-nc-nd/4.0/ dspace.vutbr.cz Quantile-based versus Sobol sensitivity analysis in limit state design Highlights - Rd as the design quantile of resistance (R) of a steel beam under lateral-torsional buckling. - Quantile-oriented sensitivity analysis (QSA) subordinated to contrasts between R and Rd. - Comparison of QSA with classic Sobol sensitivity analysis (SSA). - Identification of similar sensitivity ranking of input random variables by QSA and SSA. - Different interaction effects are obtained from QSA (strong) and SSA (weak). Abstract In limit state design, the reliability of building constructions is generally verified using design quantiles. The design resistance of a structure is explicitly expressed as a low quantile of the cumulative distribution function of resistance. The aim of this article is to show the connections and differences between quantile-oriented sensitivity analysis subordinated to a contrast and classic Sobol sensitivity analysis. Changing the fixed input variable causes synchronous change in the quantile and mean value, but how do the results of these two sensitivity analyses differ? The question is whether or not the changes around the design quantile (measured by contrast indices) are similar to the changes around the mean value, which are measured using Sobol’s indices. Comparison is performed on a case study, where the resistance of the structure is expressed by a non-linear function, the inputs of which are random material and geometric characteristics of the structure. The nondimensional slenderness is a deterministic parameter, which changes the influence of input variables on the resistance as the model output. It was concluded upon comparing the results of both sensitivity analyses that the rank of the most important variables is the same for both low and high slenderness and is similar for intermediate slenderness. However, the interaction effects are very different. The identification of insignificant variables is the same. Other significant similarities and differences between both types of sensitivity analyses are presented in the article. Keywords: Sensitivity analysis; structural reliability; Eurocodes; design quantiles; stability; buckling; steel 1. Introduction Eurocodes [1] assess the reliability of building structures using design quantiles of load and resistance, which are more practical for structural design than failure probability [2][3]. The basic reliability targets for design values in ultimate limit states recommended in [1] are based on a semi-probabilistic approach [4][5]. Structural design is reliable if the design resistance (lower quantile) is higher than the design load (upper quantile). This concept evaluates reliability using two states, i.e., satisfactory (design resistance > design load) or unsatisfactory (design resistance < design load). From a design point of view, this reliability assessment is sufficient even without the analysis of the failure probability Pf or the reliability index  . A structure is either designed according to standards and is therefore sufficiently reliable, or such a design must be ensured. Economically and reliably optimal design means that the design resistance is equal to the design load. Higher reliability than is prescribed by standards is possible, however, higher economic costs of construction would be required (material, etc.). Reviews of sensitivity methods in interdisciplinary contexts are offered in [6], [7], [8] and [9]. In terms of utilization in the construction industry, the methods of sensitivity analysis (SA) can be divided into global sensitivity analysis [10], [11], [12], [13], etc., reliability-oriented sensitivity analysis ROSA [14], [15], [16], [17], etc., and other specific types of SA, see, for e.g., [18], [19], [20] and [21]. The first types of analyses focus on the distribution or certain moments of the output function while ROSA considers Pf,  or design quantiles as the quantity of interest. Comparative case studies [22] have shown that the differences between indices in the first and second group as well as between the groups are large and the information value of ROSA indices is unclear. In civil engineering, SA methods are often applied in a conventional manner, which takes advantage of the knowledge that the methods of the first group are partially empathic to reliability, see, for e.g., [23] and [24]. It can be argued that the input variable (imperfection) could have an influence on the distribution of the limit state function but not on Pf and conversely. However, such cases are not common in engineering practice. The results of reliability analysis tend to be more or less consistent with the results of Sobol sensitivity analysis (SSA), see, for e.g., analyses of the failure mode and load of steel frames [25]. The question is, what is the information value of the indices of the first group if the reliability analysis is the main subject of interest? Intuitively, it can be surmised that when the input variable influences the variance, it can also have some influence on structural reliability. A basic comparison of SA methods in terms of Pf has been performed in [22]. In this article, we will focus on the random resistance R, for which we will evaluate the SA of the design quantile of the resistance Rd together with SSA. The structural resistance R represents the model output, which is generally non-linearly dependent on the inputs. The aim of this article is to show the connections between quantile-oriented sensitivity analysis subordinated to a contrast [14] (referred to as QSA in this article) and classic SSA. The results of QSA are compared with the results of SSA using a simple case study, where the quantity is used in limit state design (LSD). This article examines whether or not the change in R that occurs around the design quantiles is similar to the change around the mean value, which is measured by Sobol indices. The design resistance is explicitly expressed as the  th quantity from the cumulative distribution function of R. The subject of interest is the evaluation of inputs contributing to the extreme output values. It can be noted that the link between quantile-based SA and Sobol main effect sensitivity indices has already been studied in [26], however, not with the use of GSA [14], which is applied here. The GSA applied here is based on contrasts that measure the average distance between the quantile and the output in comparison to the measures proposed in [26], which use the average distance between the quantities. Additional studies of quantile-oriented SA subordinated to a contrast can be found, for e.g., in [27], [28] and [29]. 2. Sensitivity measures based on quantiles of output The model is represented by a mapping f which relates the domain of inputs to the output space: R = f(X1, X2, ..., XM), (1) where R is the model output (static resistance) and input variables (X1, X2, ..., XM) are described using probability density functions (pdfs), which reflect the uncertain knowledge of the system under analysis. The Goal Oriented Sensitivity Analysis is based on contrast functions [14]. The contrast function  associated with  -quantile can be written using parameter  as: (  ) = E(  (R,  )) = E((R -  )(  R<  )). (2) The estimator of  -quantile  * is given as  * = Argmin (  ). Based on [14] the first-order quantile contrast index Qi (first-order or main sensitivity index subordinated to the contrast) can be defined as:            ψmin ,minψmin i i XRψEE Q   (3) where the numerator is the contrast variation due to Xi. Qi is the sensitivity index of the estimator of  *. The minimum value of contrast (  ) can be calculated, for e.g., using N runs of the Latin Hypercube Sampling method (LHS) [30] and [31] as:             N nnXnXnXf MM nXnXnXf N1...,,, * 21 * 21 1...,,, 1 ψmin    (4) where  * is the estimator of  -quantile. The  -quantile  * can be estimated from LHS runs so that  runs of R are smaller than  * and (1-    runs of R are greater than  *, see the example shown in Fig. 1. For example, the estimate of 0.1-percentile (0.001-quantile) is practically calculated as the 400th smallest value in the set arranged in an ascending order with N = 400000 runs of LHS, see Fig. 1. Fig. 1. Example of the estimation of the -percentile from 400000 runs of LHS. The second member in the numerator (3) can be determined using two sets of the LHS method. K random realizations of Xi, i.e., Xi(1), ..., Xi(k),…, Xi(K) are generated in the first set. Then N random realizations of the vector X~i are generated for each realization Xi(k), k = 1,…, K (all variables but Xi). For fixed Xi we can calculate               N nknkXkXf iii ii knkXkXf N kmXRψE 1,, * ~* ~ 1,, 1 ,min    (5) where conditional  -quantile  *(n) is evaluated analogously as in (4) with the difference that the runs of R are obtained for N random realizations of variables X  i and fixed (non-random) Xi. For K runs of Xi we then obtain       K k ikm K XRψEE 1 1 ,min   (6) The second-order quantile contrast index Qij is defined as:        ji ji ij QQ XXRψEE Q       ψmin ,,minψmin (7) The higher-order quantile contrast indices can be expressed analogously. The sum of all sensitivity indices must be equal to 1. 1... ...123  M iijjk ijk iij ij i iQQQQ (8) The sum of all indices is 2M-1. If we substitute (2) with the contrast function (  ) = E(R -  )2, then the equations presented above lead to the classic Sobol decomposition [10] and [11], where the first-order sensitivity index Si is defined as:             i ii iXRERcorr RV XREV RV XRVERV S, 2    (9) where corr is Pearson correlation coefficient. Sobol higher-order sensitivity indices Sij, Sijk, etc. are determined analogously. It can be noted that (2) is not of quadratic type as the contrast associated with Sobol. The substitution of corr with Spearman's rank correlation or Kendall's  also leads to the decomposition, where the sum of all indices is 1. However, the size of indices is generally a little different and their significance may be discussed. The paper [14] is the first step towards a generalized theory of sensitivity analysis, which views Sobol indices as a special case that is based on quadratic contrast. Analogously as in SSA [32] we can introduce the socalled total index in QSA as:            ψmin ,minψmin 1~i Ti XRψEE Q   (10) The total effect index measures the total contribution to the output contrast due to variable Xi, i.e., its firstorder effect including all higher-order effects due to interactions. The total effect is given as the sum of all terms in (8) where variable Xi is considered. As will be shown below, the introduction of (10) is very useful if the number of input variables is high and the higher-order effects due to interactions are noticeable. 3. The target reliability and design quantiles The aim of this article is QSA of design resistance Rd, which is calculated as 0.1 percentile of the static resistance R. The calculation of 0.1 percentile is based on the semi-probabilistic approach of standard [1]. For a simple case of Gaussian distributed variables of F on the load action side and R on the resistance side [33], the limit state function can be written as: 0 FRG (11) where F and R are statistically independent variables for which we can write: 2 F 2 RGFRG ,   (12) where  is the mean value and  is the standard deviation. The reliability index  can be evaluated in relation to the probability of failure Pf as: G    G  with      1 f P (13) where (•) is the cumulative normalized Gauss distribution. For the reference period of 50 years, the target reliability recommended in [1] for ultimate limit state is  d = 3.8, Pf = 7.2E-5, see Table C2 in [1] and/or [34]. In order to answer the question of which value on the load action side and which value on the resistance side result in the desired level of reliability, we have to introduce weighting factors  F (for the load action effect) and  R (for the resistance): G R R G F F         , (14) The so-called design values Fd and Rd can then be calculated as: RdRR FdFF        d d R F (15) If the design is reliable, it must hold that dd RF . The probability that resistance RRd can be written as:   dR R RRdRR              d RRP (16) Upon substitution of βd = 3.8 into (16) we obtain a probability of 0.001183, which is approximately 0.1 percentile of R, see Fig. 2. An alternative equation to (15) is listed in [1] for variables, which do not have a Gaussian distribution, but have a different statistical distribution, such as Gumbel or Log-Normal pdfs. In this article Rd is calculated in all cases as 0.1 percentile non-parametrically based on LHS simulations, see Fig. 1. Fig. 2. The random and design resistance. 4. The computation of static resistance The elastic load carrying capacity MR was derived in [35], see also [36]. In this article the computational model of elastic load carrying capacity is taken from [35] with the difference that instead of MR we introduce the symbol R, so that we are consistent with the model output (resistance R) in the preceding paragraphs. The relatively fast response of the non-linear function (17) in the closed-form permits the use of high numbers of runs of the LHS method and gives accurate and highly transparent results even without using metamodels and supercomputers. The resistance R (17) of an I-beam decreases with its increasing length L due to the LTB phenomena [35]. The initial axial deformation increases when the loading bending moment M increases, see Fig. 3. Failure occurs when the longitudinal stress reaches the yield strength fy, which occurs when M approaches R at x = L/2.   zcr cr zcr crcrcr WM QMQQ WM QMQQMQQMQQQ R 4 22 4 44244 241 2 2 2 24 2 4341 2 1      where  bIW hIW LEIP hPMMea WahPQ WaPWMQ WaPWMQ WWMfQ zz yy zz zcrcrv yvz yvzzcr yvzzcr zycry /2 /2 / 2/2 22 00 0 2 4 03 02 1          (18) The significance of all variables in equations (17) and (18) are explained in detail in [35]. The input variables are referred to as imperfections in engineering terminology [37]. Generally, these deficiencies arise during the production of the structure, its assembly and subsequent use and have an influence on the resistance, especially in cases where the structure is exposed to loss of stability, such as flexural or lateral-torsional buckling [38]. The resistance R (17) is in complete agreement with the elastic load carrying capacity obtained from the computational model created using the software Ansys, which was verified for 1.20  LT  [39] and [40], where LT  is the non-dimensional slenderness [41]. The calculation of the second moment of area is based on an idealized geometry of a doubly-symmetric IPN200 section, see Fig. 3. The random variabilities of geometric characteristics h, b, t1 and t2 (geometric imperfections) of the IPN200 section are introduced according to experimental researches [42] and [43]. The equation for the elastic critical moment Mcr is introduced according to [35] and [34], where shear modulus is G = 0.5·E/(1+  ). Fig. 3. Lateral-torsional buckling and geometric imperfections. The shape of initial geometric imperfections is introduced according to the first eigenmode (bow imperfection) of LTB, see Fig. 3. The random variability of amplitude e0 is based on the assumption that 5% of random samples of e0 lie outside the tolerance limits ±L/1000, see, for e.g., [35] and [34]. Statistical characteristics of yield strength are considered according to [42] for steel S235 and [44] for steel S355, see also the discussion in [45]. The statistical characteristics of the other random imperfections in Table 1 are justified in [46]. The influence of residual stress was neglected. All random variables listed in Table 1 are statistically independent. Table 1. Statistical characteristics of input imperfections [42] [44] No. Symbol Characteristic Density Mean  Standard deviation  1. t2 Flange thickness Gauss 11.3 mm 0.518 mm 2. fy Yield strength, S235 Gauss 297.3 MPa 16.8 MPa Yield strength, S355 Gauss 393.8 MPa 22 MPa 3. E Modulus of elasticity Gauss 210 GPa 10 GPa 4. e  Initial curvature, S235 Gauss 0 L1/1960 *  Initial curvature, S355 Gauss 0 L2/1960 * 5. h Cross-section height Gauss 200 mm 0.885 mm 6. b Flange width Gauss 90 mm 0.888 mm 7. t1 Web thickness Gauss 7.5 mm 0.293 mm 8.  Poisson's ratio Gauss 0.3 0.009 where 432 139.095.175.015.2 LTLTLTLT L   for S235 [46], 432 215.095.036.07.1 LTLTLTLT L   for S355. Dependencies between L1 and LT  , L2 and LT  listed in Table 1 approximate the procedure in [41]. The analytical function between L and LT  is derived in [47]. Slenderness is one of the basic parameters in the design of steel structures. The parametric study uses LT  to evaluate the effect of parametric change of LT  on sensitivity indices, where LT  changes with the step of 0.01. The length of the member L is assigned according to steel grade as L1 or L2. 5. Sensitivity analysis results The LHS method is used to evaluate all sensitivity indices. Each index Qi (3) is evaluated using K = 4000 runs of m(j) (5). Each realization of m(j) is evaluated using a set with N = 400000 runs of R, in which  *(j) represents the conditional 0.1 percentile, which is the 400th smallest value in the set arranged in an ascending order, see, for e.g., Fig. 1. The denominator in (3) is calculated according to (4) for N = 400000 runs. Indices Si, Qij, Qijk, Qijkl and QiT are calculated analogously using the same values K and N. 5.1 QSA for 8 input imperfections from Table 1 QSA with M = 8 imperfections (Table 1) leads to 28-1 = 255 quantile contrast indices, which is too much for practical calculation. Therefore, all global effects are described using indices Qi and QTi, see Fig. 4, Fig. 5 and Fig. 6. The total effect indices QTi provide a cumulative measure of sensitivity inclusive of interaction effects of any order. The difference QTi - Qi is a measure of the extent to which the ith imperfection is involved in interactions with any other input imperfection. The plots of QTi and Qi are similar in shape, but sizes and mutual proportions are different, see Fig. 4 and Fig. 5. In contrast to Qi, the total effects QTi show greater involvement of t2 for low slenderness and E for high slenderness. The modulus of elasticity E is a completely non-influential imperfection for 4.0 LT  (see Fig. 5) and fy is a completely non-influential imperfection for approximately 5 LT  (see Fig. 6). Fig. 4. First-order sensitivity indices Si, Qi for M = 8. Fig. 5. Total quantile contrast indices for M = 8. The plots of Qi and Si are similar in shape. However, Sobol indices Si are significantly higher than Qi, see Fig. 4. This means that QSA identifies more significant interaction effects than SSA. SSA confirms approximately zero higher-order interactions 1-iSi≈0 [35], which is in sharp contrast to the results of QSA, where the difference QTi - Qi identifies strong interaction effects. QSA confirms the results of SSA [35] in the identification of marginal main effects of imperfections h, b, t1 and  . The plots of Si in Fig. 4 (IPN 200) are almost identical to previous results of SSA of the resistance of beam IPE 220 [35]. It can be noted that beams IPN 200 and IPE 220 [35] have the same variation coefficients of imperfections h, b, t1, t2 and probability models of all other imperfections are the same in the case of steel S235 [45]. Sobol sensitivity indices Si are the same for both sections under these conditions. Fig. 6. Total quantile contrast indices for M = 8. It is apparent from Fig. 4 and Fig. 5 that the plots of sensitivity indices of steel S355 (red curves) approximately follow the plots of sensitivity indices of steel S235 (black curves), however, with the difference that the red curves are slightly shifted to the right, which is particularly noticeable in beams with intermediate slenderness. The values of indices Q5, Q7 and Q  which are lower than 0.03 and are not shown in Fig. 4. The shift to the right is due to the fact that losses (contrasts) between R and Rd of members of steel S355 are more influenced by the variability of fy and less influenced by t2 and E. Indices QTi of imperfections h and  are very small and these variables can be regarded as non-influential, see Fig. 5 and Fig. 6. The values of QTi remain practically unchanged for very high slenderness 2 LT  , with the exception of QT2 and QT4, which decrease to zero, see Fig. 6. The values of QTi and Qi are shown in Fig. 7 for the limit case of  LT  . Practically, 100 LT  was considered in the calculation. However, sensitivity indices remain constant for 10 LT  . The results in Fig. 7 apply to steel of both steel grades S235 and S355. Fig. 7. QSA indices, 100 LT  , S235 and S355. 5.2 QSA for four input imperfections The first four imperfections t2, fy, E and e0 listed in Table 1 were identified in the previous chapter as crucial. The comparison of QSA and SSA shows an agreement in the identification of non-influential imperfections h, b, t1 and  , whose variability can be completely neglected. Such imperfections can be fixed at any value of their domain without influencing the average absolute loss (contrast) between R and Rd or R and  R. Let us consider the random variability of the first four imperfections t2, fy, E and e0 (M = 4) and introduce the other four h, b, t1 and  as non-random by their mean value. This makes it possible to evaluate all 24-1 = 15 quantile indices in decomposition (8), which can be clearly illustrated; see Fig. 8 and Fig. 9.