On the optimal strategy of stochastic-based reliability assessment of railway bridges for high-speed trains

The scope of this paper is to evaluate the performance and computational efficiency of various stochastic simulation methods for a stochastic based reliability assessment of railway bridges subjected to high-speed trains. Depending on the degree of sophistication, application of crude Monte Carlo simulation to a realistic mechanical model of the uncertain bridge-train interacting dynamical system can be prohibitively expensive. Thus, three alternative stochastic methods, i.e. line sampling, subset simulation, and asymptotic sampling, are tested on two example problems. These examples represent two classes of bridges with different dynamic response characteristics. While in the one class of bridges distinctive resonance peaks govern the dynamic peak response, the random response amplification of the second group of bridges is primarily induced by track irregularities. The studies are conducted on a simplified mechanical model, composed of a plain beam representing the bridge and a planar mass-spring-damper system representing the train. This modeling strategy captures the fundamental characteristics of dynamic bridge-train interaction, and thus, facilitates the desired assessment of the stochastic methods with reasonable computational effort. It is shown that both line sampling and subset simulation reduce significantly the computational expense for the first class of bridges, while maintaining the accuracy of the predicted bridge reliability. To ensure accuracy and efficiency, these methods need to be modified when applied to systems where track irregularities dominate the random response. For the latter class of bridges, subset simulation proved to be a suitable method for assessing the reliability of this dynamic interacting system when appropriately modified.


Introduction
High-speed rail traffic operates significantly faster than traditional rail transport. Its development poses new challenges to the designer of structures along these lines. For instance, some bridges are excited to excessive vibrations when crossed by trains with high speed, and the ultimate bridge deck acceleration may be exceeded. As a consequence, the ballast becomes unstable, the rail quality is impaired, and the hazard of train derailment increases. Basically, it can be distinguished between two sources that may lead to excessive bridge vibrations. Excitation to bridge resonance at critical speeds is the result of repetitive axle loading of the crossed bridge, transferred through the train wheels to the structure. In this case, at a critical speed, where the frequency of the moving repetitive axle loads or a multiple of it is close to one natural bridge frequency [14], the dynamic bridge response exhibits a distinct resonance peak. The second substantial source of large bridge vibrations are track irregularities [27]. Track irregularities increase the dynamic response not only at certain critical speeds but in a large range of train speeds. Although track irregularities are limited by regulations [11], it should be considered that the track quality decreases with time, and if insufficiently maintained, track irregularities may become large.
From this short discussion it becomes apparent that bridge and high-speed train represent a highly dynamic interacting system, whose reliability should be investigated in more detail by the design engineer. Uncertainties of the train-bridge interaction problem, which considerably influence the structural response, are dispersing material parameters (in particular damping), track irregularities, and environmental effects [16,23]. For instance, if ballast and subsoil start to freeze due to the prevailing temperature conditions, the global stiffness of the bridge, and consequently, its natural frequencies increase. A summary of the identified uncertainties of ballasted railway bridge-train interaction and their quantification and assumed distribution is provided in Salcher [20] and Salcher et al. [24]. As outlined in several studies, such as Gulvanessian et al. [13] and Ü lker-Kaustell and Karoumi [26], the bridge acceleration response is often the limiting and decisive quantity when a dynamic assessment of this interacting system is necessary, whereas the ultimate load bearing capacity of the bridge is usually not exceeded. In a comprehensive reliability analysis, a stochastic simulation method is applied, based on a sufficiently sophisticated mechanical model of the bridge-train interaction problem, and considering structural, environmental and excitation uncertainties. The computational effort of crude Monte Carlo simulation (MCS), which is robust, reliable, and conceptually the most simple stochastic method to predict failure probabilities of highly non-linear problems, becomes rapidly prohibitively large when evaluating the nonlinear limit state function of the considered problem for the required small failure probabilities. Therefore, it is desirable to use more efficient stochastic methods to deal with high-dimensional reliability problems including many uncertainties, such as line sampling (LS) [15,19], subset simulation (SS) [1], and asymptotic sampling (AS) [4], among others. In comparison to crude MCS, depending on the considered problem, those stochastic methods reduce the variance of the probability and simultaneously also the number of computationally expensive limit state function evaluations. There are, however, problems where these methods do not necessarily predict the actual probability of failure efficiently or even fail to predict the failure probability. It is, thus, the objective of the present contribution to evaluate suitability and efficiency of LS, SS, and AS for assessing the reliability of the bridge-train interaction problem.
In the simplest possible modeling approach, where dynamic interaction between bridge and train is explicitly considered, an uncertain beam represents the bridge structure, and the train is described by an uncertain planar mass-spring-damper (MSD) system. In this approach, which is employed in this study, coupling between the beam and the MSD system is achieved through the so-called corresponding assumption [29]. This modeling strategy captures all essential characteristics of uncertain dynamic bridgetrain interaction with lowest possible numerical effort, and thus, the problems can be solved for the required number of samples. A brief summary of different bridge-train models of various sophistication is, for instance, provided in Cantero et al. [6].
In the present study, two bridge models subjected to the Austrian high-speed train type Railjet serve as example problems. The dynamic acceleration response of one bridge is governed by a distinctive resonance peak induced by the repetitive loads of the Railjet train. In the second example problem, where in the considered speed range no pronounced resonance peak is observed, the dynamic response amplification due to track irregularities becomes critical. These example problems represent two classes of the bridgetrain interaction systems with fundamentally different dynamic response characteristics that lead to an exceedance of the ultimate bridge deck acceleration.
The paper is organized as follows. At first, a short description of the utilized bridge-train interaction model, and its uncertainties and random variables is provided. For the two example problems, reliability analyses are conducted utilizing crude MCS. The outcomes of these computationally expensive analyses serve as reference solution for the subsequent evaluation of LS, SS, and AS applied to predict the reliability of the example problems. The paper is concluded by a critical discussion about advantages and drawbacks of LS, SS, and AS when applied to dynamic train-bridge interaction systems.
2 Modeling approach

Mechanical bridge-train interaction model
Based on the substructure approach, the equations of motion are specified for both the bridge and the train subsystem model separately, and then coupled by imposing coupling conditions at the interface of the subsystems.
In this study, the bridge substructure is modeled as a simply supported continuous Euler-Bernoulli beam depicted in Fig. 1. Beam deflection w(x, t) of the undamped Euler-Bernoulli beam bridge with uniform bending stiffness EI and mass per unit length qA subjected to a series of N w concentrated forces F ðiÞ b ðtÞ, i ¼ 1; . . .; N w , moving with constant speed v is governed by the following partial differential equation of motion [8,28], qA € wðx; tÞ þ EIw ;xxxx ðx; tÞ Position n i ¼ vt À s i of the ith interaction force F ðiÞ b ðtÞ between track and the ith wheel pair of the train at time t is indicated by delta function dðx À n i Þ, see Fig. 1 U n ðxÞq n ðtÞ, decouples the partial differential equation of motion into N U ordinary single degree-of-freedom oscillator equations in terms of modal coordinates q n ðtÞ [30], with In Eq. 2 bridge damping f n has been added modally. The N U modal equations of motion written in matrix form, represent the mechanical model of the bridge subsystem. Herein, M b , C b , K b , p b denote the modal mass, damping and stiffness matrix, and modal interaction force vector, respectively, In the present study, the dynamic response of the considered bridge structures is approximated by seven modes (i.e. N U ¼ 7), verified in convergence studies. The N c vehicles of the train subsystem are modeled as planar MSD systems consisting of rigid bodies with mass, which represent passenger stage, two bogies and and four wheel pairs, connected by spring-dashpot elements. As such, the interaction between the bridge and the vehicle can be considered, and consequently, also the effect of track irregularities. Figure 2 shows the MSD system, representing both the power car and the passenger cars of the Austrian high-speed train type Railjet used in this study as reference high-speed  [20] train. In total, this vehicle model exhibits ten degreesof-freedom (DOFs), i.e. seven translational DOFs, The equations of motion of the jth vehicle are, for instance, derived by Lagrangian equations [30]. In a common assumption that horizontal interaction between the vehicles can be neglected, the equations of motion of the train subsystem are expressed as [20] where mass matrix M c , damping matrix C c , and stiffness matrix K c are composed of the corresponding vehicle submatrices Vehicle matrices M ðjÞ c , C ðjÞ c , and K ðjÞ c are found in Salcher [20]. Displacement vector u The Railjet train model used in this study is composed of one power car and seven passenger cars, i.e. in total N c ¼ 8 vehicles and N w ¼ 32 wheel pairs. The equations of motion of the beam and MSD subsystems, Eqs. (4) and (6), respectively, are coupled according to the so-called corresponding assumption [21]. This assumption implies that the elasticity of wheel and track at the contact point is neglected and a continuous contact of the bodies is supposed. Consequently, the vertical displacement of the ith wheel pair, u The vertical track displacement is composed of bridge deflection wðn i ðtÞ; tÞ and track irregularity described by irregularity function I r ðn i ðtÞÞ. From the corresponding assumption it also follows that the contact force between each wheel pair and track is equal,   [21] the applied CMS method, the involved system matrices and vectors are described in full detail.

Uncertainties and random variables
The uncertainties of a bridge-train interaction system can be categorized into three groups. The first group is related to the bridge, the second group to the excitation  [21] induced by the passing train, and the third group to the environmental impact [22]. In this study, geometry and material parameters of the train are assumed to be deterministic.
According to Salcher [20] and Salcher et al. [24] the most important uncertain bridge parameters are damping, dispersing structural parameters, and the properties of the ballast. In the current model, all dissipation mechanisms of the bridge including structure, subsoil and track are captured globally by modal damping coefficients f n (n ¼ 1; . . .; N U Þ, see Eq. 2, assumed to be the same for all modes, To describe the dispersion of damping, coefficient f is modeled as truncated lognormally distributed random variable. Its parameters listed in Table 1 (i.e. mean f ¼ 1% and coefficient of variation CV ¼ 0:3% [24]) refer to a ballasted steel bridge [24]. The truncation at f ¼ 0:5% excludes unlikely low damping values [24]. The Metropolis algorithm [18] is used to generate the truncated lognormally distributed samples of f. To consider the effect of naturally dispersing structural properties, moment of inertia I of the bridge girder is treated as Gaussian distributed random variable. Mean and CV are listed in Table 1. Due to water saturation and sediment pollution, the ballast is a bridge element with strongly dispersing properties. To capture this behavior, in the model ballast density q b is assumed to be a random variable with uniform distribution, see also Table 1.
At high train speeds, track irregularities represent a substantial source of excitation. In the present bridge model track irregularities are considered via random irregular track profile functions I r ðxÞ, superposed to the beam deflection, compare with Eq. (9). Random profiles I r ðxÞ are assumed to comply with a stationary Gaussian stochastic process [25], realized as a stochastic superposition of J harmonic functions with discrete frequencies X m , m ¼ 1; . . .; J, and uniformly distributed random phase angle u m in the range 0 to 2p [7], In the frequency increment DX ¼ ðX u À X l Þ=J, X u denotes the upper boundary frequency and X l the lower boundary frequency. Note that a ¼ 4 for m ¼ 1, Variable SðX m Þ represents the power spectral density, expressed as one-sided density function U v ðX m Þ [7], with X c ¼ 0:8246 rad/m and X r ¼ 0:0206 rad/m [7]. The uppermost frequency used to generate irregular profiles is X u ¼ 2:1 rad/m, and the lowermost frequency X l ¼ 0:07 rad/m. Amplitude Q, which defines the track quality, is considered as continuous uniformly distributed random variable between the boundaries 0:592 Á 10 À6 radm and 1:586 Á 10 À6 radm. The lower boundary represents a high quality rail and the upper boundary a low quality rail [7], see also  Changing environmental conditions affect the dynamic behavior of the bridge and in further consequence its response. Measurements on ballasted railway bridges have shown that their natural frequencies change almost stepwise close to the freezing point of the surrounding environmental temperature [12]. This behavior can be attributed to the change of water inside the ballast and subsoil to ice at freezing temperature, leading to an increase of the global bridge stiffness. Considering that the natural bridge frequencies are directly affected by the frost depth, in Salcher et al. [24] a stochastic model capturing the temperature dependence of the natural frequencies was proposed, which is used for the current study. In this model, the daily mean air temperature T is considered as a random variable following an extreme value distribution with its realizations T j . The unfrozen state of the ballast assumed to be attained at a daily minimum temperature at ground T g [ À 1 C and the fully frozen state at T g \ À 10 C are defined by random variables T 0 and T 1 , respectively, expressed in terms of the conditional probability distributions CÞ as a function of the daily mean temperature T. If the daily mean temperature T j is inbetween the temperature for the unfrozen and the fully frozen state, Young's modulus of the ballast is determined by linear interpolation of the Young's modulus of the unfrozen and the fully frozen ballast [24]. The parameters of the random variables of this model listed in Table 1 are based on the temperature data of Munich Airport [24].
The number of random variables, N r , without track irregularities is N r ¼ 7. When considering track irregularities, additional 1001 random variables are considered, through the superposition of J ¼ 1000 harmonic functions with discrete frequencies X m and uniformly distributed random phase angle u m . Together with random irregularity amplitude Q, N r ¼ 1008 in total.
2.3 Definition of failure, limit state function, and probability of failure Failure of a structure indicates an undesirable state of the structure, i.e. the structural response exceeds certain limits that can be expressed in terms of a limit state function (performance function) gðXÞ, where X ¼ ðX 1 ; X 2 ; . . .; X N r Þ is the vector comprising the N r random input variables. In general, the load bearing capacity of railway bridge structures is satisfied, however, the admissible vertical bridge deck acceleration may be exceeded when crossed by a high-speed train with critical speed. Large bridge deck accelerations may lead to ballast instability, and thus, it is commonly the limiting factor in the dynamic design of railway bridges [6]. In the present study, bridge deck acceleration € w at a given speed v governs the limit state in the stochastic reliability approach of the bridge-train interaction system. The corresponding safety distance reads as where c bt is the acceleration threshold for a full stochastic reliability assessment of a ballasted railway bridge. In the present study, c bt ¼ 7 m/s 2 is a deterministic quantity that is two times the acceleration limit value of 3.5 m/s 2 associated with ballast instability, recommended in design guideline Eurocode 1 [10] for a semi-probabilistic bridge assessment. In this way, the safety factor of two is eliminated [24]. Failure is predicted if the limit state function gðZÞ 0, with safety distance Z according to Eq. (13).
Probability of failure p f is defined as probability P that an undesired state is reached: p f ¼ PðfX : gðXÞ 0gÞ [5]. Since instability of the ballast is a serviceability problem, according to Eurocode 0 [9] the maximum probability of failure of this class 2 bridge-train interacting system is p f ¼ 10 À3 for the serviceability limit state (SLS). Crude Monte Carlo simulation (MCS) is conceptually the most simple and a robust numerical method to evaluate p f . When applying MCS, N MC independent and identically distributed samples x ðrÞ (r ¼ 1; . . .; N MC ) with distribution as defined before are generated. For each sample, the limit state function is computed (commonly numerically, for instance, by response history analysis). When the response of the sample is inside the failure domain X f ðxÞ ¼ fx : gðxÞ 0g, indicator function I f ðxÞ is one, and zero otherwise. The Þ represents an unbiased estimate of p f [5]. If N MC is sufficiently large, the predicted probability of failure p MC f is reliable, and the efficiency of this method does not depend on the dimension N r of the problem. However, for complex structural systems, such as elaborate three-dimensional mechanical models of bridge-train interaction including the subsoil, the computational costs of MCS are prohibitively large. Thus, subsequently, three alternative stochastic methods, i.e. LS, SS, and AS, are tested on two example problems. A short summary of these alternative stochastic methods as applied in this study can be found in the ''Appendix''.
In the present study, the outcomes of a crude MCS with N MC ¼ 10;000 samples for discrete train speeds in the range from 20 to 100 m/s (speed increment 1 m/ s) serve as reference solutions. With this sample size, the required failure probability p f ¼ 10 À3 can be reliably predicted [24].

Example problems
For the assessment of the considered stochastic methods two different bridge structures are considered. The bridges located close to Munich airport are assumed to be crossed by the Austrian high-speed train type Railjet with the common configuration of one power car and seven passenger cars. The deterministic parameters of the corresponding planar MSD vehicle model according to Fig. 2 are found in Salcher [20].
Example bridge 1 is a simply supported ballasted steel bridge of span L ¼ 16:8 m with mean bending stiffness EI ¼ 3:262 Á 10 10 Nm 2 and mean mass per unit length qA ¼ 1:220 Á 10 4 kg/m. Natural frequencies of this model with assigned mean parameters are f 1 ¼ 9:10 Hz, Figure 3 shows the peak acceleration, max € wðx; tÞ j j, of the mean parameter bridge with perfect (smooth) tracks (red line) and with one arbitrarily chosen sample of irregular rail profiles (black line), respectively, as a function of train speed v in the range 20 v 100 m/s. As observed, the spectra exhibit a distinct local maximum close to critical speed v ð1Þ 2 ¼ 75:1 m/s, i.e. the acceleration response is dominated by one significant resonance peak where max € wðx; tÞ j j! c bt . Note that at a critical train speed, defined as v ðnÞ l ¼ df n =l; l ¼ 1; 2; 3; . . . [14], the nth natural bridge frequency f n ¼ x n =ð2pÞ is excited to a state of resonance due to the repetitive axle loads with constant distance d. For the Railjet train this distance corresponds to wagon length d ¼ 16:5 m. It is also seen that irregular rails amplify considerably the peak acceleration, and close to the resonance peak the acceleration threshold c bt ¼ 7 m/s 2 is exceeded. This model is, thus, representative for bridges whose peak acceleration amplification is governed by resonance due to repetitive axle loads.
For the MCS reference solution, 10, 000 random bridge samples are generated. Figure 4a shows max € wðx; tÞ j j of those samples disregarding track irregularities plotted against speed v. Additionally, mean, 10% and 90% quantiles, minimum and maximum of these spectra are depicted. In the range around critical speed v ð1Þ 2 , acceleration threshold c bt ¼ 7 m/s 2 is frequently exceeded. The corresponding spectra of the bridge considering track irregularities modeled according to Eq. 11 are depicted in Fig. 4b. Comparing Fig. 4b with Fig. 4a confirms that track irregularities increase the peak bridge acceleration, in this example particularly at speeds larger than 60 m/s. The second bridge is a simply supported steel bridge of length L ¼ 21:36 m, mean bending stiffness EI ¼ 2:562 Á 10 10 Nm 2 and mass per unit length qA ¼ 0:813 Á 10 4 kg/m. The first seven frequencies of this model with assigned mean parameters are f 1 ¼ 6:11 Hz, f 2 ¼ 24:4 Hz, f 3 ¼ 55:0 Hz, f 4 ¼ 98 Hz, f 5 ¼ 153 Hz, f 6 ¼ 220 Hz, and f 7 ¼ 300 Hz. Figure 5 shows the peak acceleration spectrum of a bridge model with perfect tracks (red line) and of a model with an arbitrary irregular track profile (black line). In contrast to bridge 1, the 13 ¼ 31:0 m/s are observed. The maximum peak acceleration is much less than acceleration threshold c bt ¼ 7 m/s 2 . However, the results indicate that rail irregularities increase the bridge response significantly, and hence, a reliability assessment is conducted considering the random variables introduced previously.
In Fig. 6a the beam peak acceleration of N MC ¼ 10; 000 random bridge samples for train speeds in the range from 20 to 100 m/s is depicted. These results show that at bridge 2 with perfectly smooth rails acceleration threshold c bt ¼ 7 m/s 2 is exceeded so infrequent that the probability of failure is less than 10 À3 . Spectra of bridge model 2 considering track irregularities are depicted in Fig. 6b. Especially at train speeds larger than 50 m/s track irregularities increase the peak accelerations compared to the samples with perfectly smooth rails. The exceedance of c bt is primarily caused by track irregularities, i.e. the peak bridge accelerations of bridge structure 2 is significantly dominated by track irregularities. Example bridge 1 Figure 7a shows probability of failure p f of the bridge with perfect (smooth) rails in logarithmic scale, plotted against train speed v in the critical range from 60 m/s up to 90 m/s. Circular markers refer to the reference solution derived by crude MCS with 10, 000 samples. As observed, at speed v ¼ 71 m/s, p f estimated by crude MCS is 10 À4 . At lower speeds, p f is less than 10 À6 , and thus, cannot be predicted with the used sample size. With increasing v also p f increases, and for a probability of failure p f ¼ 10 À3 the acceleration threshold c bt ¼ 7 m/s 2 is the first time exceeded at v ¼ 72:0 m/s. Note that p f ¼ 10 À3 corresponds to ten failed samples out of 10, 000. The  Additionally, in this figure also the outcomes of LS based on N LS ¼ 100 lines (star-shaped markers), SS (plus-shaped markers), and AS (square-shaped markers) are depicted. Line Sampling as utilized in this example is based on simplified limit state evaluation of the root by quadratic interpolation, as explained Appendix ''Line sampling''. For SS of this particular study, an initial sample size of N SS ¼ 100 and intermediate failure probabilities of p interm ¼ 0:1 have been selected. The limit state function is evaluated for the N SS initial Monte Carlo samples. The number of samples that lie closest to the failure domain of these N SS initial samples is N SSw ¼ N SS p interm . In this example, at each of the following ðk À 1Þ substeps, the worst 10% of the samples serve as seeds for Markov chains of length C l ¼ 10 used to determine the intermediate probabilities of failure according to the conditional probability density functions qðxjF m Þ, see Eq. (18). Asymptotic sampling as applied in this contribution with N AS ¼ 100 initial samples, starts at (1=r ¼ 1), evaluating the limit state function N AS ¼ 100 times. Subsequently 1=r is successively reduced by 0.1 until failure of the structure is reached, leading to one pair of ð1=r; bÞ. Five pairs of ð1=r; bÞ are used to determine bð1Þ by a least squares fit.
It is seen that LS, SS, and AS with the selected resolution predict the general trend of p f with respect to speed v, however, at some specific speeds discrepancy and scatter compared to the crude MCS solution is quite large. As discussed in Appendix ''Line sampling'', importance direction e a of LS points towards the gradient of g(Z) at the origin in standard normal space. According to the outcomes of LS with N LS ¼ 100 lines, the desired p f of 10 À3 is exceeded at a speed of v ¼ 74 m/s. Subset simulation and AS predict for p f ¼ 10 À3 a maximum admissible speed of v ¼ 73 m/s and v ¼ 72 m/s, respectively. Figure 7b shows p f with respect to v for the bridge with track irregularities, derived by the same stochastic methods with same parameters as before. The p fpredictions are based on simplified LS, simplified SS, and simplified AS. In particular, track irregularities are neglected in the algorithms of the stochastic procedures and only considered when evaluating the limit state function. That is, the importance direction of simplified LS is determined without considering track irregularities. Concerning simplified SS, the Markov chains are evaluated leaving the random phase angles of the irregularity profiles out of consideration. The number of random variables for simplified AS is the same as for the bridge with perfect rails, i.e. N r ¼ 7.
For instance, at v ¼ 72 m/s, for the structure with irregular tracks p f estimated by crude MCS is 1:4 Á 10 À3 , while for the perfectly smooth bridge p f ¼ 1:1 Á 10 À3 . However, the maximum probability of failure is similar for both structural models, i.e. p f ¼ 0:209 at v ¼ 78 m/s (irregular rails), and p f ¼ 0:192 at v ¼ 78 m/s (perfect rails). Interestingly enough, also in the bridge with irregular profiles the target probability of failure of 10 À3 is exceeded at the same speed as in the object with perfect tracks, i.e. v ¼ 72 m/s. Also here, the difference of p f -predictions from LS, SS, and AS compared to the reference MCS outcomes are significant. In particular, at several discrete speeds SS does not predict failure with p f [ 10 À6 with the employed number of initial samples.
In a subsequent study, LS with N LS ¼ 300 initial lines, SS with N SS ¼ 300 initial samples (p interm ¼ 0:1) and AS with N AS ¼ 300 initial samples is conducted, in an effort to determine p f more reliably. Figure 8a for perfect tracks and Fig. 8b for irregular track profiles show that an increased number of lines respectively increased initial sample size significantly improves the accuracy of these methods. In the entire considered speed range, the difference between reference p f from crude MCS, and the outcomes of LS, SS and AS is small. It can be, thus, concluded that in this example problem simplified LS, simplified SS and simplified AS predict accurately the probability of the failure if the number of lines and initial samples, respectively, is sufficiently large. This outcome also confirms the chosen strategy, where the importance direction of simplified LS and the Markov chains of simplified SS have been determined without considering the random variables of the track irregularities. In the simplified AS approach, the random variables of the irregular profiles have been also neglected when initially generating the random samples. Such simplified LS, SS and AS procedures perform well, because in the considered bridge-train interaction problem failure is governed by resonance close to critical speed v ð1Þ 2 ¼ 75:1 m/s. Line sampling allows also to determine failure domain X f ¼ fZ : gðZÞ 0g [17] of the bridge-train interaction system. To this end, limit state function g(Z) is evaluated at certain distances b along the random lines in standard normal space (see Appendix ''Line sampling''). As outcome, Fig. 9 shows for train speed v ¼ 80 m/s the limit state function of the bridgetrain interaction system (a) without track irregularities and (b) considering track irregularities along N LS ¼ 100 lines plotted in standard normal space against LS parameter b, with À 5 b 5. Direction a of the walk along the lines corresponds to importance direction e a of LS with N LS ¼ 100 lines. If gðZÞ 0, the failure state of the bridge-train interaction system is reached, i.e. the actual peak bridge acceleration, max j € wðx; tÞj, exceeds the threshold of the permissible deck acceleration c bt ¼ 7 m/s 2 , i.e. c bt À max j € wðx; tÞj 0. A positive limit state function implies that the bridgetrain interaction system is safe. As can be seen, the limit state function is negative only in a certain range of b, i.e. in the range around critical speed v ð1Þ 2 . At lower and higher speeds the acceleration threshold of 7.0 m/s 2 is not exceeded, see also Figs. 4a, b. This result shows that the failure domain of the considered bridge-train problem is bounded. It is also observed that the scatter of g with respect to b is much larger for the bridge samples with irregular rail profiles than for objects with smooth track surface, compare with Fig. 9b. Reason behind this response behavior is that in the importance direction of e a for LS track irregularities are not taken into account but only the limit state function is assumed to be random.
Crude MCS with increased standard deviation does not necessarily generate more samples in a bounded failure domain as in the considered problem. Consequently, AS with N AS ¼ 100 initial samples, which has been developed for unbounded failure domains, does not accurately predict p f of the bridge. In theory, LS is also based on the assumption of an unbounded failure domain, however, in the present study LS has been proven to be more robust than AS when using N LS ¼ 300 lines.

Example bridge 2
Subsequently, the assessment of probability of failure p f of example bridge 2 by the considered stochastic methods is discussed. Figure 10 depicts p f in the critical range from 60 m/s up to 90 m/s resulting from crude MCS (circular markers), LS with 100 lines (starshaped markers), SS with 100 initial samples (crossshaped markers), and AS with 100 initial samples (square-shaped markers) for bridge 2 with perfect (left subplot) and irregular tracks (right subplot).
In the perfect bridge, according to the outcomes of a crude MCS with sample size N MC ¼ 10; 000, solely at speed v ¼ 83 m/s the object fails with the failure probability p f ¼ 2 Á 10 À4 , compare also with Fig. 6a. Since crude MCS with N MC ¼ 10; 000 samples estimates only failure probabilities up to 10 À3 with a coefficient of variation of 0.3 reliably, another computation of p f at speed v ¼ 83 m/s with the appropriate number of N MC ¼ 100;000 samples has been conducted. The latter MCS yields a failure probability of p f ¼ 1:5 Á 10 À4 , which is close to the outcome of the first MCS (i.e. p f ¼ 2 Á 10 À4 ) with N MC ¼ 10; 000 samples. At discrete speeds other than 83 m/s no failure is predicted larger than 10 À6 . Line sampling, SS, and AS with an equivalent setup as in bridge example 1 and the resolution specified above fail to  In a further effort to predict p f more precisely, the resolution of LS and SS is increased, as in the bridge example 1. That is, LS with N LS ¼ 300 initial lines, and SS with N SS ¼ 300 initial samples and p interm ¼ 0:1 is conducted. For perfect tracks, according to the results shown in Fig. 11a, at speed v ¼ 83 m/s LS and SS with enhanced resolution deliver a probability of failure very close to the MCS reference solution. In particular, LS with N LS ¼ 300 lines results in p f ¼ 1:47 Á 10 À4 , and SS with N SS ¼ 300 initial lines yields p f ¼ 2:6 Á 10 À4 , compared to p f ¼ 2 Á 10 À4 from MCS with N MC ¼ 10; 000 initial samples. In contrast to crude MCS, LS and SS also predict failure probabilities at speeds close to v ¼ 83 m/s, however, they are too small to be predicted by MCS with the considered sample size. Figure 11b shows that a larger sample size does not enhance the capability of simplified LS, simplified SS, and simplified AS to predict reliably p f of bridge 2 in the presence of rail irregularities. The outcomes of LS, SS, and AS diverge strongly from the reference crude MCS results. It can be concluded that LS, where the importance direction e a has been estimated without considering random irregularity amplitudes Q and random phase angle u m , cannot predict p f if the track irregularities are the governing source of failure. Likewise, for this type of railway bridges, SS fails to predict p f because the acceptance of new samples are disturbed if Q and u m are not considered in the Markov Chain Monte Carlo procedure used to generate the conditional samples. This is in contrast to bridge example 1 with one dominating resonance speed, where the simplified versions of LS, SS, and AS (i.e. random irregularity variables are considered only in the evaluation of the limit state function) allow to predict p f of bridges with random irregularity profiles. Subsequently, modifications of SS and AS are proposed aiming at improving their performance for a bridge-train interaction problem whose failure is governed by track irregularities. Subset simulation is modified in this respect that random amplitude Q and random phase angle u m are taken into account when determining the acceptance ratio of the Markov Chain Monte Carlo algorithm, which is necessary to generate the conditional samples. This modification does not impair the efficiency of SS because no additional evaluations of the limit state function need to be conducted. To account in AS explicitly for track irregularities, initial samples N AS of all N r ¼ 1008 random variables (i.e. also including the random variables of the irregularity profiles) are generated. This modification affects only slightly the efficiency of AS because again no additional limit state evaluations are necessary. In AS as applied in bridge problem 1, an initial vector is generated using the seven random variables of the perfect bridge, and the rail irregularities are only considered when evaluating the limit state function. The outcomes show that these modified versions of SS and AS yield predictions of p f that are still far-off the reference solution, and are, thus, not reliable.
However, repetitive prediction of p f using a small number of initial lines/samples and averaging over these multiple p f -predictions improves significantly the accuracy of LS (quadratic interpolation) and modified SS. Figure 12 depicts the mean of p f determined by LS applied five times with N LS ¼ 100 initial lines each. As observed, the results of LS based on such procedure reflect the global trend of p f of bridge-train interaction problem 2, however, at certain speeds the differences between the LS predictions and the outcomes of MCS are still large. In contrast, the average of p f as outcome of five modified SS applications with N SS ¼ 100 initial samples each also shown in Fig. 12 reveal the accuracy of this modified SS approach.

Discussion on the computational efficiency
In LS, the central difference quotient has been used to determine the gradient at the origin in standard normal space, which corresponds to importance direction e a , based on the N r ¼ 7 independent random variables in the considered example (without the random variables of the track irregularities). Thus, additionally to the computations of the Monte Carlo samples as starting points for each line (i.e. in this study N LS ¼ 100 and N LS ¼ 300, respectively), the limit state function needs to be evaluated 2N r times more, with N r denoting the total number of random variables. In a first LS approach, where the root gðZÞ ¼ 0 is found by quadratic interpolation, for each of the N LS lines the limit state function is evaluated at three points, requiring altogether additional 3N LS limit state function evaluations. Hence, for this LS approach as performed for the assessment of example bridge 1, the limit state function is evaluated in total N LStot ¼ 3N LS þ 2N r times to predict probability of failure p f for a single train speed v. That is, with N r ¼ 7 independent random variables in the considered example and N LS ¼ 100 initial lines, the total number of simulations N LStot adds up 314. When taking 300 initial lines, N LStot is 914, which is significantly less than the 10, 000 samples of crude MCS for a p f ¼ 10 À3 . Consequently, considering p f at a distinct speed v, LS reduces reasonably the number of limit state function evaluations compared to crude MCS. In the second example bridge, the first LS approach repeatedly applied five times with N LS ¼ 100 initial Fig. 12 Probability of failure p f based on crude Monte Carlo simulation, line sampling (5 Â 100 lines) and modified subset simulation (5 Â 100 initial samples). Example bridge 2 with irregular track lines each, is computationally more efficient than crude MCS. However, at certain train speeds p f obtained by averaging the individual probability of failures diverges from the crude MCS outcome because the track irregularities impair the simplified estimation of e a .
For SS the total number of limit state function evaluations, referred to as N SStot , depends on the actual probability of failure, the number of initial seeds, and the length of the Markov chains used to generate conditional intermediate samples. Since in the present study for the intermediate probability of failure a value of p interm ¼ 0:1 has been specified, 10% of N SS serve as Markov chain seeds N SSw . In the next ðk À 1Þ substeps, the limit state function is evaluated ðN SS À N SSw Þ times. The iteration is stopped when the actual intermediate probability of failure at the current step is larger than the predefined value of p interm ¼ 0:1. Thus, in total, subset simulation comprises N SS evaluations of the limit state function at the initial MCS step, and ðN SS À N SSw Þ evaluations at each of the ðk À 1Þ substeps. To predict a p f of 10 À1 , the initial sample of N SS ¼ 100 without any substeps is sufficiently large. A p f of 10 À3 requires N SS ¼ 100 initial samples plus ðN SS À N SSw Þ samples at ðk À 1Þ ¼ 2 substeps, i.e. in total N SStot ¼ 280 limit state function evaluations. When selecting an initial sample size of N SS ¼ 300, the limit state function needs to be evaluated N SStot ¼ 840 times to predict a p f of 10 À3 reliably. Repetitive application of modified SS based on N SS ¼ 5x100 initial samples is more efficient than crude MCS to determine p f of the second bridge-train interaction problems where track irregularities dominate the structural response amplification. This results in total in N SStot ¼ 5x280 ¼ 1400 limit state function evaluations to determine p f in the range of 10 À3 .
The computational efficiency of AS depends on p f . For the present bridge-train interaction problem 1 considering failure probabilities above 10 À3 in the speed range of 72 to 84 m/s failure is usually reached at (1=r ¼ 1) or (1=r ¼ 0:9) or (1=r ¼ 0:8). Since five pairs are used to determine p f , the limit state function is evaluated in total at least N AStot ¼ 500 times. The simulation time for one limit state function evaluation was on an average about 5 s but depends on the actual train speed. In the actual study, due to the large computational effort of crude MCS with N MC ¼ 10; 000 samples, the simulations were carried out on two different machines with different architecture. Crude MCS with N MC ¼ 10; 000 samples was conducted on the high performance computer cluster LEO3 of the University of Innsbruck, split into 100 packages yielding 100 response spectra each. The computation time for one package was t MC % 11:3 h, whereas sequential computing would have lasted 47 days. At the cluster, the packages were solved in parallel. The simulations for LS, SS and AS were conducted on a local machine, i.e. an Apple workstation with 6 core Intel(R) Xeon(R) CPU E5-1650 V2 @ 3.50GHz processor. On an average for one discrete speed, LS with N LS ¼ 100 lines took t LS 100 % 0:44 h and LS with N LS ¼ 300 lines t LS 300 % 1:3 h. The prediction of the 31 discrete failure probabilities in the critical speed range from v ¼ 60 m/s up to v ¼ 90 m/s added up on an average to t totLS 100 % 13:5 h (N LS ¼ 100) and to t totLS 300 % 40:3 h (N LS ¼ 300 lines), respectively. Subset simulation with N SS ¼ 100 initial samples took for one discrete train speed about t SS 100 % 0:38 h, and in total t totSS 100 % 12:1 h for the critical speed range. For the larger initial sample size of N SS ¼ 300, the computation time increased to t SS 300 % 1:2 h (one speed) and t totSS 300 % 37:2 h (critical speed range), respectively. Asymptotic sampling with an initial sample number N AS ¼ 100 yielded a minimum computation time of t AS 100 % 0:14 h, and with N AS ¼ 300 t AS 300 % 0:42 h for one discrete train speed. The prediction of failure probabilities in the order of p f 10 À3 with AS led rapidly to a significant increase of the computational effort.

Conclusions
In this paper, applicability and efficiency of line sampling (LS), subset simulation (SS), and asymptotic sampling (AS) for a stochastic-based reliability assessment of railway bridges subjected to high-speed trains were evaluated. The limit state of this problem is governed by the bridge deck acceleration. These stochastic methods were employed to predict failure probabilities of two bridge-train interaction problems with different response characteristics. In example problem 1, structural accelerations leading to the exceedance of the limit state are related to one significant resonance peak, whereas in example problem 2 track irregularities are the primary source of inadmissible large peak bridge acceleration amplifications. Random parameters were defined to capture the essential uncertainties.
Compared to crude Monte Carlo simulation (MCS), in example bridge 1 LS, SS, and AS reduce significantly the number of limit state function evaluations for a reliable prediction of the probability of failure p f in the order of 10 À3 . It was shown that LS with a resolution of 300 random lines pointing towards the importance direction delivers a failure probability close to the reference solution obtained by crude MCS based on 10, 000 samples. For bridge 1, performance and accuracy of SS turned out to be sufficient when using 300 initial samples. It was also shown that for this object accuracy and efficiency of AS is similar to LS and SS.
The dimension of the systems with perfectly smooth tracks is low, i.e. in the current study the number of random variables is seven. When taking track irregularities into account, the problem becomes of high dimension with 1008 uncertain parameters involved. This high-dimensionality puts some constraints on the efficient applicability of LS, SS, and AS to dynamic bridge-train interaction. Thus, in a simplified approach, the importance direction in LS and the acceptance ratios of the Metropolis-Hastings algorithm in SS are purely based on the limited number of random variables of the system with smooth tracks, whereas the limit state function is evaluated considering the full set of random variables including the ones of the random irregularity profiles. For example bridge 1, whose peak accelerations are dominated by resonance at a critical speed, this simplification did not impose any limitation to predict the probability of failure reliably by LS, SS, and AS. However, in the second example object 2, where track irregularities govern the acceleration amplification, such simplified approach impairs LS, SS, and AS, and consequently, these methods fail to predict the actual failure probability.
Modification of SS by considering the random variables of the track irregularities in the Markov Chain Monte Carlo procedure improves the accuracy of p f -estimations also for example bridge 2. The robustness can be further enhanced by a serial application of LS and SS with a relatively low number of lines and initial samples, respectively, and then, averaging the individual p f estimates. Conducting five times LS with 100 initial lines on bridge example 2 allowed to predict qualitatively p f with respect to speed v. Modified SS based on 100 initial samples and five reiterations yields also quantitatively the reference p f obtained by crude MCS.
From this study it can be concluded that LS, SS, and AS are reliable methods to predict the probability of failure of the bridge-train interaction problem when appropriately applied. The required sophistication of these methods, and thus, their computational efficiency, depends on the response characteristics of the considered bridge. The failure probability of railway bridge objects, where mainly track irregularities amplify the random response, can only be predicted by more robust and less computationally efficient stochastic methods, such as repeated application of modified SS and crude MCS. dimensional systems with many uncertain parameters, where the evaluation of limit state function gðxÞ is numerically expensive. Line sampling utilizes random lines instead of random points to examine the failure domain. Estimatep LS f of probability of failure p f is obtained by computing the one-dimensional failure probability along these lines pointing towards failure domain X f ðxÞ of the problem. The direction of the random lines, referred to as importance direction a, points towards X f ðxÞ. To put it another way, a is defined as the direction of greatest impact on performance function gðxÞ in standard normal space [19]. In the present study, the gradient of the performance function, serves as importance direction a because it points towards the greatest rate of increase of performance function gðxÞ. Determining single points in the random variable space leads to computationally demanding evaluations of the nonlinear limit state function for each point of interest. Calculation of gradient rgðxÞ at point x of the implicitly available limit state function gðxÞ is achieved through 2N r limit state function evaluations (such as the central difference quotient) of gðxÞ [19], compare with Eq. (14).
Then, the normalized direction e a ¼ À rgðxÞ krgðxÞk can be computed, which provides information on the parameter variation leading to the greatest impact on gðxÞ.
To estimate p f by LS, first direct MCS is carried out by generating N LS samples x ðrÞ , r ¼ 1. . .N LS . Then, for each of the N LS samples the one-dimensional reliability problem (along one line) is solved, where the probability of failure of the rth one-dimensional reliability problem p ðrÞ f is given by p ðrÞ f ¼ UðÀb ðrÞ Þ [19]. Herein, U denotes the standardized Gaussian cumulative distribution function, and b ðrÞ is the safety index of the rth one-dimensional reliability problem. Safety index b ðrÞ is the distance from the origin to the boundary of the failure domain, g ¼ 0, along one line. An approximation of b ðrÞ can be obtained by quadratic interpolation [19], where the values of the limit state function evaluated at three different b-values serve as data points for a fitted second-order polynomial, whose root approximates b ðrÞ .
After analysis of p ðrÞ f for each of the N LS samples, the LS estimator of the probability of failurep LS f is determined [19], The computational effort to numerically determine the gradient increases with the dimension of the problem and may even become prohibitively large in high dimensions.

Subset simulation
The basic concept of the subset simulation (SS) method [1]  The SS procedure is stopped when the actual intermediate probability of failure at the current step is larger than the predefined value of p interm . According to Eq. (16), k steps are needed to determine the initial probability of failure p f ¼ PðFÞ ¼ PðF k Þ by SS. Hence, the quantity of the required SS steps, k, depends also on the actual probability of failure of the problem. Computing the required conditional probabilities of failure PðF mþ1 jF m Þ, generation of samples based on the conditional distribution of x given by Au and Beck [1] qðxjF m Þ ¼ qðxÞI F m ðxÞ=PðF m Þ ð 18Þ is necessary. Sampling according to qðxjF m Þ generates samples that lie in F m . This non-trivial task can be solved by Markov Chain Monte Carlo (MCMC) simulation, based in the present study on a modified Metropolis-Hastings (MMH) algorithm [2] that overcomes simulation problems in high dimensional space. Length C l of the Markov chains is set to C l ¼ 1=p interm , and for the proposal density p Ã l a standard normal distribution N ðx s ; IÞ centered at sample point x s is used.

Asymptotic sampling
Asymptotic sampling (AS) [4] is based on the observation that probability of failure p f exhibits in standard normal space asymptotic behavior. Hence, it approaches zero as standard deviation r approaches zero. Safety index b of a structure is defined as b ¼ U À1 ð1 À p f Þ, where U À1 denotes the standardized inverse Gaussian distribution function. Due to the asymptotic behavior of p f , for linear limit state functions safety index bðrÞ for r [ 1 and bð1Þ are related as [4] bðrÞ ¼ bð1Þ r ð19Þ b(1) for r = 1 represents the safety index of the initial reliability problem. By artificially increasing variance r of the random variables, more events in failure domain X f ðxÞ are generated, and safety index bðrÞ can be determined by MCS with affordable computational costs. Extrapolation of bðrÞ by multiplying the result of Eq. (19) by 1=r yields the safety index of the original problem, bð1Þ. Therefore, an estimate of the probability of failure is obtained as [4] p f %p AS f ¼ UðÀbð1ÞÞ ð20Þ The independence from the dimension of the parameter space is one major benefit of this method [4]. The shape of the surface of limit state function gðxÞ, and the relation between the number of samples and probability of failure p f have a distinct impact on the accuracy of AS. As outlined in Bucher [4], AS expects safety index b to decrease with increasing scale factor r. This is the case for the widely spread definition of failure that loads exceed a certain limit, i.e. failure domain X f ðxÞ is unbounded. However, failure domain X f ðxÞ of nonlinear reliability problems may be finite or separated into several domains. Then, increasing of r may not lead to more samples inside the failure domain and AS can not be used to compute p f at given r efficiently. Thus, AS is not always applicable for reliability problems with bounded failure domains efficiently [4].