On Selecting Composite Functions Based on Polynomials for Responses Describing Extreme Magnitudes of Structures

Featured Application: This work may be used for selection of a structural response in, for example, reliability analysis. Abstract: The main aim of this work was to investigate a numerical error in determining limit state functions, which describe the extreme magnitudes of steel structures with respect to random variables. It was assisted here by the global version of the response function method (RFM). Various approximations of trial points generated on the basis of several hundred selected reference composite functions based on polynomials were analyzed. The final goal was to find some criterion—be-tween approximation and input data—for the selection of the response function leading to relative a posteriori errors less than 1%. Unlike the classical problem of curve fitting, the accuracy of the final values of probabilistic moments was verified here as they can be used in further reliability calculations. The use of the criterion and the associated way of selecting the response function was demonstrated on the example of steel diagrid grillages. It resulted in quite high correctness in comparison with extended FEM tests.


Introduction
Diagrid structural systems have come out as some of the most efficient, most adaptable and most innovative approaches to structural buildings of this century [1,2]. Due to an esthetic potential and to the increasing architectural popularity, it seems, therefore, essential to examine their realizations and additionally compare results with those obtained for the traditional orthogonal structures. On the other hand, a variety of natural sources of uncertainties and human activity-driven reasons necessitate a series of analyses and tests to determine the impact of random parameters on the reliability of such structures.
A solution of the structural problem including randomness can be accompanied by a verification of the reliability indices calculated, for example, according to Cornell [3]or using more sophisticated indicators [4][5][6][7][8]-by using a simple limit state function g defined as the difference between structural responses fa(b) and their given thresholds fmax ( Figure 1). Then, one may apply polynomial chaos [9], Monte-Carlo simulation [10], approximated principal deformation modes [11] and/or the probability transformation method [12], for instance, to compute the basic probabilistic moments of the limit state functions. Independently of the method choice, the determination of structural responses-with respect to some input parameter subjected to uncertain fluctuations-is required. In practical cases, it is one of the crucial issues in the reliability assessment of structures, and therefore the main aim of this work since the exact dependence is usually unavailable. It can only have discrete results with some error, obtaining which takes additional time. A solution to the above difficulties is the most accurate approximation of the real dependence with an analytical metamodel [13][14][15][16] based on a properly selected set of numerical experiments [17,18]. This approach is widely used in many technical fields [19,20]. One may apply many different established metamodeling techniques: sparse polynomial chaos expansion [21,22], artificial neural network [23], multivariate adaptive regression splines [24] and radial basis function approximations [25], support vector regression [26] and Kriging models [27], for instance. The choice is a compromise between accuracy and computational efficiency.
To increase the former while maintaining the latter, the determination of structural responses was assisted here by the global version of the response function method (RFM) [28][29][30], which is quite similar to the response surface methods (RSM) applicable frequently in reliability analysis [31][32][33][34][35][36]. It is necessary to determine such a function using multiple solutions of the given boundary, the initial value problem about the expectation of the random parameter to complete this task. The resulting discrete values obtained traditionally from the finite element method experiments series [37][38][39][40][41] are in the first phase alternatively generated here from several hundred selected polynomial-based functions. They were given some predefined analytical forms and can be used as reference responses to examine various approximations of the same structural behavior; additionally, the numerical error of such modeling was under investigation here. Moreover, the generation of discrete data allows obtaining them in a very short time with any known accuracy and mathematical complexity.
Unlike the classical problem of curve fitting, the influence of these approximation choice is verified here on the final values of the probabilistic moments as they can be used in further reliability calculations ( Figure 1). Attempts were also made to find-as much as possible-marked dependences of the above value on some criterion value of the input discrete data with the proposed approximation based on curve fitting error measures. The Comparing results for the examined structures final goal was to find a limit value of such a criterion, above which the satisfactory correctness is achieved-to use it for selecting a structural response. This performance metric of approximation should also be independent of the adopted weighting scheme as much as possible and counteract the phenomenon of overfitting at the same time.
This problem was partially examined in [28,42], while the criteria proposed there were based on the results obtained using a small number of approximating curves-30 and 68 for each reference function, respectively-and only three dependences as reference response functions (additionally simplified to a linear function). Therefore, the results of the search for the response function selection criterion based on the use of 480 reference functions and 480 approximations of each of them are presented here. In [42], a generalized stochastic perturbation method was used, and now the problem is discussed in the context of the direct symbolic integration approach, which should theoretically be a more accurate method if only probabilistic moments can be recovered by integration. The integration process was provided by the computer algebra system MAPLE [43], where the coefficients of responses were also computed from several solutions of the original problem obtained for random variables varying about its expectation.
Polynomial responses are commonly applied [32,36,[44][45][46][47][48][49], and this idea follows some other applications in computational practice [50,51]. However, such approximations do not always provide satisfactory correctness [28,52]. Moreover, their application to describe the extreme magnitudes, for example, of structures seems unfortunate because polynomial limits at infinity are always infinite. In one of them, responses should tend to 0 instead as their values are only positive. Therefore, composite functions based on polynomials computed using the least squares method (LSM) were analyzed in this article. They are more flexible-due to the wider range of available functions-and can be further expanded into the Taylor series-as the compositions of elementary functions are analytical-in the perturbation technique, for instance. In turn, the different polynomial basis should help in capturing nonlinearity between the structural responses and the basic random variables.
In the last part of this work, the use of the criterion and the associated way of selecting the response function was presented on the example of eight steel diagrid grillages. The performance functions concerning the basic eigenfrequency, the extreme reduced stresses, the maximum of the global vertical displacements as well as local deflections were considered with the multiplier of (1) Young's modulus e, (2) wall thickness t and (3) length of the manufactured elements l as the truncated Gaussian random variables. Additionally, the correctness of the response function selection was verified using extended FEM tests.

Direct Symbolic Integration Approach
In the case of some real function f(b) of the stationary input random variable b having probability density function p(b), classical definitions of the basic probabilistic moments and coefficients were introduced [53-56], for example: • the expectation: the variance: • the m th central probabilistic moment (for m > 2): The Gaussian probability density function is considered further, so that: where σ(b) ≡ σ stands for the standard deviation of b substituted further by a product of the expected value The integral definitions proposed in Equations (1)-(3) are very rarely used with infinite limits. Random variables are usually analyzed as truncated having some lower and upper bounds of the probability density function driven by the physical meaning of the specific parameter or just the experimental works. Therefore, the integration process is limited here using the well-known three-sigma rule recognized as very efficient in various computational experiments [29,57]. Once the response function represents a reference exact dependence f(b), we obtain a solution of the probabilistic analytical method (AM), but when it is only an approximation fa(b)-of the semi-analytical method (SAM) [30,58].
Due to a large number of symbolic calculations in Section 3, in order to reduce their computation time, the semi-analytical method was approximated by determining the integrals using a numerical quadrature technique-the midpoint Riemann sum method [59]. The number of its subintervals was selected experimentally (see Section 3.4).
Probabilistic characteristics calculated on the basis of the formulas provided in this subsection can only be discrete. Their values were determined each time for the input coefficients of variation αi ∈ {0.0125, 0.025, ..., αmax} with an increment of Δα = 1/80, where for αmax = 0.30, there were 24 calculation points.
In this work, the results of the probabilistic analytical method were taken as exact values. The approximate solutions of the semi-analytical method were compared to them in order to test the criteria for selecting the response function.

Response Function Method
The response function method consists in approximating the real dependence of the structure response f(b) with an analytical metamodel-the response function fa(b). To complete this task, it is necessary to determine such a function using a multiple solution of an investigated boundary value problem about the mean value of the random parameter, for bi (i = 1,...,N) belonging to the interval b 0 − Δb,b 0 + Δb. As a result, we obtain a set of discrete values fo (bi) with an error in the numerical procedure, e.g., the finite element method.
There are several ways to choose the Δb/b 0 ratio [57,60] as well as this interval's discretization both in terms of uniformity and the number of trial points [29,61]. Here, uniform interval subdivision with n = 11 trial points and the ratio Δb/b 0 = 0.05 was applied ( Figure 2) as the most frequently used. This choice was confirmed using the numerical experiments included in [28,58].
Each unknown response function was approximated here using composite functions based on polynomials (Table 1) computed using the classic least squares method [62] as well as its weighted version (WLSM) to modulate the importance of computational analysis results [63]. In the latter method, the values computed for the expectation of random variables were recognized as crucial. Their weights (10) were assumed as the sum of the rest ten weights of equivalent results (equal to 1) [57,64,65]-the Dirac-type distribution of the weights. However, in connection with the conclusions contained in [28], a uniform weight scheme was adopted as the basic one. The weighed one is only auxiliary. When considering many approximations fa(b), the selection of the final response function in the initial stage of RFM development was made arbitrarily [30,66], while now it is the subject of additional optimization [57,65,67,68] related to curve fitting error measures (see Section 2.3): correlation maximization, minimization of the root-meansquare error (RMSE) or the residues variance. In addition to or interchangeably with the criteria calculated between the discrete data and the approximation, analogous magnitudes calculated in-between the LSM and the WLSM solutions are used [28,42]the weights change in calculations does not affect how the structure behaves in reality. If the selected criteria indicate different approximations, it is assumed that the one with the smaller coefficients number is chosen [64,67], which reduces the risk of the Runge effect [69]-deterioration of the quality of polynomial interpolation, especially visible at the ends of the intervals, despite increasing their degree-and is consistent with the principle of using the simplest possible theory called Occam's razor. It is widely used, among others, to choose between statistical models with a different number of parameters: in the AIC [70], BIC [71], SABIC [72] or RMSEA criteria [73].

Curve Fitting Error Measures
The first curve fitting error measure is the linear correlation coefficient. Let fo(bi) and fa(bi) denote the obtained discrete data and the values of the considered approximations, respectively, for the given bi (i = 1,...,N). Then, the following formula applies: where the general formula for the sample mean has the following form: The closer the value of the coefficient ρ is to 1, the more the considered data series are correlated and the approximation courses are closer to the trial points. Due to the above and the narrow range of values (−1 ≤ ρ ≤ 1), greater visual differences are obtained when considering the value -log10 (1 -|ρ|).
Another measure of approximation fitting is the root-mean-square error: where the magnitude ri is called the residuum, error or remainder. The next measure, fitting variance, is also based on this value: The coefficient of determination was also considered with the given formula: while due to its range of values (0 ≤ R 2 ≤ 1), greater visual differences were obtained considering the quantity −log10 | R 2 -1|.

Assumptions and Selection of Reference Functions for Numerical Examples
It was decided to consider here both the reference dependences and approximations having the following polynomial form: where it was assumed that the quantities used above are functions of single variables, e.g., in order to use the LSM method: and This leads to the following dependence: where the function w is hereinafter called the base polynomial and the argument is a random variable x with the truncated Gaussian probability density distribution, which is characterized by the input coefficient of variation α(x)∈[0.0, 0.30] to explore a fairly wide three-sigma range (0. The use of as mathematically simple as possible selected elementary functions used in empirical formulas in Equations (11) and (12) led to obtaining 48 groups of dependences based on polynomials. Their list is included in Table 1.   Ultimately, they are to constitute performance functions describing extreme magnitudes of structures, e.g., deflections, displacements or stresses, with respect to a random parameter. This means that they only have to take positive values, which is the first limitation when it comes to choosing the reference base polynomials for numerical examples: x E x (15) and according to the three-sigma rule. It should be emphasized that the sought reference base polynomials are assumed to be used for each dependence from Table 1, therefore, taking into account all forms of the functions, Inequality (14) is reduced to the following: which, due to the properties of the considered functions φ(x), is additionally simplified in the case of x > 0: and Finding the exemplary coefficients for which Inequality (18) is satisfied was facilitated by narrowing down the interval (Xmin, Xmax) in which the base polynomial was considered. Its smallest width occurs at E[x] ≈ 1.14, while in further calculations E[x] = x 0 = 1.0 was assumed. With such a value, the variable x is considered in the interval (0.1, 1.9), which satisfies the constraint imposed by Inequality (18) so that x > 0.
In order to be able to also easily assess the accuracy of the approximation visually, all coefficients values of the sought base polynomials are assumed to have the form |Di| = 10 k , where k ∈  and i ∈ {0,...,n}. Moreover, their value should be such that after substituting for all groups of functions, it should be possible to further use the generated discrete data in the computer algebra system so that they do not exceed the limits imposed on the numbers, which was not easy to achieve, especially for the dependency y30.
It was decided to analyze the base reference polynomials up to the 6 th degree inclusive in order to obtain the appropriate level of their "complexity" and at the same time leave the possibility of approximation using the functions based on the base polynomials of much higher orders-up to the 10 th degree inclusive. Additionally, the possibility of obtaining a polynomial in which the influence of components derived from higher powers of the variable would clearly decrease-the polynomials only theoretically of a given degree-was avoided. The assumption was made that it should be similar to what can be schematically written as follows: where j ∈ {1,...,n − 1}. Its satisfactory approximation is obtained-for E[x] = 1.0-when |Dj + 1| ≈ |Dj| which, in the presence of the previous assumptions, assumes the following: The value of the coefficient D0 was obtained by considering the function After substituting the data for Inequality (17), it was observed that D0 > 1; therefore, in the presence of the earlier assumption, the smallest possible coefficient was adopted: D0 = 10.
Various values of the coefficient k occurring in Equation (22) and all the possible sign configurations in the D1, ... ,Dn coefficients were analyzed. For numerical examples, four polynomials were finally selected for each degree-only two for the first one because a larger number was not possible-for which Inequality (18) was satisfied and the maximum value of the function y = ψ −1 (w(X)) in the interval (Xmin, Xmax) was the smallest. In the case of function w, the degree of which is not greater than 4, it was also possible for half of the dependences to select them for the two different smallest values of the coefficient k. The list of the selected reference base polynomials is given in Table 2.
Step One-Accuracy of the Coefficients In the first part of the experiment, the influence of the discrete data accuracy on the coefficient values D * i of the approximation polynomials w* of the order POa was analyzed. For this purpose, cases were considered where discrete data were generated on the basis of a given group of functions as well as a reference base polynomial w and then approximated by the same type function. Due to the assumption of a different a priori error, by rounding the values of the trial points to a certain number of significant digits, it was then possible to analyze its impact on the a posteriori error assumed as the maximum from the value of the relative error (more precisely, its module-this detail is omitted in the further part of this work to shorten the descriptions). It was calculated for individual coefficients of the approximation polynomial in relation to the reference polynomial. This part of the experiment was shown schematically in Figure 3. All the groups of functions from Table 1 as well as all the base reference polynomials from Table 2   In the case of polynomial dependences (dashed lines in Figure 4), the relative error value was initially 1.0 regardless of the degree of the polynomial under consideration. This is due to the fact that the generated discrete data, when rounded to a very small number of significant digits, had identical values. The approximating function was then constant, the values of the coefficients at the nonzero powers of the variable were equal to 0, so their relative error in relation to the respective nonzero coefficients of the reference polynomials was 100%. At a slightly larger number of significant digits of rounding, the individual discrete data generated started to differ from the others, but only in the last place of their decimal notation. An attempt to approximate these values-the diagram of which is clearly a polygonal chain-with smooth polynomials caused maxΔrel,c to start exceeding 100%except for linear functions. Only when we achieved a sufficiently high accuracy of the values of the test points-which guaranteed their greater differences-it only started to decrease.
When considering all the other groups of functions (solid lines in Figure 4b), a greater precision of discrete data was required to obtain the same a posteriori error than for the polynomial dependences (dashed lines). This is equivalent to the fact that with the established accuracy of discrete data, in the case of the polynomial dependences, a smaller relative error of the coefficients maxΔrel,c was obtained than for the other considered functions (see the example in Table 3). This is due to the fact that in the case of the first group of functions, it was possible to directly calculate the coefficients in Equation (10), while in the case of the others, an additional modification of the arguments or the set of values was required in accordance with Equations (11) and (12), correspondingly.
Regardless of the choice of the functions group, with the exception of the initial perturbations, in the case of reference polynomials of higher degrees, clearly greater precision of discrete data was required to obtain a specific error rate of the coefficients. It is worth noting that the not too high critical a posteriori error value of 5% was reached approximately for the second-order base polynomial and seven significant digits of discrete data rounding (Table 3, Figure 4). The standard accuracy of the FEM results is much lower, so if they are used to obtain the structure response function, even for the base polynomials more complex than linear, there is a risk of exceeding the above critical value of the relative error of the coefficients maxΔrel,c. It is also worth emphasizing that thus far, it had been assumed that discrete data were described by a specific group of functions based on a polynomial w with a known degree, while in most practical cases, the form of the response function is unknown. Therefore, the response function with satisfactory accuracy of the coefficients can be obtained on the basis of the FEM results only in a very limited number of simple cases in which we know that the state function dependence is linear-in some cases it may be enough to reduce to it. This situation is a serious problem that could prevent further effective application of this method in practice. However, when the RFM is used for reliability analysis, more important than the accuracy of the function coefficients is the one connected with the probabilistic moments because they can be used in further calculations ( Figure 1).

Step Two-Accuracy of the Probabilistic Moments
In the formula for the Cornell reliability index, there are only the expected value and the variance. However, the description of the distribution only with their use is not sufficient because the values of these parameters may be similar for extremely different data, as in the case of the Anscombe quartet [74] (Figure 5). Taking into account the additional values of the third and fourth central moments-and the skewness and kurtosis based on themallowed indicating numerical differences between such data (Table 4). It should also be noted that reliability can be ultimately described by indicators whose formula takes into account not only the expected value and the variance [4][5][6][7][8]. This is another confirmation of the statement that the compliance of values should also apply to some additional parameters describing the distribution.  In the further part of the experiment, four probabilistic moments were considered: the first raw moment-marked as m1-equivalent to the expected value is given in Equation (1), while the central moments from the second upwards-marked from m2 to m4-are described in Equation (3). The approximation of the discrete data to D1 = 7 significant digits was adopted because it is sufficient for groups of functions based on the linear base polynomials (maxΔrel,c below 1‰), not sufficient for the second-order base polynomials (maxΔrel,c even above 100%, Table 5) and definitely insufficient for the fifth-order base polynomials (maxΔrel,c even above 5 × 10 8 ). Therefore, in the further part of the experiment, only the results obtained on the basis of base polynomials with the abovementioned degrees were considered to verify whether, regardless of the possible accuracy of the coefficient values-which is also the resultant of the significant digits number-it is possible to obtain satisfactory accuracy of probabilistic moments measured by their relative error Δrel,m. The following expression was employed on the example of the expected value: The initially obtained results confirmed that for different approximations, the shape of the diagrams of the probabilistic moments' relative error Δrel,m may also differ ( Figure  6). The higher error value was not always obtained for the higher moment. Its monotony was not always only increasing in the considered range. Therefore, it was decided to use the maximum of all the values in the interval (0.0, αmax*) as a measure of accuracy of the given approximation: The results obtained in the validation process confirmed its compliance (Figure 7). At the same time, they indicated a more general, assumed character: for some approximations, despite the significant relative error of the coefficients, the one concerning probabilistic moments was clearly smaller. One of the greatest advantages of the proposed measure of accuracy maxΔrel,m is the possibility of comparing any approximation to a given reference function despite having different mathematical forms, which was not possible in the case of the previous accuracy criterion Δrel,c.

Step Three-No Information about the Form of the Function the Discrete Values Come From
Having the new approximation accuracy criterion finally made it possible to analyze cases where no assumptions are made about the function approximating the discrete data-both the function groups and the degree of the base polynomial. However, knowing the form of the function on the basis of which the trial points had been generated, it was possible to investigate the a posteriori error for individual approximations.
In this part of the experiment, 480 reference functions were analyzed-all the 48 groups of functions were obtained on the basis of all the 10 reference polynomials-and 480 approximations of each of them-48 groups of functions with polynomials from the first to the tenth degree. Temporarily, this gave over 230,000 approximations, with this number then dropping to fewer than 175k due to the fact that the domains of the approximations did not coincide with the examined interval of x. This part of the experiment is shown schematically in Figure 8.  The first stage of the above scheme consists of determining for each approximation the relative error values for four moments and 24 discrete values of the input coefficient of variation ( Figure 6). In total, this gives almost 17 million values of probabilistic characteristics to calculate, which would be very time-consuming in the case of symbolic integration. Therefore, when determining the values for the approximation, it was decided to make approximate calculations using the midpoint Riemann sum method [59].
Calculation of the maximum relative error of the probabilistic moments maxΔ rel,m for each approximation Calculation for each approximation of the criteria value (linear correlation coefficient, root-mean-square error, fitting variance, coefficient of determination) Criterion choice

Determination of the criterion limit value
The diagrams in this part of the experiment were partially analyzed visually and the most important part of the results was related to the 1% error, so the analysis of the results focused on the range from 10 −4 to 10 0 -two orders of magnitude difference. In order to determine how large the relative error of determining the integrals for the approximation was acceptable, a line diagram with a vertical logarithmic scale and its different spreads Δappr was considered. It simulated the relative error of a certain probabilistic moment of exemplary approximation with respect to the reference value (Δrel,m). A satisfactory accuracy was obtained with the value of Δappr = 10 −5 (Figure 9). Due to the nature of the diagrams, the upper limit for considering the error Δrel,m could be increased to 10 1 without the loss of accuracy. In order to be able to present all the results in the diagrams later in this section, the upper limit was modified to represent values with the relative error greater than or equal to 10 1 , and the lower limit-those with the relative error less than or equal to 10 −4 .
(a) (b) The initial tests were carried out by analyzing Δrel,a-the relative error of determination of the SAM integrals using the approximate method in relation to the analytical method: m,AM m,approxSAM rel,a m,AM The tests were performed with a changing number of subintervals for different configurations of discrete data, approximating functions and values of the coefficient of variation α. Each time, decreasing dependences were obtained, while in the most unfavorable case, to obtain the value of Δrel,a = 10 −5 , it was necessary to divide the integration interval into 5000 equally wide subintervals ( Figure 10). This value was used for determining all the probabilistic moments for approximation with the approximate method, which was the case in the first stage of the calculation scheme presented in Figure  8.
After the calculation of the curve fitting error measures for each approximation, according to Equations (5)-(9), the third stage was started. After a long search, the most marked dependence of the maximum value of the a posteriori error of probabilistic moments on the criterion value was noted for a modified root-mean-square error ( Figure  11a) provided in the following formula: Figure 10. The relative error in determining the integrals for an approximation using the approximate method depending on the number of subintervals (the most unfavorable case of the received ones).
It assumes a weighted product where RMSE is calculated between the LSM approximation and the input discrete data, while RMSEw-between the LSM approximation and its weighted version (weights 1,1,1,1,1,10,1,1,1,1,1). It makes the result independent of the adopted weighting scheme and thus makes the solution more probable. Such a situation should occur due to the lack of influence of the weights adopted in the calculations on the actual behavior of the structure. In order to increase the readability of the graph, the logarithmic scale was also used. Moreover, it was required to use the so-called penalty term, i.e., an expression that prefers approximations with the lowest possible degree of the approximation POa. This factor counteracts the phenomenon of overfitting when the approximation has too many parameters in relation to the sample size and practically passes through all the points having a complicated shape between them. The value of the constant coefficients in Equation (28) was determined by using the trial and error method and analyzing the results obtained for αmax*(x) = 0.30 because earlier for such a value the criteria for the maximum relative error of the coefficients and moments turned out to be the closest (Figure 7). The RMSEmod criterion made it possible to obtain satisfactory accuracy of the probabilistic moments maxΔrel,m for the discrete data generated by the functions based on polynomials of various POb regardless of the possible accuracy of the values of the coefficients Δrel,c (Table 5) being the result of the number of significant digits rounding of the trial points. The validity of using the criterion was also confirmed by the graphs showing the indicated approximations in relation to the reference functions and the discrete data generated (Figure 12).
The RMSEmod rounded value above which-for αmax*(x) = 0.30-all the obtained a posteriori errors are below 1% is 36 (Figure 11a). It is worth noting that this number was almost independent-it visually imperceptibly decreased-from the interval width of the input coefficient of variation in which the accuracy of the approximation was testedequivalent to αmax*. With its decrease, the relative error value maxΔrel,m decreased ( Figure  11a-f). It should be noted that the least squares method was also based on the condition of minimizing the criterion similarly to the RMSE. In the opinion of the author, the abovementioned conclusions confirm that the proposed criterion is correct.  The limit value of the adopted criterion in relation to the form of the penalty term resulted in that only the approximations based on the base polynomial of the maximum eighth degree could exceed it. In the case of the ninth degree, despite the minimum values of the RMSE and RMSEw components (but not less than 10 −20 -the number of digits used in the experiment by the MAPLE system when performing calculations using the floatingpoint numbers was 20), the final value of the RMSEmod coefficient was still lower than the limit value equal to 36. This limitation seems rational if we take into account the number of trial points (n = 11) and the desire to avoid the overfitting phenomenon.
It is worth mentioning that all the considered criteria determining the quality of approximation fit-linear correlation coefficient, root-mean-square error, fitting variance, coefficient of determination-used both individually and in pairs did not lead to a marked dependence of the maximum value of the a posteriori error on the criterion value. Adoption of the approximation for which their extreme value is achieved in most cases also led to results of unsatisfactory accuracy-points with the ordinate above 10 −2 in Figure 13. In the abovementioned experiment, a total of 174,469 approximations were analyzed, and 17,660 of them (10.1%) had error values less than 1%. Only 1003 of them (7.91% from the previous group) had the RMSEmod value exceeding the established limit. This was due to the desire to ensure-apart from the required accuracy-that the result was independent of the adopted weighting scheme and counteract the phenomenon of overfitting. The approximation that would meet the criterion limit value was found not for all the sets of trial points in the experiment. Often, none of the 480 approximations for a given series of discrete data met the error condition-points with an ordinate greater than 10 −2 in Figure 14a. It could be caused by (1) the insufficient number of analyzed functions, (2) the insufficient number of trial points or (3) too narrow interval of their examination.
The first reason is justified by considering only polynomial approximation. The best results obtained ( Figure 14b) were significantly worse than without this restriction ( Figure  14a). Additionally, in most cases, the assumed threshold of 1% of the relative error maxΔrel,m was exceeded (Figure 14b). The second and third reasons are justified by the results contained in [28], where more accurate results were obtained when considering a larger number of discrete data for a wider range of the input variable.
Bearing in mind the application of the established criterion in practice, it should therefore be emphasized that the adopted base of approximation in Table 1 may turn out to be insufficient. In such a case, it is recommended to extend it until a satisfactory dependence is found. Increasing the number of trial points or the width of their consideration interval may then require modification of the criterion, which necessitates additional experiments. Alternatively, it is also possible to select approximations from among those maximizing the value of RMSEmod, while considering them in the full range of variability of a given random parameter and rejecting solutions that differ from the others. This way is presented in the next section.

Finite Element Analysis
In this experiment, eight grillages designed in accordance with the applicable standards for steel structures [75][76][77] were analyzed. They were computational models of a supporting steel structure for a rectangular glass floor (9.00 × 5.20 m). It belongs to the category of use C1 [76], including areas with tables where people may congregate, e.g., cafés or restaurants. It was additionally assumed that due to the visual (architect or investor) and design requirements (easier connections, as well as the internal forces transfer), the dimension of the cross-section in the normal direction must be constant within each model.
The following were assumed as random variables with the truncated Gaussian probability density distribution: the multiplier of (1) Young's modulus e, (2) wall thickness t and (3) length of the manufactured elements l. They represent the variability related to the material, cross-sections and geometry, correspondingly. The value of their coefficients of variation was considered in the interval (0.0, 0.25). All the expected values were 1.0 since e, t and l are multipliers of deterministic quantities. The whole numerical example concerning steel diagrid grillages is shown schematically in Figure 15. The FEM analysis was carried out with the use of the civil engineering system ROBOT. Its usage was driven by the fact that this system enables both obtaining discrete solutions of a given structural problem as well as shaping and dimensioning steel structures so that the analysis is fully supported by this program.
Two types of grillages were considered: orthogonal (model OB and OS) and diagrid. The latter were additionally grouped into having a mesh of right triangles (models RB and RS) and equilateral triangles arranged in the transverse (models TB and TS) as well as longitudinal directions (models LB and LS) to represent different possible architectural concepts of floor triangularization. In each of these groups, there were versions with a smaller and bigger mesh size ( Figure 16), with uniform division, so that the nodes of the real structure are as repeatable as possible. All of the supports were defined as rigid with fixed movement and free rotation in each direction. Table 6 shows the properties of the models in detail. There, it can be seen how the reduction of the mesh size directly influences the increase in the number of nodes, supports, structural members and panels. However, a thinner layer of laminated glass glazing is required. In this case, it consists of two hardened load-bearing panels and an 8 Selection of the response function according to the RMSE mod criterion mm thick anti-slip top layer-also hardened. The former were dimensioned each time according to the appropriate design standards [78][79][80].
All the loads from Eurocodes in the persistent design situation were considered for the design working life of 50 years-a common structure according to [75]. The first (marked with G) is connected with the weight of the glass covering and the supporting structure itself. The second (marked with Q) represents the imposed live loads of 3.0 kN/m 2 [76]. The surface loads were implemented. They were distributed using the trapezoidal and triangular method.
The 3D frame two-noded finite elements with six degrees of freedom in each node and rigid connection were used in each model. Based on the incremental formulation of the FEM equations [39,81,82], the following approximation of the displacement increment vector Δue (ξ) was considered in the finite element: (29) where N represents the shape functions matrix, Δqe denotes the increments of the displacements vector and ξ is the dimensionless coordinate of the section location.
The global deformation was computed in the serviceability limit state (SLS) combination (G + Q), the stresses-in the ultimate limit state (ULS) combination (1.15 × 0.85G + 1.5Q) obtained in accordance with the formulas contained in [75]. The secondorder global analysis P-Delta was used taking into account the influence of deformation on the statics of the system-geometric nonlinearity-to more realistically predict the structural behavior [77]. The incremental method was applied using a modified Newton-Raphson algorithm to solve the nonlinear problem. Therefore, the following FEM matrix equilibrium equations series [83] were used for each test necessary for the response function method recovery: where K denotes the stiffness matrix of the system changing under the influence of deformation, resulting classically from an aggregation process over all the finite elements, q is the generalized displacements vector, Δ denotes the increment of a given quantity, Q is equivalent to the external loads vector, n denotes the increment number and superscript i indicates the iteration number. The stiffness matrix K was updated only after each subdivision-not after each iteration-as a result of its correction with the algorithm of the BFGS method. The following parameters were set to complete these experiments: The finite element method-based modal analysis was performed, too. Higher modes were also computed, but only the first mode was verified. The subspace iteration algorithm was used for solving the eigenvalue problem. The maximum number of iterations was set as 40 with a tolerance factor of 0.0001.
Yield strength was assumed as 235 MPa and the structural members were made of rectangular hollow sections ( Figure 16). This type of cross-section was chosen because of its insensitivity with regard to the lateral-torsional buckling, the creation by their top walls of a flat surface to support the panels; hence, it is usually used in structures with glass covers [84]. This choice made it possible to avoid considering the warping in the FEM analysis as an additional degree of freedom in the node (the seventh in this case) because its influence is negligible for the selected cross-sections.
In each model, the sizing process was similar and deterministic, carried out entirely according to Eurocode 3-1-1 [77]. Initially, an identical RHS section of a given dimension in the normal direction was adopted for all the bars. Next, a division into groups of bars was made, taking into account the distribution of internal forces in the ULS combinations and deformations in the SLS combinations. For each of them, a cross-section with a given dimension in the normal direction was then calculated. The process of creating groups and sizing their cross-sections was iterative. Its limitations were finding no more than three profiles (for a given grillage model-economic reasons) available in Poland and their distribution in the structure, for which all the conditions of the limit states were met for no more than 90%. The above procedure was carried out for various constant dimensions in the normal direction-e.g., 250, 260, 300, 350, 400 mm-finally choosing one for which the lowest mass was obtained. The final profiles are shown in Figure 16, and the massescalculated based on the lengths of the model bars, not the real structural elements-in Table 6. The minimum steel weight was found for grillages with the orthogonal mesh (models OB and OS), the maximum-for the triangular mesh in the transverse direction (models TB and TS). The final purpose of the finite element method experiments was to obtain discrete solutions of the given structural problems for all the grillages. Therefore, the resulting values of the fundamental frequency, the extreme reduced stresses (according to the Huber-Mises-Hencky hypothesis), the maximum values of the global vertical displacement and local deflection were determined.

Using the RMSEmod Criterion
For all the FEM data series, each of the 480 previously mentioned approximations were adopted. However, those whose domains did not coincide with the examined interval of the random variables were rejected. In four out of eleven cases, the average of the RMSEmod criterion values from all the models exceeded the limit value of 36 (Table 7). In the remaining cases, several approximations with the maximum RMSEmod values were compared visually. In the range of the FEM tests, the difference was imperceptible ( Figure  17a). However, when considering them in the full range of the assumed variability of the random parameter, the diagrams of some of them significantly differed from one common course for the majority (Figure 17b). Such cases occurred even for the approximation with the maximum mean of the RMSEmod values among those obtained for the given series of discrete data (Table 7). Therefore, all of them were omitted when selecting the response function in favor of the next ones sorted in the descending order according to the criterion value. Thus, the diagram of the dependence selected in this way-in Table 7, the corresponding RMSEmod value is marked with an asterisk-coincided with those obtained for many unselected approximations having significantly different forms of formulas. In the authors' opinion, this indicates correctness of the solution.  The forms of some of the final response function formulas additionally confirmed correctness of the response function selection according to the RMSEmod criterion. When considering the variability of Young's modulus, we obtained inverse proportionality to the first power for deflections and dependence on the root of e for the basic eigenfrequency ω. This is justified by well-known dependences of the linear elastic analysis: • for elastic displacements (with the assumption of small deformations): • and for the first eigenfrequency: where C is a constant depending on the parameters of the structure and the considered performance function. When the random parameter is t, both the state function ug and ω are based on the fourth-order polynomial, which occurs in the analytical description of the inertia moment Jrect for rectangular hollow sections as a function of their walls thickness t0: where bw is the internal dimension of the shorter side and hw is the internal dimension of the longer side. However, it should be mentioned that the assumption of constant internal dimensions and simplification that ignores the occurrence of fillets were used in the example. The correctness of the response function selection was also confirmed by the course of their diagrams (Figure 18), which is compliant with the engineering intuition in this matter. With the increase in wall thickness multiplier t, the deflections and reduced stresses nonlinearly decrease, and the fundamental frequency goes up slightly. However, when the length of the elements increases, an even faster gain is observed for the state functions ug, ul and σred, but a decrease for ω. As far as the quantitative assessment is concerned, it should be emphasized that the abovementioned functions assume elasticity in the full range of variable uncertainties, so that the theoretical values of the reduced stresses sometimes exceed the adopted yield strength fy = 235 MPa (Figure 18d). This limitation should be taken into account as a limiting condition in the further reliability assessment, so the values in the part of the σred diagram above fy will not matter.
The correctness of the response function selection was confirmed by the extended FEM tests, too ( Figure 18

Concluding Remarks
The method of creating groups of composite response functions describing the extreme magnitudes of structures was applied. They are based on polynomials, the coefficients of which can be obtained using the LSM. The functions type and the number of significant digits of trial points rounding had a clear impact on the a posteriori error of the approximation coefficients. When a sufficiently high accuracy of the input discrete data was achieved, it only started to decrease. At the same time, a smaller relative error of the coefficients was obtained for the polynomial dependences than for the composite functions. In order to maintain it, in the case of the base polynomials of higher degrees, a greater number of significant digits of discrete values rounding was required. Finally, it was noticed that even knowing the form of the response function, with the standard accuracy of FEM results, for dependences based on a polynomial more complex than linear, there is a risk of exceeding the 5% relative error of the coefficient values.
An a posteriori correctness measure of approximation was proposed. It is the maximum of the relative errors of four probabilistic moments obtained for the variance coefficients in the interval (0.0, αmax*). Its most marked dependence on the criterion value-between the approximation and the input data-was noted for the modified root-mean-square error. It makes the result independent of the adopted weighting scheme, at the same time counteracting the phenomenon of overfitting. Its rounded value, above which the satisfactory correctness was achieved, equals 36. If this condition was not met, satisfactory results were obtained by selecting approximations from among those maximizing the value of RMSEmod, while considering them in the full range of variability of the random parameter and rejecting solutions that differed from the others. This was confirmed by the forms of the simplified formulas, the extended FEM tests and shapes of diagrams of the response functions thus selected for the grillages.
The analyzed responses had a meaningful influence on the probabilistic moments. Therefore, if the form of the response function is unavailable, it is suggested to search for it according to the proposed criterion and the given way related with it. The adoption of approximation on the basis of the previously commonly used criteria-e.g., fitting correlation maximization or fitting variance-did not lead to results with satisfactory accuracy.
This study may be further extended to analysis of some random variables not exhibiting Gaussian probability density functions when precise determination of higherorder statistics is of paramount importance in computational mechanics [85].