Revisiting Two Simulation-Based Reliability Approaches for Coastal and Structural Engineering Applications

Featured Application: Normality polynomials can be used to compute reliabilities for coastal and structural engineering applications, including the assessment of uncertainty in the estimated reliability index. Additionally, multi-linear regression can be applied to the simulated results to determine design points and sensitivity factors. These applications can be potentially extended to di ﬀ erent engineering (or other) ﬁelds and to system reliability (e.g., for reinforced concrete frame buildings). Abstract: The normality polynomial and multi-linear regression approaches are revisited for estimating the reliability index, its precision, and other reliability-related values for coastal and structural engineering applications. In previous studies, neither the error in the reliability estimation is mathematically deﬁned nor the adequacy of varying the tolerance is investigated. This is addressed in the present study. First, sets of given numbers of Monte Carlo simulations are obtained for three limit state functions and probabilities of failure are computed. Then, the normality polynomial approach is applied to each set and mean errors in estimating the reliability index are obtained, together with its associated uncertainty; this is deﬁned mathematically. The data is also used to derive design points and sensitivity factors by multi-linear regression analysis for given tolerances. Results indicate that power laws deﬁne the mean error of the reliability index and its standard deviation as a function of the number of simulations for the normality polynomial approach. Results also indicate that the multi-linear regression approach accurately predicts reliability-related values if enough simulations are performed for a given tolerance. It is concluded that the revisited approaches are a valuable option to compute reliability-associated values with reduced simulations, by accepting a quantitative precision level.


Introduction
Simulations are often used to estimate the probability of failure of structural elements and systems because they are a very versatile option which is not restricted by complex and implicit limit state functions (LSF), the use of sophisticated methods (e.g., finite element method), and/or highly non-linear structural behavior. However, millions of crude Monte Carlo simulations (MCS) could be required to adequately estimate structural probabilities of failure, which may not be feasible; also, the results could differ from a set of simulations to another. To cope with this issue, modified versions of the crude simulation approach, surrogate modeling, subset simulation, and other techniques have emerged and been used in the last decades not only for structural engineering but also in other fields. To mention only a few studies which employ some of these techniques, optimization using surrogate modeling, reliability analysis of deteriorating structural systems and the reliability assessment of a structure affected by chloride attack are reported in the studies by [1][2][3], respectively. Importance sampling has also been used to estimate the system reliability of deteriorating pipelines [4]. Other kind of approaches are also reported in the literature to compute reliabilities [5]. The different reliability methods can be applied not only to different fields in structural and geotechnical engineering, as for the case of sudden column removal in reinforced concrete buildings and rockfall protection structures [6][7][8], but also to many other research and engineering fields, for instance to artic oil and gas facilities [9] and to coastal engineering applications [10]. This last case is used in the present study to show the applicability of two revisited methods to obtain the reliability of coastal structures.
Other simulation-based reliability methods have not been given so much attention. In this study a couple of these alternatives are revisited to inspect their feasibility and adequacy to estimate reliability indices. One of them employs polynomial transformations of a nonnormal variable to a normal one by fitting simulated data with fractile constraints and can be referred to as normality polynomial approach [11]. The second approach was developed to derive the design point and sensitivity factors (in the FORM, first order reliability method, perspective) from simulated data [12] and it is referred to as the multi-linear regression approach in this study. A similar approach to the normality polynomial method was previously developed by Hong and Lind [13] and named normal polynomial approach (note that the names are slightly different); unlike the normality polynomial approach, the normal polynomial approach has been paid much more attention (judging by number of citations), even very recently (e.g., [14]). Both methods are based on the fact [15,16] that a fractile of a random variable can be expressed as polynomial of a fractile of a standard normal variable (normal polynomial), and that a fractile of a standard normal variable can be expressed as a polynomial of a fractile of a random variable (normality polynomial). Although the normal polynomial approach [13] is not considered in the present study (we prefer to focus in the less explored alternative), the findings here could be extended to do so.
To inspect the adequacy of the normality polynomial approach and the multi-linear regression approach (the methods revisited in the present study), three LSFs are considered. One is based on a very simple classical case; other one is based on a structural application from a previous study and the last one is related to the reliability of a coastal structure. Extensive simulations are performed to estimate the error level by using the normality polynomial approach and its associated uncertainty; neither of these was thoroughly carried out in previous studies, nor the application to coastal engineering. The design point and sensitivity factors (also the reliability index) obtained from simulated data are compared with those obtained from FORM; they are derived from a multi-linear regression of the simulated data. It is worth to mention that these methods, and others developed in the 1990s, used to state that a large number of simulations were not feasible; nevertheless, the computer power has increased substantially in the last decades, and the limitations of those days may not be as restrictive as before and thus the applicability could be currently extended. Furthermore, the commercial software available nowadays for engineering applications normally includes amenable built-in functions for linear and multi-linear regression analysis which simplifies the programming.
The main objective of this study is to define the error statistics in the reliability index by using normality polynomials and to reassess the feasibility and adequacy of this method and the multi-linear regression approach for estimating the reliability of structural and coastal engineering systems including the determination of sensitivity factors.
Appl. Sci. 2020, 10, 8176 3 of 21 This study is of significance to the coastal and structural engineering fields because the number of simulations required to compute reliabilities can be reduced by accepting defined error levels when using normality polynomials, which was not established in previous studies. This is possible because the error in the estimation of the reliability index is mathematically defined as a function of the number of simulations for the cases investigated. Additionally, other contribution is the use of multi-linear regression applied to simulated results as a mean to determine sensitivity factors, design points and the probability of failure; not only a slightly modified (improved) version of the multi-linear regression approach but also the number of simulations and tolerances required to achieve adequate results are provided for guidance.

Normality Polynomial
In this section the normality polynomial proposed in [11] is described. The mathematical form is given by where z p denotes p-fractile of a standard normal variable Z with probability density function (PDF), Φ(z), and cumulative distribution function (CDF), Φ(z); y p denotes p-fractile of a random variable Y with PDF, f Y (y), and CDF, F Y (y); a j , j = 1, 2, . . . , r, are the coefficients of a rth-order polynomial determined by fractile fitting. The fractile fitting is based on considering the following fractile constraints from a set of independent random observations of Y (i.e., y 1 , y 2 , . . . y i , . . . y n ) arranged in ascending order These fractile constraints can be mapped into a normal space by using where Φ −1 (•) denotes the inverse of the standardized normal distribution function. The rth-order polynomial (Equation (1)) with r + 1 < n is used to model the distribution of the transformed random variable Y. By considering the constraints in Equation (2), the coefficients a j in Equation (1) can be determined using the least square method by minimizing the error ε fit given by where m is the number of constraints. The probability P (Y ≤ y 0 ) is given by the CDF, i.e., F(y 0 ), and can be computed with where z p is obtained by substituting y 0 instead of y p in Equation (1). If Y is the resulting random variable of the LSF, the probability of failure, p f , is and the reliability index is [17] From Equation (7), it can be inferred that the reliability index can be readily obtained once the coefficients a j are determined (i.e., β = −a 0 ). This is so, because a 0 has a similar meaning to the so-called Appl. Sci. 2020, 10, 8176 4 of 21 generalized reliability index [11]. It is noted that the roots of the polynomial are not required [11] (unlike in the method in [13]). In the applications shown later, a third-order normality polynomial is used [11].
It is pointed out that although the mathematical background of the normality polynomials is not thoroughly described here, it is based on sound grounds [11], like the advance theory of statistics [15,16] and the fact that fractile constraints hold by a combinatorial argument [18,19].
Before going to the applications, the second revisited approach in this study is described in the following section.

Design Point and Sensitivity Factors from Simulations
The obtaining of the design point and sensitivity factors (from the FORM standpoint) based on a previous study [12] is described herein. The basic idea is that simulated data close to the limit state surface (within a prescribed tolerance, e) can be retrieved, as well as their associated sampled values for each or the random variables given by the LSF (input values combinations considered as a point in the hyperspace), and a multi-linear regression is performed to approximate the linearized limit state surface at the design point. Before the multi-linear regression is performed to fit the hyperplane, the considered points are mapped into a standard normal space. Such hyperplane is an approximation of the LSF, g, and is used to assess the design point and sensitivity factors. The mathematical formulation is given below.
A set of n independent random variables is denoted by X = (x 1 , x 2 , . . . , x n ). X j defines the j-th randomly generated value of X. For a given number of crude Monte Carlo simulations, n sim , X js which satisfy the criterion below (slightly changed from the original formulation in [12]) are selected where X k in the LSF is included to emphasize that it is a function of a set of random variables. Hong and Nessim [12] used e = 0.05 (i.e., 5%) in their study (instead of e l in Equation (8) and defined below). However, it was noticed that this value could be inadequate depending on the units and magnitude of the considered random variables. Therefore, the distance between zero and the smallest simulated value of g in absolute terms (l low ) is used to set the tolerance as the fraction given by l low multiplied by e (i.e., e l = e × l low is used instead of e in Equation (8)). The selected values of X k based on the described criterion are then mapped into a standard normal space [17] using where Z is the image of X in standard normal space, Φ is the normal standard CDF and F i , I = 1, 2, . . . , n, are the CDFs of the random variables x i . In the standard normal space, a linear function is fitted to the set of selected points using multi-linear regression. Such linear function is given by where b i and c are constants to be determined in the multi-linear regression analysis. The resulting linear regression equation is to be used in the same sense as the FORM [17] to estimate the design point and sensitivity factors. The latter are denoted by α i , i = 1, 2, . . . , n, and given by the gradient vector of Equation (10) as indicated below Note that the reliability index can also be estimated as the smallest distance between the linear surface and the origin as Appl. Sci. 2020, 10, 8176

of 21
As regards the design point (in the normal space) it is given by The design point in the original space can now be determined with the inverse transformation of z id , as The subscript i in Equations (11)- (14) is associated to the random variable x i . The inverse transformation in Equation (14) is dependent on the PDF of x id and, in the case of the Longuet-Higgins distribution [20] used for the coastal engineering example, the inverse transformations of z id requires of a numerical approach to be determined.
The formulation in this section, and the one in the previous section, are to be applied to three case studies in the following section to evaluate their adequacy for structural and coastal engineering, and to assess their deviation with respect to the exact reliability index.

A Classical Limit State Funtion
The two approaches described in the previous section are applied here to the simplest classical LSF where R and L can be considered, in a broad sense, as random variables for the capacity and demand of a system element. For the sake of simplicity and illustration purposes, units are skipped and both random variables are assumed independent and characterized with lognormal distributions, with mean values and standard deviations m R = 10, m L = 5.6, σ R = 1, and σ L = 0.75 for the capacity and demand, respectively. These values are arbitrary, except by the fact that they lead to a reliability index equal to practically 3.5, which is a common reference for code calibration and that can be computed with the following expression [21,22] where ν R and ν L are the coefficients of variation of R and L, respectively. This reliability index is shown in Figure 1a (dashed line), as a reference to inspect how close are the βs obtained by crude Monte Carlo simulations as a function of the number of simulations (shown in logarithmic scale in the horizontal axis from 2 × 10 1 to 2 × 10 7 ) to the exact value. The reliability index for Equation (15) using Monte Carlo simulations (dashed-dotted line in Figure 1a) is obtained by plugging into Equation (7) the following probability of failure p f = n f ail /n sim (17) which is simply the ratio of number of failures, n fail , to the total number of simulations. The latter (i.e., n sim ) is also the number of fractile constraints when the reliability index is computed with the normality polynomial approach, also depicted in Figure 1a (solid line). To quantitatively inspect these deviations, βRE is to be used as benchmark to assess the exactness of the used methods (i.e., normality polynomials and MCS) in terms of the relative percentual error given by where βbench denotes the reliability index considered as benchmark and βoth is the reliability index computed with any other method. A total of 1000 runs are performed for each nsim, and ε is computed for each of the runs. Then, the mean values of the error and their uncertainties (the unbiased standard deviation) are computed and plotted for the whole range of nsim in Figure 2 (solid lines and dashed lines for normality polynomial and MCS mean errors, respectively), where the mean errors ± one standard deviation are also depicted in grey lines. Figure 2 shows that is not always possible to compute the statistics for the MCS; this happens when not a single failure is reported in one or more of the 1000 runs for a given nsim (at least 5 × 10 4 are necessary). This is not a problem when using the fitted polynomials. Additionally, the MCS approach always tend to larger mean errors and standard deviations for a decreasing number of simulations. This makes the normality polynomial approach more adequate for estimating the reliability indices; however, for few simulations the errors are too large. Nevertheless, the designer could decide which precision level (quantitatively) is willing to accept using information like the one in Figure 2 as an aid (and reduce the number of required simulations as a function of such an accepted error). Additional runs are shown in Figure 1b-d, which indicate that the results are different and dependent on the generated random numbers of each run, but they stabilize if enough simulations are performed or enough fractile constraints are used. Other observations from Figure 1 include that β cannot always be computed with the Monte Carlo simulations (MCS) (not a single failure is obtained), while the opposite occurs when using normality polynomials, although significant deviations are observed for a limited number of simulations, that the fitted normality polynomials tends to deviate less from the exact reliability index for fewer n sim , and that such error in the precision may not be large for a relative small n sim , (e.g., 1 × 10 3 ).
To quantitatively inspect these deviations, β RE is to be used as benchmark to assess the exactness of the used methods (i.e., normality polynomials and MCS) in terms of the relative percentual error given by where β bench denotes the reliability index considered as benchmark and β oth is the reliability index computed with any other method. A total of 1000 runs are performed for each n sim , and ε is computed for each of the runs. Then, the mean values of the error and their uncertainties (the unbiased standard deviation) are computed and plotted for the whole range of n sim in Figure 2 (solid lines and dashed lines for normality polynomial and MCS mean errors, respectively), where the mean errors ± one standard deviation are also depicted in grey lines.  If desired, the curves in Figure 2 could be casted as a mathematical expression. For instance, the following power equations fit very well the mean (με) and standard deviation (σε) of the error in Figure 2 (by fitting from 2.5 × 10 2 simulations and over) If the power in Equations (19) and (20) is assumed as −0.5 in both expressions, the coefficient of variation of the error is νε ≈ 150/250 ≈ 0.6, i.e., it is constant and roughly independent of the number of simulations; the actual νε obtained from the 1000 runs does exhibit such a roughly constant behavior for this case, except that it is ≈0.7 (difference related to the actual different powers in Equations (19) and (20)). Power equations like Equations (19) and (20) can be linearized by taking logarithms on both sides. Therefore, if the involved variables are transformed into the logarithmic space, a linear fitting can be performed. In this study, we simply used a built-in function in the commercial software for the fitting.
Coefficients for the fitted normality polynomial in Figure 1a are shown in Table 1 (upper set of values) for selected values of nsim. The computing of these coefficients is based on minimizing the error in Equation (4) and was implemented in the coded program by using a built-in function of the programming language employed (MATLAB). As mentioned before, the coefficients a0 can be linked to the generalized reliability index. Coefficients for results in Figure 1b-d (or those associated to Figure 2) were computed but not shown for brevity.
As regards the reliability index, design point and sensitivity factors derived from the simulations (i.e., those obtained with Equations (11), (12), and (14)), they are compared with those obtained by applying the FORM to Equation (15). They are summarized in Table 2 for selected values of nsim.
Results listed in Table 2 indicate that when the points obtained by applying the criterion in Equation (8) were enough to successfully perform a multi-linear regression (also with a built-in function, as in the case of the normality polynomial fitting), β did not deviate from the exact value, but marginally, just like the reliability index obtained with FORM. However, a relatively large nsim was required, usually at least 1 × 10 5 simulations (by inspecting all the 1000 × nsim cases used to derived Figure 2); this depends on each run (implicitly the generated random numbers in each simulation), and sometimes less than 1 × 10 5 simulations are required. For the considered runs, 2 × 10 5 simulations seem to guarantee the obtaining of the values reported in Table 2. In any case, when the muti-linear regression is successfully performed, the results are quite adequate and invariant for increasing number of simulations; this is also the case for the coefficients of the regression (i.e., they remain independent of nsim), which are c = 0.5837, b1 = 0.0998, and b2 = −0.1333.  Figure 2 shows that is not always possible to compute the statistics for the MCS; this happens when not a single failure is reported in one or more of the 1000 runs for a given n sim (at least 5 × 10 4 are necessary). This is not a problem when using the fitted polynomials. Additionally, the MCS approach always tend to larger mean errors and standard deviations for a decreasing number of simulations. This makes the normality polynomial approach more adequate for estimating the reliability indices; however, for few simulations the errors are too large. Nevertheless, the designer could decide which precision level (quantitatively) is willing to accept using information like the one in Figure 2 as an aid (and reduce the number of required simulations as a function of such an accepted error).
If desired, the curves in Figure 2 could be casted as a mathematical expression. For instance, the following power equations fit very well the mean (µ ε ) and standard deviation (σ ε ) of the error in Figure 2 (by fitting from 2.5 × 10 2 simulations and over) If the power in Equations (19) and (20) is assumed as −0.5 in both expressions, the coefficient of variation of the error is ν ε ≈ 150/250 ≈ 0.6, i.e., it is constant and roughly independent of the number of simulations; the actual ν ε obtained from the 1000 runs does exhibit such a roughly constant behavior for this case, except that it is ≈0.7 (difference related to the actual different powers in Equations (19) and (20)). Power equations like Equations (19) and (20) can be linearized by taking logarithms on both sides. Therefore, if the involved variables are transformed into the logarithmic space, a linear fitting can be performed. In this study, we simply used a built-in function in the commercial software for the fitting.
Coefficients for the fitted normality polynomial in Figure 1a are shown in Table 1 (upper set of values) for selected values of n sim . The computing of these coefficients is based on minimizing the error in Equation (4) and was implemented in the coded program by using a built-in function of the programming language employed (MATLAB). As mentioned before, the coefficients a 0 can be linked to the generalized reliability index. Coefficients for results in Figure 1b-d (or those associated to Figure 2) were computed but not shown for brevity. As regards the reliability index, design point and sensitivity factors derived from the simulations (i.e., those obtained with Equations (11), (12), and (14)), they are compared with those obtained by applying the FORM to Equation (15). They are summarized in Table 2 for selected values of n sim . Table 2. Design point, β and α i by using multi-linear regression and first order reliability method (FORM). Results listed in Table 2 indicate that when the points obtained by applying the criterion in Equation (8) were enough to successfully perform a multi-linear regression (also with a built-in function, as in the case of the normality polynomial fitting), β did not deviate from the exact value, but marginally, just like the reliability index obtained with FORM. However, a relatively large n sim was required, usually at least 1 × 10 5 simulations (by inspecting all the 1000 × n sim cases used to derived Figure 2); this depends on each run (implicitly the generated random numbers in each simulation), and sometimes less than 1 × 10 5 simulations are required. For the considered runs, 2 × 10 5 simulations seem to guarantee the obtaining of the values reported in Table 2. In any case, when the muti-linear regression is successfully performed, the results are quite adequate and invariant for increasing number of simulations; this is also the case for the coefficients of the regression (i.e., they remain independent of n sim ), which are c = 0.5837, b 1 = 0.0998, and b 2 = −0.1333.

Design Point Sensitivity Factors
If the tolerance e l in Equation (8) is increased, the minimum number of required simulations can be decreased. For instance, if e = 0.25 were used (instead of the actual used e = 0.05), 5 × 10 4 simulations would be enough for a successful multi-linear regression. Moreover, the design point, sensitivity factors and reliability index would be the same as those reported in Table 2. The opposite would occur if e = 0.005 were used (instead of the actual used e = 0.05), i.e., a much larger number of simulations would be required to determine the reliability parameters from the regression. Therefore, this approach based on multi-linear regression can be quite an adequate alternative by itself to compute β, conditioned on the feasibility of performing enough simulations, the number of which can be decreased by using a large e; with the additional advantage that the design point and sensitivity factors can also be determined.
In the following section a more realistic LSF for a structural application is used to further investigate the revisited methods.

Reliability of Reinforced Concrete Beam under Flexure Moment
In this section the approaches described before are applied to a reinforced concrete beam (RCB) subjected to flexure moment. The example is the same investigated in a previous study [23] but focused only in one design code [24] and three ratios of the mean live to the mean dead load effect for the beam. The rectangular beam section information, LSF, and statistics are succinctly reproduced below. The LSF is where B is the modeling error, f' c is the concrete compressive strength, A s is the reinforcement steel area, f y is the yielding stress of the reinforcement steel, b is the section width, h is the effective depth, and D and V are the dead load and live load effect, respectively (flexure moment). The information of all the independent variables in Equation (21) is summarized in Table 3. A s is assumed deterministic and equal to 3000 mm 2 . The PDFs of the random variables in Table 3 are based on previous literature, which in turn reflects results from experimental projects, field information, observed phenomena, and even the engineers experience to characterized these variables properly, since such PDFs have a direct impact in the computed reliabilities, code calibration tasks, and ultimately in the safety of real structures. More details can be found in [23] and the references therein. Mean values of D and V are not defined in Table 3, but they are derived by considering given mean live load effect (m V ) to mean dead load effect (m D ) ratios (r V/D = 0.4, 1.0 and 2.0 are used in this study) and the assumption that the RCB just meets the code requirement; thus, the following expression is used to determine the mean values.
where m denotes the mean values of the variables in the corresponding subscripts, and Φ = 0.9. Using the previous information, the normality polynomial approach is applied to Equation  (18) β bench as the average of the 1000 runs for n sim = 2 × 10 7 (assumed as the exact value). MCS results are depicted with dashed-doted lines. To perform the MCS, m D and m V are defined using Equation (22) and r V/D as mentioned before. Once they are determined, MCS can then be performed to obtain samples of D and V (and all other random variables in Table 3) and the probability of failure as per Equation (17) can be computed. As an example of the simulated bending moments, histograms of D and V are shown in Appendix A ( Figure A1) for 1 × 10 6 MCS and r V/D = 1.0; it can be observed that the values are comparable in average (because r V/D = 1.0 is considered), and that histograms for D and V clearly resemble normal and Gumbel distributions, respectively, which is expected given that these variables were sampled from such PDFs.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 21 The error for β obtained with the normality polynomial approach exhibits an asymptotic behavior towards approximately με = 1% for large nsim. As before, power laws fit adequately the error and its uncertainty and are defined as where δ = 656, 1103, and 1349, γ = 0.7401, 0.8493, and 0.9062, κ = 129.8, 201, and 179.1, τ = 0.4937, 0.5609, and 0.5703 for rV/D = 0.4, 1.0, and 2.0, respectively. In Equation (23) the constant unity is included to shift the curve upwards to reproduce the asymptotic behavior mentioned; nonetheless, it could be skipped, and the equations will still fairly adequately describe the mean error. The fitting for με was performed for the whole range of nsim, while for σε, nsim from 250 and over was employed. Note that although the range for the fitting could be established based on practical grounds and fitting improvement, in any case the errors and their uncertainties follow a power law; this is the case for the three case studies carried out in this study. It is noteworthy that the normality polynomial approach leads to comparable με and σε for Figures 2 and 4, considering that the LSF for the RCB is a more complex (non-linear) function, and that it has much more random variables and several PDFs.   Table 1 (middle set of values). If normality polynomials of order higher than 3 are used, no further accuracy is gained (or even higher inaccuracies could be obtained; [11]). This is confirmed by carrying out a single case for rV/D = 0.4 using a 4th-order polynomial, since the results are comparable to those of the 3rd-order polynomial case ( Table 1, lower set of values). Note that the order of the coefficients of the normality polynomials for the RCB problem can be very small (compared to the classical LSF problem and to a coastal engineering application shown later); this  Figures 1 and 2 can be extracted. Some additional observations worth to mention are that n sim = 1 × 10 4 seems a reasonable number of simulations for the normality polynomial approach, if a compromise between n sim and error (in terms of mean and standard deviation) is envisaged, that the fitted polynomials lead to better results than the FORM for increasing n sim and that larger errors are obtained for smaller r V/D ; this latter aspect could be attributed to a better approximation of the failure surface for larger r V/D , since the first order approximation to the failure surface by the other method shown in Figure 3 (i.e., the FORM), also deviates more from the exact value for decreasing r V/D . By observing Figure 4, it is pointed out once more that error and its uncertainty is less for normality polynomial when decreasing n sim than for MCS for this case too and, as previously mentioned, is not always possible to estimate the error for MCS for decreasing number of simulations. The error for β obtained with the normality polynomial approach exhibits an asymptotic behavior towards approximately µ ε = 1% for large n sim . As before, power laws fit adequately the error and its uncertainty and are defined as  (23) the constant unity is included to shift the curve upwards to reproduce the asymptotic behavior mentioned; nonetheless, it could be skipped, and the equations will still fairly adequately describe the mean error. The fitting for µ ε was performed for the whole range of n sim , while for σ ε , n sim from 250 and over was employed. Note that although the range for the fitting could be established based on practical grounds and fitting improvement, in any case the errors and their uncertainties follow a power law; this is the case for the three case studies carried out in this study. It is noteworthy that the normality polynomial approach leads to comparable µ ε and σ ε for Figures 2 and 4, considering that the LSF for the RCB is a more complex (non-linear) function, and that it has much more random variables and several PDFs.
The fitted coefficients of the polynomials for Figure 3 (corresponding to r V/D = 0.4) and selected n sim are listed in Table 1 (middle set of values). If normality polynomials of order higher than 3 are used, no further accuracy is gained (or even higher inaccuracies could be obtained; [11]). This is confirmed by carrying out a single case for r V/D = 0.4 using a 4th-order polynomial, since the results are comparable to those of the 3rd-order polynomial case ( Table 1, lower set of values). Note that the order of the coefficients of the normality polynomials for the RCB problem can be very small (compared to the classical LSF problem and to a coastal engineering application shown later); this could be due to the units employed, and should not be understood as if the order of the polynomial could be decreased, whereas obtaining comparable precision, because the use of at least third-order normality polynomials was illustrated and found adequate in [11].
To end this section, the results of using the multi-linear regression approach for the LSF defined in Equation (21) are listed in Table 4 for r V/D = 1.0. The subscripts in Table 4 (and the units of the design point) are associated to the random variables in Table 3. The reported values correspond to the last of the 1000 runs used to develop Figure 4. As an example of the coefficients obtained by multi-linear regression, the ones from the last of the 1000 runs (for deriving Figure 4 Table 4. Design point, β and α i by using multi-linear regression and FORM for the RCB. The previous information indicates that a similar conclusion to that of the previous example (i.e., for the case of Equation (15)) can be drawn, i.e., at least a sufficiently large number of simulations is required for a successful multi-linear regression. However, unlike in the previous example, the design point and sensitivity factors are not invariant by varying n sim . The differences are not so significant though; therefore, once a minimum number of simulations is ensured (around 8 × 10 4 simulations) a very precise β is obtained; it is also observed that the required number of simulations for adequately carrying out the multi-linear regressions decreases with increasing r V/D (this could be attributed to the same reason argued before about the larger errors obtained for smaller r V/D ). If the tolerance e is increased, the number of simulations can be reduced, but not as significantly as for the classical LSF case (i.e., Equation (15)). For instance, an increment to e = 0.35 reduces n sim to around 5×10 4 ; this also changes the values of the design point and sensitivity factors, but not substantially. From Table 4, it is also observed that the values are in very good agreement with the FORM results, with even higher precision from the regression approach for the reliability index. Therefore, it is concluded that the multi-linear regression by itself can be a very attractive alternative to compute β, if a minimum n sim (similar to those mentioned above) is feasible; it is emphasized once more that an additional advantage is that the design point and sensitivity factors are also determined.

Design Point/(Sensitivity Factors)
One final application of the revisited described methods is performed for a coastal structure in the following section.

Overtopping Reliability of a Breakwater
In this section, we consider for the coastal engineering application the example reported in [10], where certain conditions are assumed and where the reader is referred to for further details and used references. It is a breakwater with deterministic slope, tan τ = 1/1.5, and freeboard, F b = 10 m. When the water runs up the breakwater, overtopping could occur (i.e., the water surpasses the freeboard), which is considered as a failure. This is defined by the LSF given by where A u and B u are coefficients characterized as independent normally distributed random variables, with mean values equal to 1.05 and −0.67, respectively, and coefficients of variation both equal to 0.2 [10]; H denotes de wave height and T represents wave period. H and T are random variables probabilistically characterized by the joint Longuet-Higgins distribution [20] with parameter ν = 0.25. The joint PDF of the Longuet-Higgins distribution is given by where H n = H/H s and T n = T/T z are normalized wave heights and periods by considering H s = 5 m and T z = 10 s, which are the significant wave height and the zero up-crossing mean period, respectively, which define the sea state [10]; L(ν) is a normalization factor implying only positive values of T n and defined by First, the FORM is applied to Equation (25). Salient points of performing the FORM to this overtopping LSF are briefly described in the following. First, it is noted that since H n and T n are not independent, the Rosenblatt transformation is performed for the joint distribution to map the equivalent distribution parameters into the normal space by using [17] where Φ −1 (•) denotes the inverse of the CDF of a standard normal variable, F H n (H n ) and F T n |H n (T n |H n ) are the marginal distribution of H n and the conditional distribution of T n given H n for Equation (26), respectively, and defined by where the error function is given by In Equation (29) the equivalent version reported in [25], rather than the original version in [10], is considered. This is so, simply because the error function used in [25] is more readily available in current software. Then, to derive the conditional probability distribution, we divided Equation (26) by Equation (29) yielding Equation (30) given above. Since the CDFs of Equation (29) and Equation (30) are also required to obtain the equivalent parameters mapped in the standardized normal space, other point to highlight is that they were obtained numerically at the design point, unlike for the normal distributed random variables, where simple analytical expressions can be used (which is also possible for other common PDFs).
Additionally, it is noted that as part of the procedure to obtain the reliability index in each iteration of the FORM, usually a vector obtained by multiplying each partial derivative of the LSF (i.e., Equation (25)), evaluated at the design point, by the equivalent second moment in the normal space (for the corresponding random variable) is enough. However, this approach is not possible for the joint random variables in this example. Therefore, the Jacobian (and its inverse) is required [17]; once the inverse of the Jacobian is computed, it is multiplied by the vector of partial derivatives evaluated at the design point mentioned above, and the reliability index can then be obtained in each iteration in the regular way for the FORM (i.e., as when the variables are independent). This approach is followed in the present study. Note that for a set of jointly distributed random variables x i (z i in the normalized space), the inverse of the Jacobian is a lower-triangular matrix determined (often numerically) as [17] where Φ (z i ) is the PDF of a standard normal random variable, with the argument z i obtained in an analogous way to Equation (28); f i and F i refers to the PDF and CDF for the variable with subscript i, respectively. It was noticed that for the present example, disregarding the elements outside the Jacobian diagonal does not impact very significantly the computed reliability indices. A few final important aspects regarding the FORM worth to mention, include that the order of the variables in defining Equation (28) does matter, although similar results may be expected [17]. For instance, in [10] the marginal distribution of T n and the conditional distribution of H n given T n are used to define Equation (28) (i.e., the order of the variables is inverted as compared with this study), which results in a reliability index, β, equal to 2.01 for the problem in question, whereas β = 2.10 is obtained in this study with the FORM formulation described earlier, and adopted in the following; β = 2.10 is also closer to the exact value to be discussed later. Another slight difference between [10] and this study when applying the FORM, is that in the present work, when assuming initial design points, one is determined by setting g bkw = 0, to ensure that the design point is on the failure boundary (e.g., [26]).
To inspect the variation of β for different F c values, the FORM is performed by varying the freeboard between 9 m and 12 m and the resulting reliability index is shown in Figure 5a with a black dashed line. As expected, it can be observed that β increases for increasing freeboard; if the slope of the breakwater is increased to tan τ = 1/2 and the FORM is carried out for the same range of F c , it further increases reliability levels, as shown by the dashed grey line in Figure 5a. These results are used as reference and for comparison purposes, with respect to the results from the normality polynomial and multi-linear regression approaches revisited in this study.  The simulations for this coastal engineering application, used as the basis of the revisited methods, are much more computationally intensive than for the classical and structural examples, because of the dependency between the wave height and period and inclusion of the Longuet-Higgins distribution, which imposes numerical computing for the probability levels (e.g., values from CDFs) and a different method for the sampling. This latter aspect, i.e., the generation of jointly distributed random numbers when a set of xi variables are dependent, is based on expressing the joint PDF as [22] with the corresponding CDF given by   The simulations for this coastal engineering application, used as the basis of the revisited methods, are much more computationally intensive than for the classical and structural examples, because of the dependency between the wave height and period and inclusion of the Longuet-Higgins distribution, which imposes numerical computing for the probability levels (e.g., values from CDFs) and a different method for the sampling. This latter aspect, i.e., the generation of jointly distributed random numbers when a set of x i variables are dependent, is based on expressing the joint PDF as [22]

Freeboard (m)
with the corresponding CDF given by Using the previous concepts, and considering a set of values U generated from n independent standard uniformly distributed random variables, the set of dependent random variables can be determined as where F −1 (•) denotes the inverse of the CDF. The obtaining of this inverse of the CDF can be relatively straightforward for some common probability distributions, where an analytical expression can be used for Equation (35). This is not the case for the Longuet-Higgins distribution. In this case, the jointly distributed random wave height and period must be determined numerically. Figure 6 shows samples of jointly generated random values of wave height and period in the normalized space (for n sim = 1 × 10 3 , 5 × 10 3 ; 1 × 10 4 and 5 × 10 4 ). A few contours of the theoretical Longuet-Higgins distribution (i.e., Equation (26)) are also shown in Figure 6; it can be observed that they are in good agreement. The values in the non-normalized space can be obtained simply by recognizing that H n = H/H s and T n = T/T z .
Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 21 jointly distributed random wave height and period must be determined numerically. Figure 6 shows samples of jointly generated random values of wave height and period in the normalized space (for nsim = 1 × 10 3 , 5 × 10 3 ; 1 × 10 4 and 5 × 10 4 ). A few contours of the theoretical Longuet-Higgins distribution (i.e., Equation (26)) are also shown in Figure 6; it can be observed that they are in good agreement. The values in the non-normalized space can be obtained simply by recognizing that Hn = H/Hs and Tn = T/Tz. As mentioned before, the sampling procedure is significantly more time-consuming than for the LSFs in previous sections. Therefore, MCS are sampled only up to 1 × 10 6 simulations for all the variables of the LSF represented by Equation (25), and the reliability index is determined by employing Equation (17) and Equation (7), only for the case reported in [10] (i.e., Fc = 10 m and tan τ =1/1.5).Nonetheless, results depicted in Figure 5b indicate that the reliability index stabilizes, from approximately nsim = 1 × 10 5 and over, to a reliability index practically equal to 2.2 (gray solid line). Therefore, β = 2.2 is adopted as the exact value of the reliability index for the breakwater under overtopping. This value is to be used to assess the error by estimating the reliability index with the normality polynomial approach, and to compare versus the results obtained with the multi-linear regression approach. In fact, the results from these two approaches are also shown in Figure 5b (black As mentioned before, the sampling procedure is significantly more time-consuming than for the LSFs in previous sections. Therefore, MCS are sampled only up to 1 × 10 6 simulations for all the variables of the LSF represented by Equation (25), and the reliability index is determined by employing Equation (17) and Equation (7), only for the case reported in [10] (i.e., F c = 10 m and tan τ =1/1.5). Nonetheless, results depicted in Figure 5b indicate that the reliability index stabilizes, from approximately n sim = 1 × 10 5 and over, to a reliability index practically equal to 2.2 (gray solid line). Therefore, β = 2.2 is adopted as the exact value of the reliability index for the breakwater under overtopping. This value is to be used to assess the error by estimating the reliability index with the normality polynomial approach, and to compare versus the results obtained with the multi-linear regression approach. In fact, the results from these two approaches are also shown in Figure 5b (black solid line for the normality polynomial; grey dashed line for the multi-linear regression approach), where it is observed that the normality polynomial approach converges to a stable value (approximately β = 2.12) from about 2 × 10 3 simulations on, leading in average to a slightly smaller reliability index (i.e., in the conservative side) but closer to the exact value than by using the FORM. The multi-linear regression approach (like the MCS and unlike the normality polynomial) requires a minimum number of simulations to be carried out, being this number 1 × 10 3 for Figure 5b, but sometimes more simulations are required; nevertheless, when a sufficient large number of simulations is performed (e.g., about 3 × 10 4 or more in Figure 5b), the results of the multi-linear regression leads to practically the exact β, and the design point and sensitivity factors can also be determined. These values of sensitivity factors and design points are also very similar with those reported in [10].
It is noted that Figure 5b corresponds to only one set of simulations for every n sim , which may vary for different sets of generated random numbers (as shown in Figures 1 and 3), implying an uncertainty in the deviation from the exact value for different number of simulations. This uncertainty is assessed as for the classical and structural LFSs in previous sections, i.e., by computing the errors in the reliability index as per Equation (18) and fitting them to power laws with the mathematical functional form represented in Equations (19) and (20) (or Equations (23) and (24)), but with different values of the parameters. To do so, and unlike the case of the classical LSF and the reinforced concrete beam under flexure moment, not 1000 but only 100 sets of simulations are computed for each n sim , due to the more extensive required time and computational resources referred to earlier (A comparison in terms of computing time (CPU time), a description and a discussion are given in Appendix B and Figure A2 of the appendix). This leads to the mean errors shown in Figure 5c with a black solid line for the normality polynomial case (including mean values ± one standard deviation indicated in black dashed lines), and with a grey dashed line for the MCS case (including mean values ± one standard deviation indicated in grey dotted lines).
Even though errors reported in Figure 5c exhibit not as a smooth behavior as those observed in Figures 2 and 4 (obtained in an analogous way but for 1000 sets of n sim ), the qualitative trend is fairly similar, especially for mean values and not so small n sim . Indeed, power laws can be adequately fitted to µ ε and σ ε , as shown in Figure 5d by fitting the computed errors from 1 × 10 2 simulations and over; the mathematical functional form is analogous to that of Equations (23) and (24), except that the constant 1 in Equation (23) is omitted. The obtained fitted parameters are δ = 23.62, γ = 0.2342, κ = 47.94 and τ = 0.4254. As observed in Figure 5d the fitting is very adequate for σ ε and adequate for µ ε , albeit only 100 sets of n sim were employed for the statistics.
From Figure 5c,d similar conclusions to those found before can be drawn, namely, that for decreasing n sim the MCS tends to deviate more from the exact value than the normality polynomials (in terms of ε), that for decreasing number of simulations the error for the MCS can be unknown and that power laws are adequate to mathematically defined µ ε and σ ε for the normality polynomial approach. Therefore, a designer could for instance use the normality polynomial method to compute the reliability index for a reduced number of simulations, whereas accepting an error in the estimation. However, such an error could be estimated if expressions of µ ε and σ ε (like the power laws determined in this study) are known.
As an example, in Figure 5a n sim = 7 × 10 2 is used for the normality polynomial approach (black solid line) and, as it is shown, this leads to reasonably adequate results (using a fairly small number of simulations) when compared with the FORM and the MCS (also included in Figure 5a with a black dotted line). Moreover, the fitted equations shown in Figure 5d can be used to quantitatively compute the associated error and its standard deviation with respect to the exact reliability index, that is µ ε = 5.09% and σ ε = 2.95%. This is strictly applicable only to F c = 10 m; however, comparable errors may be expected for a range of freeboard values by inspecting Figure 5a. Naturally, the contents of this paper could be extended to investigate how the error changes by varying one or more parameters of the LFSs. In such a case, one would expect that functional forms like those reported in this study can be used to assess ε, but possibly with higher mean errors (and/or standard deviations) for higher reliability levels, because usually more simulations are required for lower probabilities of failure. This could be inferred from Figure 5a, where a last set of calculations is shown by increasing the breakwater slope to 1/2 (dashed and dotted grey lines for the normality polynomial and MCS techniques, respectively), where higher variations of the normality polynomial in relation to the FORM are observed; this higher reliability levels also have the effect of decreasing the ability of the MCS to capture the probability of failure, as also observed in Figure 5a for a wide range of F c . Additionally, although not shown in Figure 5a, it was observed that the minimum number of simulations required to adequately performed the multi-linear regression increases for higher reliability levels (e.g., larger breakwater slopes).
Overall, results in this section indicate that the revisited simulation-based methods can be also effective for coastal engineering applications.

Discussion
Results in the previous section suggest that the two revisited approaches based on simulations, namely, the normality polynomial and the multi-linear regression approaches, are effective in reducing the number of required simulations while adequately computing the reliability index, design point and sensitivity factors. It could be argued that still a relatively large number of simulations are required. However, the computing power is becoming higher every year, and these methods proposed at the end of the 1990s could become a feasible alternative for some complex models two or three decades later.
The power law (Equations (19), (20), (23), and (24)), which describes the precision in computing β for the normality polynomial approach, was found very adequate for the three LSFs considered. Albeit it cannot be concluded that the underlaying error is based on the power law for every possible LSF, since one of LSFs studied here was a simple classical case using only one type of PDF, and the other LSFs were more complex (non-linear functions), with much more random variables and included several PDFs, as well as dependency between variables and the joint Longuet-Higgins distribution, it is reasonable to believe that the error for other LSFs for a wide range of coastal and structural engineering applications could follow the power law. A designer could opt to reduce the number of simulations while accepting an error level (including its uncertainty) by using the power laws as an aid.
The multi-linear regression approach was originally developed to derive the design point and sensitivity factors not obtained when performing MCS; however, it is considered that it can be an alternative by itself to compute β in an accurate way, conditioned on performing enough simulations for a successful regression; the number of simulations can be reduced by increasing the tolerance e. Values reported in this study can be used as a guide.
It is acknowledged that some differences with respect to the present study could be found when other LSFs and applications are used. However, the values reported in this study could be used for guidance, and it is believed that the power law may hold in many coastal and structural engineering applications, since the normality polynomials are based in strong mathematical foundations, as referenced before; nevertheless, future research to further inspect the findings in this study is recommended by using mathematical LSFs considered as benchmark in the literature, but also more ultimate and serviceability LSFs for other coastal and structural engineering applications. It is also believed that the revisited methods and the findings in this study can be exported to other engineering fields if practical applications can be posed as a capacity-demand problem and when extensive simulations are required, including system reliability (for instance for reinforced concrete frame buildings, among many other possibilities). Future research could also include a systematic study for the multi-linear regression approach by varying the tolerance, for given number of simulations, so that the number of MCS can be limited to a minimum while guaranteeing the obtaining of adequate reliability-related values.
If more LSFs are investigated in future studies, perhaps it could be possible to infer general bounds for a wider applicability of the findings in the present study.

Conclusions
Two reliability methods based on simulations are revisited. One method fits normality polynomials to the simulated data with fractile constraints, and the other approximates the linearized limit state surface at the design point using multi-linear regression; for the latter, a slight modification is proposed. Three limit state functions, a very simple one, other for a structural engineering application and another for a coastal engineering application, are employed.
The most relevant findings of this study are that for the normality polynomial approach, a power law was found to adequately represent the mean and standard deviation of the error in the estimated reliability index as a function of the number of simulations. It could be used as an aid for decision makers to select a precision level (quantitatively) associated to a selected reliability index, thus reducing the number of required simulations by expressively accepting an error level. Additionally, it is found that the multi-linear regression approach is an excellent option to obtained accurate reliability levels, although a sufficiently large number of simulations is required (not prohibitive though). It also has the advantage that the design point and sensitivity factors are determined.
Other findings in this study are: 1.
When the normality polynomial approach is used, the reliability index is dependent on the generated random numbers of each run, but it becomes stable for a large number of simulations.

2.
The reliability index cannot always be determined with the Monte Carlo simulations, while the opposite occurs when normality polynomials are used, although significant deviations from the exact value are observed for small numbers of simulations. In general, for an intermediate number of simulations (e.g., 1 × 10 3 ), the fitted normality polynomials lead to a better estimate of the reliability index than the Monte Carlo simulations. 3.
When the mean relative error and its standard deviation are computed for the reliability index (compared to the exact value), for decreasing number of simulations the Monte Carlo simulation approach tend to larger mean errors and standard deviations than the normality polynomial approach. 4. 3rd-order normality polynomials were mostly used; when 4th order ones are used, the fitting leads to comparable results.

5.
When the multi-linear regression approach is considered, a minimum number of simulations is required for successfully performing the regression (in the order of 10 4 to 10 5 simulations), but once this is ensured, a very precise reliability index is obtained (more precise than by using the first order reliability method (FORM)), and the design point and sensitivity factors are also determined and in good agreement with those determined with the FORM. 6.
If the tolerance for the multi-linear regression approach is increased (i.e., if a wider range in the nearest of the failure surface is stipulated to gather the vectors of simulated data), the number of simulations can be reduced. Funding: The financial support from the Erasmus Mundus Coastal and Marine Engineering and Management (CoMEM) programme for one of the authors of this study and from Universidad de Guanajuato (División de Ingenierías and Campus Guanajuato) is gratefully acknowledged. Figure A2 shows the computing (CPU) time required for the classical and structural engineering LSFs ( Figure A2a) and for the coastal engineering LSF ( Figure A2b). The CPU times include computing of MCS and probabilities of failure, fitting of the normality polynomials and obtaining of the reliability parameters by multi-linear regression. The employed processor is an Intel(R) Core(TM) i7-9750H CPU @ 2.60 GHz with RAM of 16.0 GB and operating system of 64 bits.
It is highlighted that CPU times are for 1000 sets of nsim and for 100 sets of nsim (as indicated in the horizontal axes of the figure) for the classical and RCB LFSs and overtopping LSF, respectively. It can be observed in Figure A2 that the CPU times are significantly larger for the coastal engineering application, as shown by the different ranges used in the vertical and horizontal axes and by the fact that only 100 sets of nsim are used for this case (compared to 1000 sets for the others, as mentioned before). This significant larger computing time is imposed by the joint distribution of wave heights and periods and the Longuet-Higgins distribution used to represent them, which must be solved numerically and for which a different sampling technique is required, as indicated in the main body of this article. Figure A2 could assist the readers to establish feasible simulation schemes. It was noticed that efficient computing time is obtained if 10,000 simulations at a time are considered for the overtopping LSF. See for instance that for nsim = 100 in Figure A2b (which translates into 100 × 100 = 10,000 MCS), a reasonable CPU time is required; in fact, it is shown in Figure A2b that over this threshold the CPU time starts to increase to a much faster rate. Therefore, once the random numbers are simulated (this  Figure A2 shows the computing (CPU) time required for the classical and structural engineering LSFs ( Figure A2a) and for the coastal engineering LSF ( Figure A2b). The CPU times include computing of MCS and probabilities of failure, fitting of the normality polynomials and obtaining of the reliability parameters by multi-linear regression. The employed processor is an Intel(R) Core(TM) i7-9750H CPU @ 2.60 GHz with RAM of 16.0 GB and operating system of 64 bits.

Appendix B
It is highlighted that CPU times are for 1000 sets of n sim and for 100 sets of n sim (as indicated in the horizontal axes of the figure) for the classical and RCB LFSs and overtopping LSF, respectively. It can be observed in Figure A2 that the CPU times are significantly larger for the coastal engineering application, as shown by the different ranges used in the vertical and horizontal axes and by the fact that only 100 sets of n sim are used for this case (compared to 1000 sets for the others, as mentioned before). This significant larger computing time is imposed by the joint distribution of wave heights and periods and the Longuet-Higgins distribution used to represent them, which must be solved numerically and for which a different sampling technique is required, as indicated in the main body of this article. Figure A2 could assist the readers to establish feasible simulation schemes. It was noticed that efficient computing time is obtained if 10,000 simulations at a time are considered for the overtopping LSF. See for instance that for n sim = 100 in Figure A2b (which translates into 100 × 100 = 10,000 MCS), a reasonable CPU time is required; in fact, it is shown in Figure A2b that over this threshold the CPU time starts to increase to a much faster rate. Therefore, once the random numbers are simulated (this is not a problem in terms of CPU times for millions of random numbers), a programming scheme dividing the computing in 10,000 MCS can be used to improve the efficiency (e.g., subdividing the tasks within the same program, running several windows simultaneously and/or using several computers).