Model order reduction in subset simulations using the proper orthogonal decomposition

The crude Monte Carlo method is computationally expensive. Hence, incorporating model order reduction methods enabling reliability analysis for high‐dimensional problems is necessary. However, this strategy may result in an inaccurate estimation of the probability of failure for rare events for two reasons. First, the model order reduction, represented by the proper orthogonal decomposition (POD) here, requires response information in the form of snapshots a priori. To capture the essential nonlinear response behavior, we propose to update the proper orthogonal modes using extreme events. Second, the crude Monte Carlo simulation requires many samples to estimate low failure probabilities reliably. To this end, subset simulation found wide application in reliability analysis to reduce computational effort. Following this strategy, the proposed samples gradually move toward the failure region. Thus, incorporating updates of the modes is particularly promising in evaluating samples from the current subset region. This contribution shows the computational efficiency of POD within subset simulations. We then propose to leverage the estimation of the probability of failure by updating the modes within each subset.

methods, machine learning techniques enable rapid estimations.With a special focus on rare event estimation, dataset extensions [6] and selection strategies [7] have been proposed.2. Reducing the number of sample evaluations, that is efficient sampling: Instead of sampling from the whole distribution as performed in the crude Monte Carlo simulation, strategies are employed to achieve reliable failure estimations with a reduced number of samples.Importance sampling [1], stratified sampling [8], and subset simulation [9], for example, are strategies to achieve failure estimations with reasonable computational cost.
These strategies attract increasing attention in structural reliability analysis [10].Addressing both, the incorporation of the proper general decomposition and adaptive importance sampling has shown promising results [11].Furthermore, Kriging was used as a surrogate model and combined with importance sampling for the estimation of small failure probabilities [12,13].Furthermore, the computational advantages of neural networks have been used to accelerate subset simulations [14].In this contribution, first, to reduce the cost of evaluating each sample, we propose to use the POD.This method requires the solution of the system prior to the response evaluation.Therefore, the POD modes have to be evaluated from the full system.Afterward, the reduced system offers significant speed-up capabilities.Thus, POD can be used to solve many time history records based on the result of one sample, for example, performing a crude Monte Carlo simulation.However, for extremely low probabilities of failure, the huge amount of samples to be evaluated keeps the whole simulation slow.Secondly, to reduce the number of samples, we propose to incorporate subset simulations.
A general issue of data-driven methods is that samples from low-probability regions are not included prior to the response evaluation [15].In this case, events close to or within the region of failure might not be well covered by the POD modes.For example, if plastic deformations did not occur in the observations used to generate the POD modes, they cannot be covered in the approximation.To this end, we will demonstrate the benefits of including POD mode updates in each subset level.
This paper presents the idea of POD in subset simulations and demonstrates its effectiveness in two case studies.The second example focuses on the proposed strategy of updating the POD modes.The broader motivation is to upscale the method to complex and realistic models and loads to estimate the failure probability of buildings based on earthquake excitation.

Probability of failure
To determine the occurrence of structural failure, it is necessary to establish a limit state based on engineering decisions.Assuming that the state of the structure depends on the random variables represented by , a limit state function () is defined.In this definition, failure occurs when the value of the function () is less than or equal to zero, that is, () ≤ 0.
The probability of failure, denoted as (), can be evaluated by integrating the probability density function   () over the failure region where () ≤ 0: In this equation,  represents the failure region defined by the limit state function () ≤ 0. Alternatively, the equation can be reformulated using an indicator function   ( 1 , … ,   ), which equals one if  ∈  and zero otherwise: However, for complex systems, the integral can, generally, not be solved analytically.Therefore,  random realizations (samples) are evaluated such that the probability of failure can be expressed in terms of an expected value: This procedure is commonly referred to as the standard or crude Monte Carlo simulation.While this method yields accurate estimations around the mean of the distribution, the confidence of the estimator depends on the number of samples evaluated.As the estimation approaches the tail end of the distribution, that is, for low failure probabilities, the standard deviation of the estimate increases significantly: In simple terms, obtaining a confident estimate of very small failure probabilities requires a disproportionately high number of samples, resulting in the main drawback of the crude Monte Carlo simulation: the low computational efficiency.

Subset simulation
The objective of subset simulations is to reduce the number of samples included in predicting the failure probability without losing confidence in the estimate.The main idea is to create a decreasing nested sequence of failure regions . The probability of failure is then evaluated using conditional probabilities: Using the conditional probabilities of each subset leads to the estimate of the overall failure probability.To generate samples, we initially draw from a Gaussian distribution.Markov Chain Monte Carlo methods can then be employed to iteratively generate a more representative sample set from the target distribution.After a sufficient number of Markov Chain Monte Carlo steps, the generated samples effectively capture the target distribution and are utilized to determine the first intermediate failure region, denoted as  1 .In this step, a specified portion of the samples, for instance, ( 1 ) = 0.1, which is closest to the failure region, is selected.These selected samples are then utilized to construct the intermediate limit state function   1 ().
The generation of samples in the conditional distribution is a more demanding task.Markov Chain Monte Carlo samplers are required to create the subsets [9].Most commonly, a modified version of the Metropolis-Hastings algorithm [16] is used.However, if complex distributions are involved, using more sophisticated samplers is beneficial.For example, the Hamiltonian Monte Carlo sampler has been shown useful within the framework of subset simulations [17], which can be speed-up by incorporating Hamiltonian neural networks [18].Furthermore, the samples  of the subsets are constrained to stay in the region   , ( |   ) = ()    ()

𝑃(𝐹 𝑖 )
. This can be ensured by rejecting samples that exit the defined intermediate failure region  1 after the Markov Chain Monte Carlo step, for example.

Proper orthogonal decomposition
In this contribution, we focus on problems within structural dynamics.Accordingly, we introduce the POD based on the equation of motion for a nonlinear system, written as: The solution  of the system can be approximated using a reduced subspace by [19]: where   is the reduced basis to be evaluated by the POD.
The POD is an approach that enables finding optimal basis vectors for the transformation into the reduced system.It provides a description of a high-dimensional correlated process in an optimal low-dimensional, uncorrelated space, that is, a random vector  evaluated from a subspace with  dimensions.The approximation using  modes observed using the POD has a lower error than any other orthonormal basis [4]: The objective function for optimal POD basis vectors   = [ 1 ,  2 , … ,   ] is then defined by: To solve this minimization problem and detect   , we introduce an observation matrix   which contains  observations, also called snapshots, of the displacement history () evaluated using the full model: From these observations, a covariance matrix is evaluated: In this case, the internal restoring force must still be evaluated on the full model in each time step.Therefore, the model order reduction is based on a significant reduction of the critical time step if explicit time integration is performed.

RESULTS
The numerical results of this paper were created within an in-house finite element Python code.Furthermore, the subset simulation framework from UQpy was used [20].

Two-story one-bay frame structure subjected to white noise
In the first example, we evaluate our proposed strategy on a two-story one-bay frame structure subjected to white noise.In this example, the randomness is realized within the load, as shown in Figure 1.We generate random white noise excitations in the frequency domain: where  is the spectral density of the white noise and Δ is based on the cut-off frequency.Here, we generated time history data for 5 s using an uncorrelated gaussian distribution, including  = 200 variables  (, ).Therefore,  The structure for this example is a two-stories one-bay building with a total height of 7,m and a width of 5,m, as shown in Figure 1(A).The columns are modeled with a bi-linear elastoplastic material.If   = 235 N mm 2 is exceeded the initial Young's modulus  1 = 210 kN mm 2 is reduced to  2 = 0.1 1 .The material law used within the beams is linear.The largest plasticity occurs in the hinges of the beam-column connections which can be related to the displacement of the adjacent story.Thus, the limit state function is evaluated based on the peak story drift (PSD) of the first story, that is, the vertical displacement of the node in the corner ℎ(): () = PSD  − max |ℎ()| .The result is presented in Figure 1(B), illustrating the maximum displacement max |ℎ()|, that is, the PSD, in the form of the complementary cumulative distribution function.The crude Monte Carlo simulation is performed with the full model using 10 4 samples, which resulted in the estimation of the failure probability of   () = 7 ⋅ 10 −4 and a standard deviation of    = 2.6 ⋅ 10 −4 .Furthermore, two runs of the subset simulation strategy are shown in the diagram.We observe that the results are similar to the crude Monte Carlo simulation.The estimated probability of failure using the subset simulation is   () = 8.47 ⋅ 10 −4 .

3.2
Three-story one-bay frame structure subjected to white noise In this example, the number of stories is increased by one, while the excitation and limit state function from the previous example is kept.However, we also changed the intensity, the material parameter, and the PSD  such that plasticity occurs more rarely.Thus, the failure probability drops to () = 1 ⋅ 10 −5 in this example.
Based on these changes it is likely that plastic deformations are not included within the first realizations, such that they are not covered within the POD modes.As mentioned in the introduction, we, therefore, update the POD modes after reaching a new subset level.The comparison of the displacement history based on two different samples from the last subset is shown in Figure 2.
We observe that the initial modes of the POD cannot approximate the solution of the full model well once large plasticity occurs.In this case, the reduced model gives a lower displacement and, therefore, underestimates the probability of failure if used in the subset simulation or in the crude Monte Carlo method.On the other hand, Figure 2 reveals that the solution of the full model is approximated way better if the updated POD modes are used.
The used modes within the subsets are shown in Figure 3.We choose the number of POD modes included in the reduced system to cover 99% of the energy.Therefore, the number of included modes can vary along the subsets.The first mode is similar throughout all subsets which contains most of the linear displacement.However, the higher modes vary to account for the increasing plasticity within the higher subsets.

Computational efficiency
The incorporation of POD and Subset simulation has led to promising results in terms of statistical results as well as the displacement of the time history, as shown in Figure 2. The main goal of utilizing this strategy is to benefit from the computational efficiency of both methods.We will compare the computational efficiency based on the first example shown in Section 3.1.For our study we split the efficiency into two metrics: (1) the time required to evaluate one sample, and (2) the required number of samples evaluated for the estimation of the failure probability.The speed-up of each sample is significant if explicit time integration, that is, the central difference scheme, is used for both the full and the reduced model.The evaluation of one sample of the full model requires 137.1 s while the reduced model requires 5.56 s only.If the Newmark integration scheme is used for the full model, the computation time can be reduced to 7.48 s.However, the POD is still 1.35 times faster, and for the central difference scheme, it is more than 24 times faster.
Furthermore, we reduced the number of samples by using the subset simulation.For the shown study, we used 10 4 samples within the crude Monte Carlo simulation, while the subset simulation required only 4 000 realizations, to give a reasonable estimate.Therefore, the computational efficiency of the approach is improved by a factor of 2.5.
The overall computational efficiency of estimating the probability of failure is improved by a factor of 3.4 in case the Newmark integration scheme is used for the full model and by 61.6 if the central difference scheme is used for both models Table 1.

CONCLUSION
In this contribution, we provided a model order reduction strategy embedded in a subset simulation framework.Here, we showed the efficiency of the POD.In both numerical examples, we achieved promising results, in terms of the probability F I G U R E 3 Proper orthogonal modes within each subset for the second example, three-story-one-bay structure.The modes are extracted from one random sample within each subset.
of failure, see Section 3.1, and in terms of accuracy of the time history displacement, see Section 3.2.In the latter example, we draw attention to the updates of the POD modes in each subset, which enables accurate model order reduction in low-probability regions.As discussed in Section 3.3, the proposed strategy is more than 60 times faster, if the same integration scheme is used.We conclude that the POD performs well in the subset simulation framework and provides significant speed-up for the chosen problems and settings.We chose a conditional probability of 0.1 for the transition to the next subset and establishment of intermediate failure domains.However, for other models, the reduced model might not adequately approximate the solution for the most severe events, that is, 10%, leading to potential inaccuracies.
In future studies, it is particularly interesting to explore alternative approaches for subset selection based on accuracy while incorporating more complex models and loads.Additionally, analyzing accuracy and speed-up for high-dimensional models will reveal further benefits and limitations of the proposed approach.

A C K N O W L E D G M E N T S
Open access funding enabled and organized by Projekt DEAL.

)
The eigendecomposition of the covariance   gives the POD modes   and POD values   .Another possibility to observe the POD modes and their values is the singular value decomposition of the snapshot matrix:   = SVD(  ), where  contains the POD modes  = [ 1 ,  2 , … ,   ] and  contains the POD values  = [ 1 ,  2 , … ,   ].If nonlinear material behavior is involved -we use a bi-linear elastoplastic model in the following numerical examplesthe nonlinear equation of motion of the reduced system is written as: red  red üred +   red  red ured +   red  int () =   red ().
points are used in the generation of the white noise.
First numerical example: (A) two-story, one-bay structure subject to white noise excitation ü ; (B) complementary cumulative distribution function of the maximum displacement of node two max |ℎ()|.
Comparison of the displacement history based on two different white noise excitations (A) and (B) of the last subset.The evaluation of the full model is shown in blue, while the approximation of the displacement of the POD is shown in orange if the initial POD modes are used and in green if the updated POD modes are used.Computational efficiency of the proper orthogonal decomposition (left table) and the subset simulation (right table).