{'fraction_non_alphanumeric': 0.02855497595013888, 'fraction_numerical': 0.009924801842693584, 'mean_word_length': 4.591477272727273, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multilevel Monte Carlo for uncertainty quantification in structural engineering

Practical structural engineering problems often exhibit a significant degree of uncertainty in the material properties being used, the dimensions of the modeled structures, etc. In this paper, we consider a cantilever beam and a beam clamped at both ends, both subjected to a static and a dynamic load. The material uncertainty resides in the Young's modulus, which is modeled by means of one random variable, sampled from a univariate Gamma distribution, or with multiple random variables, sampled from a Gamma random field. Three different responses are considered: the static elastic, the dynamic elastic and the static elastoplastic response. In the first two cases, we simulate the spatial displacement of a concrete beam and its frequency response in the elastic domain. The third case simulates the spatial displacement of a steel beam in the elastoplastic domain. In order to compute the statistical quantities of the static deflection and frequency response function, Multilevel Monte Carlo (MLMC) is combined with a Finite Element solver. In this paper, the computational costs and run times of the MLMC method are compared with those of the classical Monte Carlo method, demonstrating a significant speedup of up to several orders of magnitude for the studied cases.


Introduction
There is an increasing need to accurately simulate and compute solutions to engineering problems whilst taking into account model uncertainties. Methods for such uncertainty quantification and propagation in structural engineering can roughly be categorized into two groups: non-sampling methods and sampling methods. Examples of non-sampling methods are the perturbation method and the Stochastic Galerkin Finite Element method. The perturbation method is based on a Taylor series expansion approximating the mean and variance of the solution [1]. The method is quite effective, but its use is restricted to models with a limited number of relatively small uncertainties. The Stochastic Galerkin method, first proposed by Ghanem and Spanos [2], is based on a spectral representation in the stochastic space. It transforms the uncertain coefficient partial differential equation (PDE) problem by means of a Galerkin projection technique into a large coupled system of deterministic PDEs. This method allows for somewhat larger numbers of uncertainties and is quite accurate. However, it is highly intrusive and memory demanding, making its implementation cumbersome and restricting its use to rather low stochastic dimensions.
Sampling methods, on the other hand, are typically non-intrusive. Each sample corresponds to a deterministic solve for a set of specified parameter values. Two particularly popular examples are the Stochastic Collocation method [3] and the Monte Carlo (MC) method [4]. The former samples a stochastic PDE at a carefully selected multidimensional set of collocation points. After this sampling, a Lagrange interpolation is performed leading to a polynomial response surface. From this, the relevant stochastic characteristics can easily be computed in a post-processing step. However, as is also the case for Stochastic Galerkin, the Stochastic Collocation method suffers from the curse of dimensionality; the computational cost grows exponentially with the number of random variables considered in the problem. The MC method on the other hand, selects its samples randomly and does not suffer from the curse of dimensionality. A drawback is its slow convergence as a function of the number of samples taken. The convergence of Monte Carlo can be accelerated in a variety of ways. For example, alternative non-random selections of sampling points can be used, as in Quasi-Monte Carlo [5,6], and Latin Hypercube [7] sampling methods. Also, variance reduction techniques such as Multilevel Monte Carlo (MLMC) [8] and its generalizations, see, e.g., [9,10], can speed up the method. Note that there also exist hybrid variants which exhibit both a sampling and non-sampling character. This type of methods combines, for example, the Stochastic Finite Element methodology with Monte Carlo sampling or a multi-dimensional cubature method, see, e.g., [11,12].
Monte Carlo methods have since long been used in the field of structural engineering, for example in problems of structural dynamics [13] or in elastoplastic problems where the structure's reliability is assessed [14]. The focus of this paper is to apply the MLMC method to a structural engineering problem, discretized by means of the Finite Element method. A comparison with the standard MC method will show a significant reduction in computational cost.
The structure of the paper is as follows. In section 2, we formulate the mathematical model, introduce the problem statement and describe how the uncertainty is modeled. Section 3 recalls the MLMC method, and provides some additional algorithmic implementation details. In Section 4, numerical results are presented. First, a performance comparison is made between standard MC and MLMC for a static elastic and elastoplastic load case. Next, results showing the uncertainty propagation for the static elastic, static elastoplastic and the dynamic elastic response are shown for when the Young's modulus is modeled by means of a random variable coming from a univariate Gamma distribution and as multiple random variables sampled from a Gamma random field. The fifth and last section offers concluding remarks and details some paths for further research.

Beam models and material parameters
The considered engineering problem is the response of a cantilever beam clamped at one side and a beam clamped at both sides as seen in Fig.1, assuming plane stress. The following responses are considered: the static elastic response, the dynamic elastic response and the static elastoplastic response. The first two respectively simulate the spatial displacement and the frequency response function of a concrete beam model solely in the elastic region. The third one simulates the spatial displacement of a steel beam in the elastic as well as in the plastic region. The dimensions of the beam are 2.5 m (length) by 0.25 m (height) by 1 m (width) for the elastic responses and 10 −3 m (width) for the elastoplastic response. The material parameters of the concrete are as follows: a mass density of 2500 kg/m 3 , a Poisson ratio of 0.15 and a Young's modulus subjected to some uncertainty with a mean value of 30 GPa. The material parameters of the steel are as follows: a yield strength of 240 MPa, a Poisson ratio of 0.25 and a Young's modulus subjected to a certain degree of uncertainty with a mean value of 200 GPa. In order to model the material uncertainty, two uncertainty models will be considered. The first model consists of a homogeneous Young's modulus characterized by means of a single random variable. The second model is a heterogeneous Young's modulus represented as a random field. Both uncertainty models will be used to compute the stochastic characteristics of the static elastic and elastoplastic response of a beam clamped at both sides as seen in Fig.1 (right) and the dynamic response of a cantilever beam clamped at one side as seen in Fig.1 (left).

The homogeneous model
Following [15], we opt to describe the Young's modulus in the homogeneous model by means of a univariate Gamma distribution. This distribution is characterized by a shape parameter α and a scale parameter β: The corresponding mean value and variance can be computed respectively as µ = αβ and σ 2 = αβ 2 . In this paper, we will select α = 7.1633 and β = 4.1880 × 10 9 in order to model the material uncertainty for the concrete beam, which are based on values coming from [16]. This leads to a mean of 30 GPa and a standard deviation of 11.2 GPa. For modeling the material uncertainty in the steel beam, we select α = 934.2 and β = 0.214 × 10 9 , see [17]. This gives a mean of 200 GPa and a standard deviation of 6.543 GPa. The Gamma distribution for both materials is plotted in Fig. 2  Shown also is the mean µ and the standard deviation σ.

The heterogeneous model
The Young's modulus with spatially varying uncertainty will be represented by means of a (truncated) Gamma random field. The construction of this random field is done by means of a classical, two step process.
First, a (truncated) Gaussian random field is generated, using a Karhunen-Loève (KL) expansion [18]. Next, this Gaussian random field is transformed into a Gamma random field with a memoryless transformation [19].
Consider a Gaussian random field Z(x, ω), where ω is a random variable, with exponential covariance kernel, We selected the 1-norm (p = 1), a correlation length λ = 0.3 and a standard deviation σ = 1.0. The corresponding KL expansion can then be formulated as follows: Z(x, .) denotes the mean of the field, and is set to zero. The ξ n (ω) denote i.i.d. standard normal random variables. The symbols θ n and b n (x) denote, respectively, the eigenvalues and eigenfunctions of the covariance kernel (2). For the parameter values selected above, those eigenvalues and eigenfunctions can be computed analytically, see [2,20,21,22]. For a one-dimensional domain [0, 1], they are given by The normalizing constants A n are chosen such that b n 2 = 1. The constants w n represent the real solutions, in increasing order, of the transcendental equation For the two-dimensional case, the eigenvalues and functions are obtained in a tensorproduct way, In an actual implementation, the number of KL-terms in (3) is truncated to a finite value. This value depends on the magnitude and, more precisely, on the decay of successive eigenvalues. Those eigenvalues are plotted in Fig. 3 (left). Inclusion of the first 101 KL-terms is sufficient in order to represent 90% of the variance of the random field. The percentage of the variance that is accounted for as a function of the number of included eigenvalues, corresponds to the cumulative sum of the eigenvalues, and is also illustrated in the figure. A cumulative sum of 1.0 corresponds to 100% of the variance of the field being accounted for. Normalized PDF value Once the Gaussian field has been generated, a memoryless transformation is applied pointwise, in order to obtain the Gamma random field [19]. Here, F denotes the marginal cumulative density function

The Finite Element method
The Finite Element method will be used to compute the responses of the beam assuming plane stress.
An equidistant, regular rectangular mesh is applied with bilinear quadrilateral elements. The implemented solver for the elastic and the elastoplastic response is different. The underlying equations and solution methods are reviewed hereunder.
For the static elastic response, the system equation is of the form with K the global stiffness matrix, f the global nodal force vector and u the displacement. The global stiffness matrix and nodal force vector are obtained from the element stiffness matrices K e and the element force vectors f e . These are computed analytically by evaluation of the following integrals: The element nodal force vector f e is modeled as a Neumann boundary condition, where t n stands for the surface traction specified as a force per unit area and N the element shape function matrix, integrated over the free element boundary Γ t . The element stiffness matrix K e is obtained by integrating the matrix B T DB over the element's surface Ω. Matrix B is defined as LN with L the derivative matrix specified below, and with D the elastic constitutive matrix for plane stress, containing the element wise material parameters, For the dynamic response case the following equation is obtained, Matrix M denotes the system mass matrix obtained from the assembly of the element mass matrices M e . f denotes the frequency and ρ the volumetric mass density of the material. The multiplication of the system stiffness matrix K with the imaginary unit ı and the constant η, denotes the damping matrix.
The approach for solving the static elastoplastic case differs due to the nonlinear stress-strain relation in the plastic domain. The plastic region is governed by the von Mises yield criterion with isotropic linear hardening. An incremental load approach is used starting with a force of 0 N. The methods used to solve the elastoplastic problem are based on Chapter 2 §4 and Chapter 7 §3 and §4 of [23]. For this case, the system equation takes the following form: where ∆u stands for the resulting displacement increment. The vector r is the residual, where f stands for the sum of the external force increments applied in the previous steps, ∆f for the applied load increment of the current step and q for the internal force resulting from the stresses First the displacement increment of all the nodes is computed according to (12), with an initial system stiffness matrix K resulting from the assembly of the element stiffness matrix K e , computed by means of a Gauss quadrature where D ep denotes the elastoplastic constitutive matrix. The initial state of D ep is the elastic constitutive matrix from (10). Secondly, the strain increment ∆ε is computed, Thirdly, the nonlinear stress-strain relationship, is integrated by means of a backward Euler method. The backward Euler method essentially acts as an elastic predictor-plastic corrector; an initial stress state that is purely elastic is computed and then projected in the direction of the yield surface as to obtain the plastic stress state. Due to the implicit nature of the integrated stress-strain relation, this equation must be supplemented with the integrated form of the hardening rule and the yield condition. This system of nonlinear equations is then solved with an iterative Newton-Raphson method. Afterwards, the consistent tangent stiffness matrix is computed [24]. This matrix is then used to compute the updated element stiffness matrix (15), resulting in an updated system stiffness matrix K.
The inner iteration step of solving the stress-strain relation and the updated system stiffness matrix is repeated for each outer iteration step which solves equation (12). The outer step consists in balancing the internal forces with the external ones as to satisfy the residual, which in our case equals 10 −4 times the load increment. The procedure used is incremental-iterative, relying on the iterative Newton-Raphson method.
This process is repeated for each load increment.
Each finite element is assigned a value of the Young's modulus. This is accomplished by means of the midpoint approach, i.e., the value is taken constant within each individual element and equal to the value of the realization of the random field at the center point of the element [25].

Method overview
The Multilevel Monte Carlo method (MLMC) is an extension of the standard Monte Carlo (MC) method, see, e.g., [8,26]. The method relies on a clever combination of many computationally cheap low resolution samples and a relatively small number of higher resolution, but computationally more expensive samples.
The application of MLMC to our problem will be based on a hierarchy of nested finite element meshes.
These meshes will be indexed from 0 to L, with 0 indicating the coarsest mesh and L the finest mesh. An example of such a hierarchy is shown in Fig. 5. It is common in the PDE setting to use a geometric relation for the number of degrees of freedom between the different levels. We set the number of finite elements for a mesh at level proportional to 2 d , where d is the dimension of the problem.
Multilevel Monte Carlo, on the other hand, starts from a reformulation of E[P L ] as a telescoping sum. The expected value of the quantity of interest on the finest mesh is expressed as the expected value of the quantity of interest on the coarsest mesh, plus a series of correction terms (or differences): Each term in the right-hand side is then estimated separately by a standard Monte Carlo estimator with N samples, i.e., where Q MLMC with V [x] denoting the variance of a random variable x. The MLMC estimator in (20) can be written as a sum of L + 1 estimators for the expected value of the difference on each level, i.e., where we defined P −1 := 0.
Because of the telescoping sum, the MLMC estimator is an unbiased estimator for the quantity of interest on the finest mesh, i.e., Denoting by V the variance of the difference P − P −1 , the variance of the estimator can be written as In order to ensure that the MSE is below a given tolerance 2 , it is sufficient to enforce that the variance where C denotes the cost to compute a single realization of the difference P −P −1 , subject to the constraint Treating the N as continuous variables, we find Note that if E[P ] → E[P ], then V → 0 as increases. Hence, the number of samples N will be a decreasing function of . This means that most samples will be taken on the coarse mesh, where samples are cheap, and increasingly fewer samples are required on the finer, but more expensive meshes. In practice, the number of samples must be truncated to N , the least integer larger than or equal to N .
Using (27), the total cost of the MLMC estimator can be written as This can be interpreted as follows. When the variance V decreases faster with increasing level than the cost increases, the dominant computational cost is located on the coarsest level. The computational cost is then proportional to V 0 C 0 , which is small because C 0 is small. Conversely, if the variance decreases slower with increasing level than the cost increases, the dominant computational cost will be located on the finest level L, and proportional to V L C L . This quantity is small because V L is small. For comparison, the computational cost of a Monte Carlo simulation that reaches the same accuracy is proportional to V 0 C L .
In our numerical results presented next, we will compare the cost of the MLMC estimator to the cost of a standard MC simulation, both in actual runtime (seconds) and in some normalized cost measure. This cost is chosen such that the time needed to obtain a sample on the coarsest mesh (level zero) is equal to a cost of one unit. When using the standard geometric mesh hierarchy, the cost for a single sample on level is proportional to 2 γ . The factor γ is determined by the efficiency of the solver. The advantage of using this normalized cost is that it decouples the cost of the simulation from the specific computer hardware.
The second term in (21) . Then, we can use the heuristic and check for convergence For completion, we now mention the central MLMC complexity theorem. We refer to [26] for a proof.
Theorem 1. Given the positive constants α, β, γ, c 1 , c 2 , c 3 such that α ≥ 1 2 min (β, γ) and assume that the following conditions hold: Then, there exists a positive constant c 4 such that for any < exp(−1) there exists an L and a sequence {N } L =0 for which the multilevel estimator, Q MLMC L has an MSE ≤ 2 , and

Implementation aspects
The MLMC method is non-intrusive. Only an interface between the Finite Element solver routine and the multilevel routine is necessary. All modules used in the simulation are written in Matlab.
The elastic computations are run in sequential mode, without making use of parallelization. The elastoplastic computations are run in parallel mode because of their significantly larger computational cost. In particular, the computation of the individual samples is parallelized in order to get a speedup. This is possible because of the embarrassingly parallel nature of the MLMC and the MC methods. For the aforementioned configuration, a number of 12 samples can be computed concurrently. For more details on load balancing of MLMC samplers, we refer to [27] Computation of the optimal number of samples per level according to (27)

Numerical Results
In this section, we discuss our numerical experiments with the MLMC method. We consider the static and dynamic response, using both a homogeneous and a heterogeneous uncertain Young's modulus. First a comparison will be made between the MC method and the MLMC method. Secondly, the rates from

Comparison between MC and MLMC
First, we compare the efficiency of MLMC and MC for the static elastic and elastoplastic response of the beam clamped at both sides and loaded in the middle as shown in Fig. 1 (right).  Fig. 1   Number of samples

Elastoplastic Response
For the elastic response, the MLMC simulations are level adaptive. The number of extra levels to be added depends on the condition of the bias, (29).
However, for the elastoplastic response we chose to manually set the maximum level. This level is chosen based on a mesh convergence analysis. This is done because for the tolerances considered in Tab. 4, samples would have to be taken on levels up to 5 or more. This would constitute a considerable time cost, which results from the iterative nature of the Finite Element implementation.
The results of a mesh convergence study are shown in Fig. 9, respectively for the elastoplastic response   Fig. 1   As stated before, the parallelization enables the algorithm to compute 12 samples concurrently. This is highly efficient for levels with a large number of samples, as is the case for coarser levels, and less efficient for finer levels. This is due to the load imbalance in the latter case.

Rate verification
In this section, we estimate the rates α, β and γ from Theorem 1, first for the static elastic response, then for the static elastoplastic response.

Elastic Response
We list the online estimated rates, i.e., the rates estimated during the run of the algorithm, for various tolerances on the RMSE, see Tab. 6, corresponding to the tolerances considered in Fig. 6. Following the results in this table, it is possible to estimate the asymptotic cost of the MLMC estimator, using the different regimes from Eq. (30). Since β > γ in all cases, we expect an optimal cost proportional to −2 . This is indeed what we observed in the numerical experiments in Fig. 6.

RMSE [/]
Homogeneous Heterogeneous  In case where the Young's modulus is homogeneous, Fig. 12 (left), the variance corresponding to the finest level is lower than the variance on the previous levels. This phenomenon can be explained by the number of samples on that level, which is used to estimate the variance. Only 7 samples were taken. Taking a higher number of samples on that level would bring the variance to equal the variance of the other levels.

Elastoplastic Response
As for the elastic response, we list the online rates for various tolerances on the RMSE in Tab Fig. 13 shows the behavior of the expected value and the variance of the quantity of interest P , and of the difference P − P −1 , in case of a RMSE equal to 2.5E-6. Comparing Fig. 12 with Fig. 13, one notices that for the elastoplastic response the mean is not perfectly constant over the levels. This slight increase stems from the number of elements used for the discretization; the displacement is indeed dependent upon the discretization of the mesh, Fig. 9, for the mesh sizes considered here.

Static Elastic Response
In this part we show the displacement of the concrete beam under a static load including uncertainty bounds in the elastic domain. The beam configuration is a beam clamped at both ends, loaded at mid span, as shown in Fig. 1 (right). The load is modeled as a distributed load acting on each of the vertical middle nodes of the beam. The force on each node is such that the sum of these forces equals 10000 kN, regardless of the refinement of the mesh.       Fig. 17 show a cut-through in order to better illustrate how the shades of blue represent the PDF.

Static Elastoplastic Response
The solution is presented as a force deflection curve, Fig. 18. The beam configuration is the same as for the static elastic response but with a different material; a steel beam clamped at both sides loaded at mid span, Fig. 1 (right). The load is modeled in the same way as for the static elastic response. The beam is loaded from 0 N to 13.5 kN in steps of 135 N. The line style and color convention is the same as for the static elastic response. As can be observed, the uncertainty bounds in case of a homogeneous Young's modulus, Fig. 18 (left) are wider and more spread out than in case of a heterogeneous modulus (right). This behavior corroborates with the ones from the static elastic response. In Fig. 19, ten individual samples are shown for a homogeneous Young's modulus (left) and a heterogeneous Young's modulus (right). Comparing these two sets of figures, it shows that, the deflection corresponding to the highest applied force in Fig. 19 (left and right) is less than the deflection shown in Fig. 18 (left and right). This stems from the fact that these ten samples per Young's modulus are generated on the coarsest mesh. These results illustrate even better the power of the MLMC method; even whilst taking most of the samples on coarser meshes, one still gets the result as if the simulation was run on the finest mesh.

Dynamic Elastic Response
For the dynamic elastic response, the solution is a frequency response function (FRF). The beam configuration we consider is a concrete cantilever beam clamped on the left side and loaded on the right side, as shown in Fig. 1 (left), with a dynamic load. The load is modeled in the same way as for the static elastic response with the same magnitude. In order to ensure the correct representation of the dynamic response, and to determine the minimum number of finite elements required, the bending wavelength, λ min , must be evaluated as with E the mean Young's modulus, I the moment of inertia, A the area, ρ the density, f max the highest simulated frequency, and λ min the smallest obtained wavelength for the highest input frequency.
For the considered beam configuration, it has been checked that at least six elements are used to represent the wavelength on the coarsest grid for the highest simulated frequency, which in this case is 400 Hz. The FRF results are presented in Fig. 20. As was the case for the static elastic and elastoplastic response, the shades of blue represent the PDF, with the blue line being the most probable value, and the orange line the average value. As can be observed, the uncertainty bounds for the FRF are wider and more spread out when the Young's modulus is homogeneous, Fig. 20 (left), as opposed to a heterogeneous Young's modulus, Fig. 20 (right). This discrepancy is due to the fact that in case of a homogeneous Young's modulus, the resonance frequency will be shifted for each different sample. Averaging all these samples gives rise to a broad and wide uncertainty bound. In case of a heterogeneous Young's modulus, the different samples compensate each other, in analogy with the explanation given in §4.3.1. This gives rise to much smaller uncertainty bounds. Fig. 21 shows the resulting FRF for ten realizations. An important inequality that must hold for MLMC to work well is It has been observed empirically that (32) is not necessarily fulfilled for resonance frequencies and frequencies close by those resonance frequencies for the dynamic elastic response in case of a heterogeneous Young's modulus. The reason for this phenomenon appears to be closely linked to the resolution of random fields on the different levels, as discussed in §4.1. In order to remedy to this, we use the following strategy. Using a small number of samples, the magnitude of V[P 1 ] − V[P 1 − P 0 ] is estimated. When the estimation is below a certain threshold T , the coarsest level is discarded and the algorithm is started on a finer mesh. We write this condition in a log 2 base, in order not to lose consistency with Fig. 12, which depict the rates in a log 2 base. The estimation after simplification is written as For the experiments reported here, we selected the right-hand side of (33) heuristically to be equal to 2.3.

Conclusion
In this work, it has been demonstrated that the Multilevel Monte Carlo method provides a significant computational cost reduction and speedup compared to the standard Monte Carlo method for computing the stochastics of the response in a structural engineering problem. This has been shown by means of actual computing times and normalized computational costs. Speedups are typically of the order of a factor ten for moderate accuracies and can be even much higher for higher accuracies. In addition, figures depicting the static elastic, the static elastoplastic and the dynamic elastic response have been presented with uncertainty bounds. To model the uncertainty in the Young's modulus, we considered both a homogeneous model, represented by a single random variable, and a heterogeneous model, modeled by means of a random field.
We observed that the nature of the uncertainty in the Young's modulus has a big impact on the uncertainty characteristics of the simulation results. For the homogeneous Young's modulus, the variance in the solution is larger than for the heterogeneous Young's modulus. This is apparent from the wider and more spread out uncertainty bounds in §4.3.1. Hence, for realistic applications, the appropriate choice of the uncertainty model is of great importance.
Further paths of research will focus on ways to exploit the similarities between the responses for neighboring frequencies and to couple the Finite Element solvers with an existing Multilevel/Multi-index framework [9,10].