Modiﬁed, Reliability-Based Robust Geotechnical Design Method, in Accordance with Eurocode 7

: In this paper a modiﬁcation of the reliability-based robust geotechnical design (RGD) method is proposed. The intention of the proposed modiﬁcations is to simplify the method, make it less computationally expensive, and harmonise of the results with Eurocode 7. The complexity of the RGD method mainly stems from the calculation of the design’s robustness measure, which is the feasibility robustness index ( β β ). Due to this fact, the replacing of the existing robustness measure with a generalised reliability index ( β ) is considered. It was demonstrated that β ﬁts into the robustness concept, and is traditionally used as a construction reliability measure, making it intuitive and “user friendly”. It is proposed to conduct a sensitivity analysis using Soboli indices, with the aim of freezing the variables whose contribution to the system response variance is negligible, which will further simplify the method. By changing the robustness measure, the number of the required reliability analyses is signiﬁcantly decreased. Further reduction is achieved by conducting analyses only for the designs chosen in the scope of the genetic algorithm. The original RGD method is used as an extension of traditional reliability-based design. By applying the proposed modiﬁcations, the RGD method can be used as an alternative to the classic and reliability-based design method.


Introduction
A new geotechnical design approach, called "reliability-based robust geotechnical design method" (RGD) was suggested by Juang et al. [1] in 2012. The term "reliabilitybased" stems from the fact that the aim of the method is "to ensure robustness of reliabilitybased design even if uncertainties exist in the estimated statistical moments of noise factors". They prescribed the methodology of defining non-dominated and robust designs of shallow foundations in cohesionless soil. The aim of the RGD method is to calculate a set of non-dominated designs (Pareto front) through multi-objective optimisation, which is conducted using the genetic algorithm NSGA-II. The Pareto front is composed of robust designs that meet the prescribed safety requirements. The final design is chosen from this set, and it is a trade-off between the cost and the robustness of the foundation.
The concept of robustness dates from the 1950s, when Genichi Taguchi set the basic principles of the robust method, based on research in Quality engineering [2]. In this context, a robust product is defined as "a product whose performance is minimally sensitive to factors causing variability (at the lowest possible cost)" [3]. The application of the robustness concept to the problems in structural and geotechnical engineering has resulted in an adjustment of terminology and methodology, and an array of different definitions connected to the usage of the term itself. A review of relevant literature has shown differing approaches to the structural robustness, with some of them being fundamentally different.
In the RGD method, a design is considered robust if the system response is insensitive to the variation of noise factors and the correlations between them [1].
The RGD procedure is generally complex and computationally expensive. Different authors therefore suggest its modifications, most commonly referring to the robustness measure. Tan et al. [4] keep the existing robustness measure (feasibility robustness index β β ), but suggest its calculation using the sensitivity of reliability index (SRI) incorporated with FORM. The original method uses FORM with integrated PEM. In this case, the total number of analyses is n = 7Mn c , where M is the number of possible designs, and n c is the number of noise factors. The modified method has an M number of analyses, which accelerates the calculation significantly.
Gong et al. [5] presume that the existing robustness measure is not fundamentally sound and intuitive, and that it is often isolated from the analysis of design constraints. Consequently, they suggest a gradient-based robust design approach, in which noise factors are treated as random variables, and the design constraints are analysed using advanced First Order Second Moment (AFOSM) method. The design robustness measure is a gradient-based sensitivity index (SI), calculated within probabilistic analyses. Treating noise factors as random variables eliminates standard deviations of their own standard deviations from the procedure, which, in turn, simplifies input parameters significantly. The total number of reliability analyses is reduced to the value of M, which improves the computational efficiency of the method. The sensitivity index, as the design robustness measure in the RGD of soldier-pile-anchor shoring system for deep excavation, is also suggested by Gong et al. [6]. They simplify RGD by reducing multi-objective optimisation to a sequence of single-objective problems. Yu et al. [7] apply the RGD for slope stabilisation using piles. They use signal to noise ratio (SNR) for the design's robustness measure, which is calculated from the standard deviation and the mean value of the factor of safety (FOS). The same robustness measure is applied by Zhou et al. [8] in geotechnical slope design. Fan et al. [9] define a new robustness measure for application in the design of rock wedge slopes-the sensitivity index of variability (SIV). The main shortcoming of SIV, considering its more general application in geotechnics, is the assumption that all random variables are independent and normally distributed, which is often not the case in practice. Ravichandran and Shrestha [10] suggest a robust procedure of optimisation for a typical foundation for wind turbines. The robustness of the design is measured using a standard deviation of differential settlement and calculated using the Monte Carlo method. The Monte Carlo method is generally simpler, but in a scenario with a large number of random variables and complex objective functions, it becomes computationally expensive. Peng et al. [11] simplify the RGD procedure by reducing the problem of multiobjective optimisation to a sequence of single-objective optimisations, and by substituting the FORM method with integrated PEM with the Monte Carlo method. This makes the RGD procedure more intuitive and more user friendly, while maintaining the robustness measure. An example of a rock slope design is given. Since the objective functions in the example are mathematically simple, the suggested method is more efficient than the original one. The authors suggest an assessment of the method's efficiency in examples with more complex objective functions and a greater number of random variables. The summary of reviewed literature is presented in Table 1.
In this paper the RGD method is simplified, made less computationally expensive, and its results are harmonised with the criteria prescribed in the Eurocode 7. The simplification of the method's application stems from the change of the robustness measure, where a generalised reliability index (β) is chosen as a substitute of β β . It has been used for years in structural reliability analyses to express construction reliability. Due to this fact, its value can be considered intuitive. Furthermore, the values are within the approximate 0-20 range, which means that it can also be considered "user friendly". In this paper we demonstrated that β, by definition, fits into the robustness concept. Based on the chosen robustness measure, the required number of reliability analyses has been reduced. A modified PEM method (IPEM) has been suggested for the calculation of β, which enables a use of simpler mathematical operations while maintaining a high degree of accuracy. The total number of reliability analyses is further reduced by conducting them only for the designs chosen in the process of optimisation, using a genetic algorithm. Additional safety constraints are introduced into the optimisation process, to adjust the designs to the ultimate limit state (ULS) and serviceability limit state (SLS) criteria, as defined by Eurocode 7.
The suggested modifications can be applied to geotechnical problems in general; however, for illustration purposes, the methodology in this paper is elaborated using the example of shallow foundations in cohesionless soils.

Original Method-Reliability-Based, Robust and Optimal Design (RGD) Method of Shallow Foundations in Cohesionless Soils
The procedure of applying the RGD method as defined by Juang et al. [1] is shown in Figure 1. The first step requires a characterisation of noise factors, decision variables and the design domain. Since there are uncertainties in the second statistic moment of geotechnical random variables as well, they are considered to be "noise factors", i.e., "hard-to-control" parameters. Consequently, their standard deviations are treated as random variables, i.e., data on their statistical distributions, mean values and standard deviations is required. Since engineering practice often lacks data for a quality statistical analysis, an estimation is suggested using the statistic bootstrap method. Unlike the statistic characteristics of geotechnical parameters, literature holds no information on typical statistic characteristics of their standard deviations, and their estimation via the bootstrap method is not reliable. This introduces additional uncertainties into the method, which impacts the final results, attained conducting reliability analyses. The middle part of the procedure features an inner and outer loop, which include the calculation of the design robustness measures using the FORM (First Order Reliability Method) with integrated PEM (Point Estimate Method). Design robustness is measured by the index od feasibility robustness (β β ). Determining β β requires the information on the statistical distribution of β. The authors of the RGD procedure believe it is reasonable to assume a normally distributed β, in which case β β can be expressed as follows: In Equation (1), µ β and σ β are the mean value and the standard deviation of the reliability index.
We would like to point out that a false estimate of the statistical distribution of β β may also result in a significant error in the estimation of its value.
The robustness measure is calculated for every design in the design space, followed by a constrained, multi-objective optimisation using the NSGA-II genetic algorithm. The objectives in the optimisation procedure are the maximisation of design robustness, and the minimisation of the foundation construction price; foundation width and depth are optimised. Within the procedure, the feasible region is limited by applying constraints defined as follows: β ULS ≤ β T ULS and β SLS ≤ β T SLS , where β ULS and β SLS are reliability indexes, and β T ULS and β T SLS target reliability indexses for ULS and SLS, respectively. The optimisation result is a non-dominated set of designs which is called a Pareto front, or Pareto-optimal set. The choice of the final design is a trade-off between robustness and design costs.

A Modification of the RGD Procedure
As suggested in other studies [5,7,8], "noise factors" are treated as random variables in this paper; the number of input parameters is therefore reduced and simplified. In the case of insufficient data for a quality statistical analysis, the typical statistical characteristics of geotechnical parameters can be found in literature. The input data is further simplified by applying a sensitivity analysis, based on which the random variables, whose contribution to the considered system response variance is negligible, are "frozen". As a robustness measure, we use a generalised reliability index β, which is traditionally used as a reliability measure in classic reliability analyses, and by definition fits into the robustness concept. Introducing Eurocode 7 into the procedure additionally ensures the robustness of the designs. Unlike the original procedure, the results are directly applicable in everyday engineering practice. The key terms used in the procedure are defined as follows: Constraints are the criteria that divide the design space into feasible and non-feasible regions.

Harmonising the RGD Method Results with Eurocode 7
In the existing RGD method, a feasible region encompasses all designs that meet the pre-defined target failure probability requirement (P T f ), i.e., the associated value of the target reliability index (β T ). The connection between P f and β is mathematically defined, and these can therefore be used interchangeably. In this paper we use β because it is more user-friendly than P f . The associated range of P f for 0.5 ≤ β ≤ 10 would be 0.3 ≤ P f ≤ −7.6 × 10 −24 .
Eurocode 0 [12], in Annex C, gives a recommendation on limit values for β T , in accordance with the considered limit state, reference period and reliability class. Thus, for a 50-year reference period and reliability class RC2, the following recommendations for β T were given: β T ULS = 3.8, β T SLS = 1.5. The same values are used in the foundation optimisation using the original RGD procedure.
Meeting the criteria for the aimed reliability index value does not necessarily ensure meeting the criteria for ultimate and serviceability limit states defined in EC7. Forrest and Orr [13] study the reliability of shallow foundation designed according to Eurocode 7, in the case of ultimate limit state (ULS). In their study, they conduct a series of analyses, which consider foundation's reliability indexes according to EC7, and according to the allowable stress design format, with overall Factors of Safety FOS = 2, and FOS = 3. Different design situations are defined by the variance of the resultant force position (centric and eccentric load), soil type (granular and fine grain soils), relative load sizes (large and small loads), and the correlation tan φ − c (with and without correlation). They conduct multiparameter analyses for all the listed design situations, with the variance of geotechnical parameters: vertical scale of fluctuation (δ V ), tangent of the angle of internal friction (V tanφ ), characteristic values of φ (φ k:mean and φ k:low ). The presented results indicate that there are designs that meet the target reliability index, but do not meet the ultimate limit state criterion defined in EC7.
An example of such designs is shown in Table 2, which lists selected final designs attained through the application of the RGD procedure from the illustrative example given in [1]. The designs marked grey fail to meet the ULS criterion defined by EC7; therefore, additional assessments of limit states according to EC7 are necessary for the chosen design, along with the RGD procedure. By introducing additional constraints into the optimisation process, it is possible to exclude all designs that do not meet the ULS and SLS criteria defined in EC7 from the feasible region. Figure 2 schematically illustrates the design space and the existing and suggested feasible regions in the case of ULS. According to the existing RGD procedure, the feasible region is composed of regions 1 and 3 from Figure 2, and is defined as follows: Region 3 includes designs that meet the criterion concerning the aimed reliability index value, while at the same time failing to meet the criteria defined by EC7. By introducing additional constraints, region 3 is excluded from the feasible region, whose definition is then altered as follows:

Simplification of the RGD Procedure by Substituting the Robustness Measure
The most demanding part of the RGD procedure is the outer loop from Figure 1, which is composed of M repetitions of FORM with integrated PEM, where M is the total number of the designs in the design domain. Every design requires seven FORM and one PEM analysis, done separately for ULS and SLS. The number seven refers to the number of estimating points in PEM [14]. The aim of the mentioned procedure is the calculation of the mean value and standard deviation of β, which are used to determine the design robustness measure β β . Such a procedure is extremely computationally expensive, and includes differentiations of limit states' functions, which are mathematically complex in the case of shallow foundations. Furthermore, another shortcoming of applying β β is the fact that we need to estimate the probability density function (PDF) of β, which is unknown. Making a wrong assumption considering the PDF of β can lead to a significant mistake in the estimation of β β . This paper explores the possibility of substituting the robustness measure β β with a generalised reliability index β [15], which is defined as follows: In Equation (4), Φ −1 is the inverse Gaussian distribution and P f is the probability of failure.
Defining β is a part of the RGD procedure within the constraint verification part; consequently, by applying it as the robustness measure, the whole procedure is significantly simplified. We consider β to be a good indicator of the design's robustness, since its value is directly related to the standard deviation of the considered system response. Designs with a smaller standard deviation of the system response are, by definition, more robust. Figure 3 illustrates the relation of β to different values of the standard deviation of system response in shallow foundations in cohesionless soil. In the given example, a factor of safety (FS) was chosen for system response, with the assumption that it is lognormally distributed. The statistical distribution of the factor of safety in shallow foundations was studied by Dodigović et al. [16]; based on comprehensive statistical analyses, they concluded that such an assumption was justified. Figure 3 clearly shows that designs with higher β favour lesser values of the standard deviation of system response, making them less sensitive to the variations of input parameters-i.e., more robust.
In order to compare the results yielded from the application of the original robustness measure, we analysed the relation β β − β for different foundation widths and coefficients of variations of the soil's internal friction angle (COV φ ). The results of the analysis are shown in Figure 4. Due to a linear relation with a high correlation coefficient between β β and β, we conclude that the design within the Pareto front-but the shape of the Pareto front as well-will be very similar for examples where the robustness measure is β or β β .
It is possible to calculate the reliability index by using the PEM method, which is significantly simpler and computationally less demanding than the other reliability methods (FORM, SORM, FOSN, Monte Carlo). Zhao and Ono [14] suggest PEM with seven estimating points, used in the outer loop of the RGD procedure, whose accuracy is significantly greater than the earlier instances of the method using 2, 3 and 5 estimating points. The main shortcoming of PEM with seven estimating points is the fact that the random variables are transformed into standard normal space using the Rosenblatt transformation. In the example of correlated random variables, the transformation of random variables requires having total probabilistic information data, which is almost never the case in engineering. By substituting the Rosenblatt transformation with the Nataf transformation in PEM, we enable the consideration of correlated variables with the knowledge of the correlation matrix, their marginal distributions, mean values and standard deviations. The data on correlations between geotechnical parameters, their statistical distribution and standard deviation are available in literature, and the mean values can be determined from the geotechnical investigation results. The modification of PEM by introducing the Nataf transformation was suggested by Yu et al. [16], and they named the modified method IPEM. The Nataf transformation and its inverse are expressed in the following equations: where F X (.) is the marginal cumulative probability density function (CDF) of X, Φ(.) is the standard normal CDF of X, L 0 is the lower triangle matrix yielded from the Cholesky decomposition of the correlation matrix R 0 = ρ 0,ij . The procedure of determining the correlation matrix R 0 is composed of a series of complex function integrations [17], but its approximation is possible, based on a set of semiempirical equations suggested by Kiureghian and Liu [18]. R 0 is approximated based on the known correlations of ρ ij and the ratio of F in the following way: The value of the ratio of F is given for different sets of marginal distributions, divided into two groups. According to the distribution pair, tables are given which provide the values/equations of F.
The IPEM procedure is conducted using simpler mathematical operations. The inverse Nataf transformation, which is a part of IPEM, is required only at the estimating points, and can be done simply, e.g., by using the built-in function of Microsoft Excel, or the Python program language used with the "SciPy" package. Since the accuracy of the calculations of β is highly dependent on the method used [19], the accuracy of IPEM for ULS and SLS of shallow foundation was examined. The Monte Carlo method was chosen as the control method, and the results are given in Tables 3 and 4.  The average error in the estimation of the reliability index using the IPEM method for ULS and SLS of shallow foundations, for different variation coefficients of the internal friction angle, is approximately 1%-which is negligible.

The Optimisation of the Number of Random Variables Using a Sensitivity Analysis
Obtaining non-dominated designs by applying the genetic algorithm in the RGD procedure can be optimised by reducing the number of random variables included in the reliability analysis. In the reliability calculation, failure is a probabilistic event, and its probability is given by: where X is a random vector and g(.) is a limit state function.
In the case of a shallow foundation loaded by permanent and variable vertical action in cohesionless soil, the function of the limit state for ULS can be expressed as follows: (9) where N q i N γ are load capacity coefficients, s q and s γ are the coefficients of foundation shape, γ is the effective unit weight of soil, B is the effective foundation width, q is the effective overburden pressure at the foundation base level, V g is the permanent vertical load, V q is the variable vertical load. The function of the serviceability limit state can be defined by applying a formula developed by Akbas and Kulhawy [20], which is derived from load-settlement behaviour of a shallow foundation under axial compression loading where s t is the tolerated settlement, B is the foundation width, and coefficients a and b are the parameters of the hyperbolic model that fit the normalised-settlement curve. Random vectors from Equation (8) for the ultimate limit state can be expressed as X ULS = φ , γ , V g , V q T ; for the serviceability limit state, the following applies: The Sobol sensitivity analysis is conducted, with the aim of determining the contribution of the variance of each individual random variable to the total variance of limit state functions for both ULS and SLS. The Sobol sensitivity analysis is a variance-based method which determines the influence of uncertainties in the model input factors on the uncertainty in the output of the model [21]. By determining the first order Sobol indices, we determine the effect of the variation of individual random variables on the variance of system response. The Sobol indices were calculated in the Python program, using the open source package "OpenTURNS" [22]. The variables whose contribution to the system response variance is negligible can be "frozen", which will make the RGD procedure less computationally demanding. Tables 5 and 6 show the results of sensitivity analyses for different COV φ values, for both ULS and SLS. The presented results demonstrate that the soil's internal friction angle is the dominant variable and contributes most to the total variance of system responses. The value of the first order Sobol index is similar for ULS and SLS, and for different COV φ values it is within in the ≈ 0.84 − 0.99 range, depending on the value of the internal friction angle. Due to the high values of first order Sobol indices of the soil's internal friction angle, reliability analyses of the ultimate and serviceability limit states were conducted for two groups of random vectors, which differ in the number of random variables. The first group of random vectors is composed solely of the soil's internal friction angle, and can be expressed as follows: X ULS = (φ ) T and X SLS = (φ ) T . In the second group the random vectors include all the random variables defined by limit state functions for ULS and SLS: X ULS = φ , γ , V g , V q T and X SLS = φ , γ , a, b, V g , V q T . The aim of the analysis is to determine the influence of reducing the number of random variables on the reliability analyses results, in order to optimise the RGD procedure. Reliability analyses were conducted using the IPEM method, and the results are presented in Tables 7 and 8.  The average difference in reliability indexes β 1 and β 2 of the results presented in Tables 7 and 8 for ULS is 1.3%, and 1.42% for SLS, with maximum deviations of 2.36% and 3.25%. For the value of COV φ = 0.1, as suggested in literature, deviations are smaller-0.56% for ULS and 0.41% for SLS. Such small deviations are in line with the results of sensitivity analyses presented in Tables 3 and 4. We therefore conclude that the errors in calculating the reliability index, stemming from freezing random variables γ , V g and Vq in the ULS reliability analysis and γ , V g , Vq, a and b in the SLS reliability analysis, are negligible. Consequently, we propose conducting reliability analyses for ULS and SLS in which the internal friction angle is the only random variable. An analysis of the influence of ODF values on first order Sobol indices was also conducted. The results of this analysis are not presented, but it was confirmed that ODF has no significant influence on the value of the first order Sobol indexes for both ULS and SLS.
We would like to point out that the limit state functions for shallow foundations can also be determined by using other calculation models. This paper used the analytical expression for estimating the bearing capacity of a shallow foundation which is recommended in Eurocode 7, Annex D [23] (Equation (9)). This expression is based on a theory extrapolated from the results of laboratory tests on the behaviour of small footings in dense sand and under 1 g condition [24]. Altaee and Fellenius [25] are researching soil behaviour in such conditions and conclude that it has little relevance to the behaviour of a full-scale prototype. The main reasons for this are inadequate stress levels in soil and non-linear stress-strain soil behaviour. With small-scale tests, the depth of model influence is relatively small, so soil behaves as if unconsolidated-after an initial volume assessment, soil dilates and then contracts. This results in a stress-strain curve which suggests ultimate resistance; this has not been proven in full-scale tests to this day. Due to all of the above, Fellenius [26] recommends bearing capacity analysis to be used only as a simple estimate for comparing the design with previous designs.
Despite the mentioned shortcomings and simplifications, Eurocode 7 recommends using the expression from Equation (9) for determining the bearing capacity of a shallow foundation; it is consequently used in everyday engineering practice. Since the aim of this paper is to present the procedure of the modified RGD method whose results have been harmonised with EC7, the bearing capacity of a shallow foundation is determined according to the recommendations from this technical standard. Limit state functions can also be defined using other calculation models following the same principle, which will result in different mathematical expressions of limit state functions.

Optimisation Using the Non-Dominated Sorting Genetic Algorithm II (NSGA-II)
Genetic algorithms (GA) are the heuristic optimisation and search techniques that mimic nature's evolutionary principles. The concept of the genetic algorithm was developed by John Holland, and it is used in solving problems from various problem domains, including sciences, commerce, and engineering [27]. The non-dominated sorting genetic algorithm II (NSGA-II) was first suggested by Deb et al. [28] with the aim of improving the existing multi-objective evolutionary algorithms that use non-dominated sorting. A flowchart of the NSGA-II algorithm is shown in Figure 5. In the original RGD procedure, objectives are calculated for every design in the design space, followed by non-dominated sorting to determine the Pareto front ( Figure 1). Such an approach results in a high number of reliability analyses, which depends on the number of designs in the design space, as well as on the number of noise factors. In the example from [1], the number of designs is 450, and the number of noise factors is 4 and 5 for ULS and SLS, respectively. In that case, the RGD procedure includes 7 × 4 × 450 = 12,600 FORM and 450 PEM analyses, and 7 × 5 × 450 = 15,750 FORM and 450 PEM analyses for SLS (the number 7 refers to the 7 estimating points according to the PEM procedure).
We suggest calculating objectives only for the designs chosen in the procedures included in the NSGA-II algorithm. In that case, the number of computations for every objective depends on the given termination criterion, which can be, e.g., the number of function evaluations, the number of iterations, or an advanced criterion based on the defined performance metric.
Due to a number of reasons, applying any of the algorithmic methods will rarely enable determining the true Pareto front, with only its approximation being possible. The approximation quality can be measured using various indicators, among which the hypervolume indicator (HI) [29] is highly important. An increase in the value of HI is an indicator of convergence towards the true Pareto front. When its value ceases to change significantly, in relation to the number of evaluations of objective functions, the algorithm is considered to have converged. An example of the hypervolume indicator for a twoobjective case is shown in Figure 6, and includes a hatched region bound by the reference point r. We analysed the performance of the NSGA-II algorithm in the case of the RGD procedure, modified as described in Sections 2.2.1-2.2.3. An optimisation was conducted, with three objectives defined as follows: minimise foundation cost, maximise ULS and SLS robustness with the decision variables being foundation width (B) and depth (D).
Since the choice of genetic algorithm parameters directly influences the quality of solutions and convergence [30], we conducted a parametric analysis of crossover and mutation of operator parameters, with the aim of determining their optimum values. The analysis was conducted for different COV φ values. The work included the usage of a simulated binary crossover (SBX) operator, which was shown to be efficient for real variables [31]. The parameters of the SBX operator are crossover probability (p c ) and distribution index (η c ). Crossover probability is the number of realised crossovers in one generation. If its value is 0%, then the entire new generation equals the preceding one; if the crossover rate is 100%, the entire generation is substituted with new offsprings, yielded from the crossover of units in the previous generation. The distance of the offsprings from the parent solution will depend on the value of η c : if η c is large, the resulting offsprings will be near the parent solution, with the opposite being the case for smaller values. The mutation operator ensures the maintaining of genetic diversity in the genetic algorithm population. Deb and Deb [32] analyse the usage of polynomial and Gaussian mutation operators for real-parameter genetic algorithms. They conclude that both operators are equally efficient, and this work used the polynomial mutation operator with a defined mutation probability (p m ) and distribution index (η m ). The value of the p m parameter was chosen according to the proven efficient expression p m = 1/n, where n is the number of decision variables [30]. The parametric analysis yielded optimal parameters of the NSGA-II algorithm: p c = 0.9, η c = 30, p m = 0.5 and η m = 20. The relation between normalised hypervolume and the number of evaluations of objective functions, for the mentioned parameters, is shown in Figure 7 and Table 9.   Figure 7 and Table 9 both show the quick convergence of the NSGA-II algorithm. Compared to the original RGD procedure, which requires 450 evaluations for the same input parameters, the suggested approach reduces that number by 40-50%.
By adopting the modification suggested in the Sections 2.2.1-2.2.4, the RGD procedure is simplified, and its flow diagram is shown in Figure 8. During the first step, it is necessary to determine the statistic characteristics (mean value, standard deviation, statistical distribution) for all random variables and their mutual correlations. The number of random variables depends on the chosen geotechnical model of ultimate and serviceability limit state. The next step is conducting a sensitivity analysis with the aim of freezing the variables whose contribution to the variance of system response is negligible. When all the random variables are defined, we need to choose a decision variable, which makes the optimisation problem completely defined. The final step is composed of conducting a multi-objective optimisation using the NSGA-II algorithm. Unlike the original procedure, the inner loop ( Figure 1) is not conducted and determining the robustness measures from the outer loop is done within the NSGA-II algorithm, which reduces the required number of objective functions evaluations. The result of the optimisation is a set of non-dominated solutions, out of which the final design is chosen.

Illustrative Example 1-Design of Shallow Foundation
The problem of optimising a shallow foundation on non-coherent soil was considered, similar to the one presented by Juang et al. [1], from which we took the characteristics of the used variables. As shown in Figure 9, the foundation is subjected to vertical permanent (V G ) and variable (V Q ) loads. Aside from the soil's internal friction angle, all other random variables are "frozen", since their contribution to the total system response variance is negligible, as shown in Section 2.2.3. A three-objective, constrained optimisation was conducted (ULS, SLS, foundation cost) with two decision variables: foundation width (B) and depth (D). A square foundation layout is assumed (B = L).  (5) and (6), and maximum tolerated settlement is set at 25 mm. The estimation of foundation cost is conducted using the expression defined by Wang and Kulhawy [33]: where Q e , Q f , Q c , Q b signify the quantities of excavation, formwork, concrete and compacted backfill, respectively, and c e , c f , c c , c b the associated unit prices. More on the cost calculation for a shallow foundation can be found in [33]. The foundation soil is dry sand with the following statistical characteristics of the internal friction angle: µ φ = 34.4 • and σ φ = 0.59. The characteristics of the deterministic variables are shown in Table 10. The problem of shallow foundation optimisation in the case in question can be written out as follows: The solutions for the optimisation are shown in Figures 10 and 11. A set of nondominated designs (Pareto front) in objective space is shown in Figure 10 on the left, and in decision space in Figure 10 on the right. Aside from the Pareto front (magenta), also shown are the dominated feasible designs (green) and infeasible designs (black). The design space is composed of a total of 451 designs, 311 of which are feasible and 140 infeasible. The optimisation resulted in 71 non-dominated designs which make up the Pareto front. The minimal foundation dimensions from the Pareto front are B × D = 2.1 × 1.8 m, with the maximum being 5.0 × 2.0 m, at the costs of USD 1160 and USD 5687, respectively. Figure 10, on the right, shows the grouping of non-dominated designs in the vicinity of the maximum prescribed depth of D = 2.0 m, which is in accordance with the original method's solutions presented by Juang et al. in [1]. Figure 11 shows the projections of the Pareto front onto axes β ULS and β SLS . The minimum values of the reliability index for ULS and SLS are 9.4 and 3.2, respectively, which is significantly higher than the prescribed 3.8 and 1.5, respectively. This fact is a consequence of a relatively high value of the internal friction angle's standard deviation, along with the fact that the prescribed condition ODF ULS ≥ 1 is more critical than the criteria associated with the reliability indexes.  We analysed the influence of individual constraints on the position of non-dominated designs in the decision space, and the results are shown in Figure 12. The four broken lines represent the boundary between feasible and infeasible designs for different constraints. Left of the line are infeasible designs, with feasible designs being right of the line and on the line itself. Figure 12 clearly shows that the constraint ODF ULS ≥ 1 defines the position of the Pareto front in the design space. Comparing to the other constraints, to meet ODF ULS ≥ 1 constraint, a wider foundation width at the same depth is needed. Such a relation among constraints is the consequence of a relatively low value of the standard deviation of the soil's internal friction angle (σ φ ). Reliability indexes are more sensitive to a change in standard deviations of random variables than ODF. To illustrate this point further, Figure 12-on the right-shows lines that represent the boundaries of the distribution of design space into feasible/infeasible designs, in the case of a larger standard deviation of the internal friction angle, σ φ = 3.44. In that case, the number of feasible designs is greatly decreased (from 311 to 236), and the main reason behind the decrease is a significant reduction of the reliability index, which led to an additional filtering of the designs. In this case, the position of the Pareto front is defined by the constrains β ULS ≥ 3.8 and β SLS ≥ 1.5, while the constraints relating to ODF ULS and ODF SLS are not as critical. A comparison was made between the optimisation results yielded from the original and the modified RGD procedure. The feasible region of the original RGD procedure is composed of 331 designs, while the modified RGD procedure includes 310 designs. The results of the comparison are presented in Figure 13. The modified RGD procedure resulted in more conservative solutions, as seen from the boundary between the feasible and infeasible region, which is shifted slightly towards the right (Figure 13). A reduced number of designs in the feasible region in the modified RGD procedure, as opposed to the original procedure, is the consequence of additionally introduced constrains, according to the criteria of limit states prescribed in Eurocode 7.

Illustrative Example 2-Design of Axially Loaded Pile
The modified RGD method is applied for the design of an axially loaded pile in a stiff clay, whose geometry is shown in Figure 15. The soil parameters and external action are taken from the example given in [34]. It is necessary to optimise pile length (D) and diameter (B), with the aim of maximising robustness while minimising pile cost. The "α method" [35] is used for determining the bearing capacity of the pile, and the ULS limit state function is defined as follows: (12) where R b = (N c c u + σ v0 ) D 2 π/4 and R s = α u c u DπL are the pile base and shaft resistances [35], W p is the weight of the pile, N c = 9, c u is undrained shear strength, σ v0 is total vertical stress at the depth of the pile base, α u is the coefficient used to relate c u to the adhesive stress along the pile shaft.
The SLS limit state function is defined by using the equation which connects the bearing capacities of ULS and SLS, attained through the statistical analyses of the results of pile load tests [36]: where Q ULS and Q SLS are the bearing capacities for ULS and SLS, s a is allowable settlement (20 mm), a and b are hyperbolic curve-fitting parameters for the normalised load-settlement curve. Using Equation (13), the SLS limit state function can be expressed as follows: The cost of the pile (C) is related to the concrete volume. In the presentation of applying of the RGD method to the case of an axially loaded pile, Equation (15) was used, in which the cost is estimated by the following expression [34]: The pile is subjected to axial permanent (V G ) and variable (V Q ) loads. Undrained shear strength c u is determined from a triaxial test, and its mean value increases linearly with depth, and according to the following expression: c u = 13.0z, where z is depth. The characteristics of random variables are shown in Table 11. In order to determine the influence of individual random variables on the variance of ULS and SLS limit state functions, a sensitivity analysis was conducted using Soboli indices. The analysis results are given in Table 12. Based on the conducted sensitivity analysis, the variables with small contributions to the system response variances are frozen. The reliability analyses are conducted using random vectors with a reduced number of random variables, defined as follows: X ULS = (c u , α) T and X SLS = (c u , a, b) T . The next step in the implementation of the RGD method is the identification of decision variables and their lower and upper boundaries. In this case, the following decision variables were chosen: -length of the pile: 9 ≤ L ≤ 18 m; -diameter of the pile: 0.5 ≤ D ≤ 0.8 m.
After defining limit state functions, random and deterministic variables, and decision variables, the problem of optimisation according to the modified RGD method can be expressed as follows:  Figure 16 illustrates the convergence of the NSGA-II algorithm by measuring hypervolume indicators. After conducting 337 evaluations, the algorithm converged. The yielded Pareto front is composed of a total of 107 non-dominated solutions. Since the feasible region is defined according to Equation (3), all non-dominated solutions meet the ULS and SLS criteria prescribed in EC7. The optimisation was also conducted by applying the original RGD method. Since the feasible and infeasible regions match the modified RGD method's results completely, they are not shown separately. Figure 17 shows optimisation solutions in objective and decision space. Figure 17, on the right, illustrates a tendency of grouping non-dominated solutions with bigger pile lengths for all considered diameters.
Interruptions in the Pareto front, visible in Figure 17 on the left and in Figure 18, are the consequence of the discretised space of the decision variables, and the significant influence the diameter has on pile cost. Each of the curves show designs with the same diameter. The lowest curve matches pile diameter D = 0.4 m, while the highest matches pile diameter D = 0.8 m.

Discussion
This paper presented the modified RGD method, which is used in geotechnical engineering to determine optimal, robust designs which meet the safety and cost requirements. The method is computationally expensive and its application in everyday engineering practice is complex. This fact is mainly the consequence of a large number of reliability analyses within the inner and outer loop of the RGD method, which needs to be performed to calculate the feasibility robustness index β β . The procedure consists of the FORM method with integrated PEM, as shown in Figure 1. To calculate β β , we need to know the standard deviations of "noise factors". In the RGD method, aside from geotechnical parameters, their standard deviations are also treated as "noise factors", so a knowledge of their statistical characteristic is required. In everyday practice we often do not have sufficient data for a quality statistical analysis, so their values can only be roughly estimated. Tan et al. [4] simplify the calculation of β β by applying the sensitivity of reliability index. Doing so reduces the number of reliability analyses for M (the number of designs in the design space), but still uses mathematically complex operations, and input parameters remain unaltered. Different authors introduce changes into the design robustness measure [5][6][7][8][9][10], which also greatly reduces the number of reliability analyses for M, and simplifies the input parameters; however, complex mathematical procedures are still applied, unlike the ones used in everyday engineering practice.
Unlike previous modifications, we-for the first time-suggest a modification of the RGD method that harmonises the optimisation results with the criteria prescribed in Eurocode 7. By doing so, we eliminated the need for additional assessments of the designs, i.e., the RGD method can be used as an alternative to the traditional geotechnical analysis. As the design robustness measure, we suggest using a generalised reliability index (β), which is traditionally used in engineering practice as the measure for structural reliability. We demonstrated that β by definition fits into the reliability concept, and its relation to β β is linear, as shown in Figure 4. Due to the linear relation β − β β , the feasible and infeasible regions, as well as the Pareto fronts yielded from the original and from the modified methods, are similar-as illustrated in Figures 13 and 14. We suggest calculating the generalised reliability index using the modified PEM method (IPEM), which is simpler to use, i.e., it does not include mathematically complex operations, while maintaining a high degree of accuracy-as presented in Tables 2 and 3. We introduce a sensitivity analysis into the RGD method, with the aim of "freezing" the random variables whose contribution to system response is negligible. For this purpose, we suggest calculating first order Soboli indices, which can also be estimated using the IPEM method. To further reduce the number of conducted reliability analyses, we suggest conducting them only for the designs chosen via the genetic algorithm NSGA-II. This reduced the number of analyses in the example of a shallow foundation on coherent soil from M to 0.5 − 0.62M, as presented in Figure 7 and Table 9. The given illustrative example 1 shows a good match of results yielded from the original and modified methods, illustrated in Figures 13 and 14. According to expectations, and due to the introduction of additional constraints, the feasible region of the modified method was moved right in relation to the original method.
In the case of the axially loaded pile shown in the illustrative example 2, the feasible and infeasible regions calculated using the original and modified methods were the same. The reason behind this is a relatively high variance of the random variable (c u ), which is why the position of the Pareto front in decision space is defined exclusively by the constraint related to β ULS .
A satisfactory result match was achieved, with significantly simpler input parameters, a smaller number of conducted reliability analyses, the usage of simpler mathematical operations, and greatly improved computational efficiency. The suggested modifications simplify the RGD method, drawing it closer to engineers and its usage in everyday engineering practice.
The paper gives two simple examples with the aim of presenting the procedure of conducting the modified RGD method. In case soil conditions are more complex, the procedure remains the same, provided there are analytical expressions of limit state functions. When these expressions are unknown, i.e., when it is impossible to express them mathematically, their estimation is made possible by using more sophisticated methods, such as the finite element method. This can be conducted using the Plaxis software, which enables conducting analyses by running scripts written in the Python program language. For further research, we recommend exploring the possibilities of integrating the finite element method into RGD, which could be used for more complex geotechnical tasks.
To further improve the RGD method, we suggest researching the influence of substituting a discretised design space with a continuous one. Since the designs included in the Pareto front are grouped close to the maximum prescribed foundation depth (D max ), there is a possibility that the proposed substitution results in a Pareto front consisting only of designs with the foundation depth of D max . This could potentially reduce the number of decision variables. Furthermore, we noticed that the variation of the crossover and mutation parameters significantly influences the convergence of the NSGA-II algorithm. We suggest further research into the application of the improved NSGA-II algorithm with self-adaptive crossover and mutation parameters, with the aim to further reduce the number of the required reliability analyses.