A hierarchical kriging approach for multi-fidelity optimization of automotive crashworthiness problems

Multi-fidelity optimization schemes enriching expensive high-fidelity functions with cheap-to-evaluate low-fidelity functions have gained popularity in recent years. In the present work, an optimization scheme based on a hierarchical kriging is proposed for large-scale and highly non-linear crashworthiness problems. After comparison to other multi-fidelity techniques an infill criterion called variable-fidelity expected improvement is applied and evaluated. This is complemented by two innovative techniques, a new approach regarding initial sampling and a novel way to generate the low-fidelity model for crash problems are suggested. For the former, a modified Latin hypercube sampling, pushing samples more towards design space boundaries, increases the quality of sampling selection. For the latter, a projection-based non-intrusive model order reduction technique accelerates and simplifies the low-fidelity model evaluation. The proposed techniques are investigated with two application problems from the field of automotive crashworthiness—a size optimization problem for lateral impact and a shape optimization problem for frontal impact. The use of a multi-fidelity scheme compared to baseline single-fidelity optimization saves computational effort while keeping an acceptable level of accuracy. Both suggested modifications, independently and especially combined, increase computational performance and result quality in the presented examples.


Introduction
As computational power has increased exponentially in recent years, also Finite Element (FE) models reached a higher level of detail and complexity-e.g. modern day car models for crash simulations may contain more than ten million elements. This balances out such that simulation times are not significantly decreasing. Especially in multi-query analysis such as optimization or robustness applications a high number of evaluations is required, which increases the computational effort to an infeasible level.
One possible remedy is the use of specifically designed optimization approaches such as Efficient Global Optimization (EGO) that was first proposed by Jones et al. (1998).
The idea is to build a surrogate model from an initial design of experiments (DoE) and adaptively improve it utilizing a so-called infill criterion (Jones 2001;Forrester and Keane 2009). In this context, mostly kriging models (Krige 1951;Matheron 1963;Sacks et al. 1989) are exploited as surrogate models as their inherent error approximation features are especially beneficial.
More recently, these types of surrogate models were integrated in a multi-fidelity scheme, whereby the high-fidelity FE-analysis is complemented with an additional low-fidelity model: The corresponding multi-fidelity kriging schemes can be categorized into two variants. One class of techniques considers correction-based methods, where a "bridge function" or "scaling function" models the differences between high-and low-fidelity models (Chang et al. 1993;Gano et al. 2006;Han et al. 2013). The second type of multi-fidelity approaches are named cokriging. The idea of the latter is to enhance the low-fidelity surrogate model by utilizing the covariance matrix between low-and high-fidelity model. Originally proposed in the geostatistics community (Journel and Huijbregts 1978), this approach was extended to 1 3 114 Page 2 of 15 computer experiments by Kennedy and O'Hagan (2000) and called KOH autoregressive model. Due to its success, numerous extensions and modifications have been added to cokriging since its introduction:  proposed an alternative approach for creation of the covariance matrix between low-and high-fidelity models. Moreover, Gratiet and Garnier (2014) reformulated the cokriging algorithm in a recursive manner to reduce computational complexity. An extension considering complex cross-correlations between varying fidelity models can be found in Perdikaris and Karniadakis (2016). The present work is based on hierarchical kriging (HK), suggested by , whereby the low-fidelity surrogate model represents the trend term in the multi-fidelity predictor. It is beneficial in the context of multi-fidelity optimization as it provides better error estimation capabilities compared to other cokriging models.
As surrogate models were adapted to multi-fidelity applications, so were infill criteria. A criterion named Augmented EI, capable of adaptively sampling low-and high-fidelity models by considering coefficients for cross-correlations and cost ratios between models was suggested by Huang et al. (2006). Moreover, Zhang et al. (2018a) proposed the variable-fidelity expected improvement (VF-EI) criterion that implements a similar idea but with an analytical derivation and free from external coefficients. Therefore, the latter is used in the present work.
A common approach in multi-fidelity optimization is to combine FE models with varying levels of mesh sizes for high and low-fidelity models, such as realized by Zhang et al. (2018a) for an airfoil shape optimization. In combination with a cokriging adaptation presented by Gratiet and Garnier (2014), a hydrofoil shape optimization with varying mesh size levels was performed by Bonfiglio et al. (2018). Similar approaches are investigated in the applications of crashworthiness for honeycomb structures and sheet metal forming in  and , respectively. Alaimo et al. (2018) proposed a multi-fidelity approach where an adaptive functional principal component analysis (PCA) model is utilized with a simulated annealing (SA) algorithm applied to linear elastic structural topology optimization. Anselma et al. (2020) published a multifidelity scheme for the crashworthiness discipline inside a multidisciplinary optimization. The authors use analytical equations as a low-fidelity model and propose to only evaluate the FE high-fidelity model if the former predicts infeasible results. Also a cokriging-based multi-fidelity version of EGO was exploited for inverse problems in haemodynamics (Perdikaris and Karniadakis 2016).
In automotive crashworthiness, optimization has been performed for many years (Redhe et al. 2004;Duddeck 2008;Hunkeler et al. 2013). More recently, multi-fidelity schemes have also been applied in this field . Acar et al. (2020) investigated a multi-fidelity optimization for a frontal impact problem of a bumper system with the multifidelity surrogate modeling approach suggested by Zhang et al. (2018b). Results show that multi-fidelity approaches are capable of yielding significant time-savings while maintaining acceptable accuracy. Other mechanics-based lowfidelity models available for crashworthiness applications are listed in Lange et al. (2018). The authors begin with lumped mass-spring models and subsequently motivate the introduction of the so-called component solution space approach that can be applied in early phase component development for crashworthiness analyses.
Recently, model order reduction (MOR) techniques have been introduced also for non-linear problems (Guo and Hesthaven 2017;Swischuk et al. 2019) and applied in crashworthiness (Kneifl et al. 2021). The non-intrusive approaches are based on the results of training simulationshere named snapshots-which are utilized to compute a reduced subspace. In addition, a regression model is trained that combines the basis vectors of the subspace to represent the physical behavior of the system (Guo and Hesthaven 2019). The non-intrusive MOR has been integrated into a multi-fidelity training scheme by Kast et al. (2020) and related projection-based approaches for crashworthiness applications and optimization have been conducted (Le Guennec et al. 2018;Assou et al. 2019;Gstalter et al. 2020;Ren et al. 2020). A summary of recent developments in the field of non-intrusive MOR is presented in Yu et al. (2019) for fluid mechanics application. Moreover, principal component-based surrogate models can also be found in the field of structural topology optimization (Alaimo et al. 2018;Xiao et al. 2020;Choi et al. 2019).
In the present work we aim to develop enhanced multifidelity optimization schemes in crashworthiness applications. To that end, we propose to integrate an incremental projection-based MOR approach as low-fidelity model into a multi-fidelity EGO algorithm. In a second step to reduce computational effort, our recently developed isovolumetric sampling approach placing samples closer to design space boundaries is adapted (Kaps et al. 2021). When assessing algorithm performance, two main criteria can be established. The primary goal is to find an optimization approach with reduced computational effort produced by the high number of expensive evaluations of the objective function during the optimization. Secondly, an acceptable level of accuracy must be maintained; i.e. a multi-fidelity optimization scheme should not lead to inferior results compared to an optimization based using only high-fidelity simulations.
This work is structured as follows. Initially, the novel design of experiments approach is introduced in Sect. 2, followed by the optimization scheme based on HK and VF-EI in Sect. 3. The MOR approach used for low-fidelity model generations is presented in Sect. 4. Subsequently, the proposed optimization scheme and its implementation Page 3 of 15 114 are explained in Sect. 5. The performance of the complete set of methods is illustrated by a lateral impact example and a crashbox design problem in Sect. 6, and final results are summarized together with an outlook into future work in Sect. 7.

Isovolumetric design of experiments
The first step of any population-based optimization is to determine an initial set of sample points by means of DoE. Covering the full design space with a small amount of samples is the general aim. As there is no unique criterion for this vaguely formulated goal, DoE is still an active field of research. An overview of popular criteria and approaches is given in Garud et al. (2017).
In the present work, a modified optimal Latin hypercube (OLH) approach is used. The standard construction of a Latin hypercube design (LHD) for N samples in d dimensions is described in the following. In each dimension the space is divided into N strata of equal probability, i.e. the design space then consists of N d cells. Randomly, N cells are selected such that each stratum of a dimension may only contain a single sample. Each sample can be placed in the center of its cell or randomly located within it (Rajabi et al. 2015), whereby the centered case is considered in the present work.
Initial Latin hypercube samples are incrementally optimized with a Simulated Annealing algorithm that consists of random pairwise-and coordinate-wise swaps. In an iterative process new samples are accepted if they improve a space-filling criterion or accepted with a certain probability if they do not bring an improvement (Morris and Mitchell 1995). Other approaches with deterministic sample selection, e.g. (Ye et al. 2000) or more elaborate optimization schemes such as the Enhanced Stochastic Evolutionary algorithm (Jin et al. 2003) have been suggested.
Recently, we proposed an adaptation to Latin hypercube sampling, named isovolumetric LHD, that places samples closer to design space boundaries (Kaps et al. 2021). The idea is to rethink the uniform strata that are created in each coordinate dimension for standard Latin hypercube sampling as nested hypervolume shells in the design space. These are created as shown with the colored regions in Fig. 1. Here, all shells are required to have identical volume, which is especially advantageous for higher dimensions, i.e. higher number of design variables. Applying this condition to a d-dimensional unit hypercube, with strata boundaries p i = i N and sizes a j = 1 N for standard LHS, yields the following new equations Here, N is the number of samples to be drawn and N v = N 2 the number of nested hypervolume shells. An exemplary isovolumetric Latin hypercube design for six samples in two dimensions is depicted in Fig. 1. The adaptation is easy to implement and does not increase computational requirements; however, it has been shown to increase the quality of DoEs in popular space-filling criteria such as potential energy (Audze and Eglais 1977) in mid-to high-dimensional situations, i.e. five and more dimensions (Kaps et al. 2021). One expected advantage in the current application is that training points associated to the surrogate models are closer to the design space boundaries. Thereby, the prediction of surrogate models will be more based on interpolation of samples as opposed to extrapolation. This approach, named optimal isovolumetric Latin hypercube (OIVLH), is transferred to the context of optimization in crashworthiness applications in the present work.

Multi-fidelity efficient global optimization
In the following, a multi-fidelity EGO approach is introduced. In general, EGO techniques are based on kriging models, which have first been introduced by (Krige 1951; (1) (2) a j = p j+1 − p j , j ∈ {1, 2, ..., N}.  Matheron 1963;Sacks et al. 1989) and are nowadays a popular choice for surrogate models. The idea of adaptively improving an initially created kriging model by means of an infill criterion was proposed by Jones et al. (1998). Common early variants of the approach are compared in Jones (2001).
Over the years, many surrogate models of variable fidelity have been suggested. Overviews can be found for example in Forrester et al. (2007) or Park et al. (2016). In the present work, hierarchical kriging is applied as the approach has been shown to yield better error estimations than other cokriging methods . The relevant aspects of creating an HK model are summarized below, readers are referred to the original publication for detailed derivations .
The idea is to first create a kriging model of the lowfidelity function that is subsequently used in the hierarchical kriging model for high-fidelity prediction. Based on these surrogate models, new sample points can be evaluated. The models are then adaptively improved. To that end, an infill criterion is introduced that can be maximized to determine the best location and fidelity level for a new adaptive sample.
Following, all steps of the outlined process are explained in details. The low-fidelity function is represented by a normal kriging model. Consider a random process for the low-fidelity (LF) function with 0,lf an unknown constant and Z lf (x) a stationary random process. For a d-dimensional problem, the lowfidelity kriging model is based on a set of sampling data (S lf , y S,lf ) consisting of m lf samples with input variable data S lf ∈ ℝ m lf ×d and corresponding output y S,lf ∈ ℝ m lf .
To be able to construct the kriging predictor for new points, a correlation function, also named kernel, is needed to model the correlation between given sample points and new points. An overview of popular choices is for example given in Rasmussen and Williams (2005). In the present work, a squared-exponential kernel, also called Gaussian radial-basis function (RBF) kernel is utilized due to its smoothness and infinite differentiability: Here, k denotes the kernel length scale. The kernel is called anisotropic, if there is a separate length scale for each design space dimension as in the equation above. In the present work, an isotropic kernel is chosen, where the hyperparameter k = is a scalar, i.e. independent of coordinate dimension. In fact, many different RBF kernels with similar properties than the Gaussian kernel exist. The focus in the present work lies solely on the latter for the sake of clarity.
The low-fidelity predictor for a new point x can then be written as where r lf is the correlation vector between the sample data and the new point, R lf ∈ ℝ m lf ×m lf the correlation matrix representing correlation between sample data points and 1 ∈ ℝ m lf a column vector filled with ones.
Following the calculation of the initial sample data set, the kriging model is fitted by separately optimizing the kernel hyperparameter . Differential evolution (Storn and Price 1997) is selected for the optimization of the hyperparameters in the present work due to its simplicity and its good global search characteristics. More on hyperparameter optimizations can be found in Toal et al. (2008).
Next, the hierarchical kriging model can be constructed including the predictor of the low-fidelity model. Therefore, consider a random process corresponding to the high-fidelity function Here, the low-fidelity predictor scaled by an unknown constant 0 is a trend term and Z(x) is a stationary random process. Given a d-dimensional sampling data set (S, y S ) containing m samples with input variable data S ∈ ℝ m×d and corresponding output y S ∈ ℝ m , the HK predictor for the high-fidelity function can be written as where 0 is a scaling factor indicating the correlation between high-and low-fidelity functions and F = [ŷ lf (x (1) )...ŷ lf (x (n) )] T , ∀x (i) ∈ S represents the low-fidelity prediction at high-fidelity sample points. r ∈ ℝ m and R ∈ ℝ m×d are defined the same way as for the low-fidelity predictor above. The factor R −1 (y S − 0 F) , named V HK in the original publication, does not depend on the untried point x and can thus be calculated at model fitting time. The mean squared error (MSE) of the HK prediction is given with 2 , the process variance of Z(x) Once the initial HK model is built, it is adaptively improved using an infill criterion to determine the ideal position of new samples. In multi-fidelity applications two options exist: a 'classic' single-fidelity infill criterion or a multifidelity criterion. Among the former, expected improvement . (Jones et al. 1998) is the most popular method. The interested reader is referred to Jones (2001) for other common techniques. The disadvantage of a single-fidelity criterion is that only high-fidelity samples can be considered adaptively. Hence, they are not further discussed here. In the present work, the so-called variable-fidelity expected improvement criterion (Zhang et al. 2018a) is utilized. It is defined at location x and fidelity level L as s(x,L) and y min as the currently best observed feasible high-fidelity function value. The term s(x, L) denotes the uncertainty of the HK model. The previously introduced scaling factor between fidelity levels 0 is used here to model the uncertainty in high-fidelity prediction caused by the lowfidelity predictor Here, MSE(ŷ(x)) and MSE(ŷ lf (x)) are the MSEs of the highand low-fidelity kriging predictors, respectively. (•) and (•) in Eq. (9) represent the cumulative distribution and probability density functions of the standard normal distribution, respectively. The two summands in Eq. (9) can be identified with exploration and exploitation. The first term y min −ŷ(x) (u) is dominated by improving the solution ŷ(x) and thus represents exploitation, while the second term s(x, L) (u) represents exploration because it is dominated by the solution uncertainty s(x, L).
Due to the highly multi-modal nature of the EI functions, differential evolution (Storn and Price 1997) is selected for optimization of the infill criterion in the present work.
Notably, the VF-EI formulation is similar to the original EI definition (Jones et al. 1998); however, in addition the dependency on the fidelity level is introduced. In terms of multi-fidelity optimization the VF-EI criterion is comparable to the augmented EI criterion proposed by Huang et al. (2006). Both describe the expected improvement of the high-fidelity function with respect to adaptive samples on both fidelity levels. To that end, augmented EI contains two factors that are multiplied with a standard EI for the high-fidelity function. A more detailed discussion about the differences between VF-EI and augmented EI can be found in Zhang et al. (2018a). Here, we use VF-EI as it is free of empirical parameters and is as such more intuitive.

Model order reduction
The proposed multi-fidelity approach exploits a data-driven model order reduction (MOR) technique to create the lowfidelity model. As an analytical simplification is not available for non-linear problems, a data-driven approach is commonly used (Sirovich 1987). Therefore, an online and offline phase are introduced, whereby the offline phase can be understood as the counter part to the DoE. In particular, the surrogate model is created during the offline, also called training phase. Afterwards the simplified model can be evaluated for multi-fidelity analysis in the online phase.

Training phase
Within the training phase, a set of full order simulations is created, whereby all resultants are stored as snapshots x i ∈ ℝ N with N degrees of freedom. Combining the training simulations to a so-called snapshot matrix A ∈ ℝ N×n , with n the number of collected snapshots, a reduced subspace and its projection matrix can be computed. Through the Singular Value Decomposition (SVD), also referred to as thin SVD (Golub and Van Loan 2013) the snapshot matrix A can be represented by the left-singular vectors U ∈ ℝ N×n , the diagonal matrix Σ ∈ ℝ n×n containing non-negative singular values i in descending order and the right-singular matrix Z ∈ ℝ n×n . Thus, the columns of the matrix U are the eigenvectors of AA T .
Moreover, the matrix A is approximated by truncating its parts to a rank k ≤ n, m , such that U k ∈ ℝ N×k , Σ k ∈ ℝ k×k and Z k ∈ ℝ k×N , respectively. To define the reduced basis of the subspace, V ∶= U k ∈ ℝ N×k is further utilized as the projection matrix. In practice, the optimal rank k is not known beforehand and k = min k with can be found for an error threshold . In other words, the matrices are truncated by k such that an energy cutoff ratio is maintained: For large-scale matrices the evaluation of the full SVD is cost intensive as its complexity is in the range O(n 2 ) with n as the number of snapshots. Therefore, multiple approaches (Bach et al. 2019;Phalippou et al. 2020) to efficiently compute the truncated projection matrix, such as randomized or incremental SVD techniques, e.g. (Oxberry et al. 2017), can be applied. Here, an incremental SVD algorithm (Baker et al. 2012) 1 is utilized. To save computational resources, the snapshot matrix A is built up incrementally. Therefore, A is divided into batches, which are added to the projection matrix V . The SVD is computed within every iteration via QR decomposition and its truncation rank k is evaluated. The batch is added to the global V and Σ before the algorithm steps into the next iteration with new snapshots. Readers are referred to Baker et al. (2012) for a detailed algorithm description. In summary, a complexity of O(mnk) with m full order unknowns can be reached.
After the construction of the reduced subspace a regression or interpolation approach is introduced. The metamodel represents the unknown in the reduced space and is therefore restricted to the physical solution spanned in the subspace. In addition, the number of unknowns n is drastically reduced, as k << n . Within this work the k-nearest neighbor (kNN) approach is utilized as an interpolation technique, but any other machine learning approach such as polynomial regression function, Gaussian process regression or neural networks could be used (Swischuk et al. 2019). k-nearest neighbor is mainly known as a classification technique, but can also be applied as a regression model. The function y = f (x) is interpolated by the k-nearest neighbors of x as shown in Fig. 2 for a one-dimensional example. Here, three neighbors of x and their distances are evaluated to estimate f(x).

Online phase
After the construction of the low-fidelity model is completed, it is evaluated in the online phase. Recalling the truncated singular value decomposition, the columns of the matrix V ∈ ℝ N×k can be interpreted as the basis vectors v of the subspace.
The full displacement vector is estimated by the linear combination of the basis vectors v , whereby every basis is weighted by a scalar value i .
The kNN approach provides an estimation for each i , the degrees of freedoms of the subspace.

Proposed optimization scheme
In the present work the performance of multi-fidelity optimization schemes in crashworthiness applications is investigated. To that end, a multi-fidelity optimization method that integrates a non-intrusive reduced order model into the hierarchical kriging surrogate is proposed. The schematic process of this approach is shown in Fig. 3 following the baseline multi-fidelity scheme proposed by Zhang et al. (2018a). Initially, DoE is performed as described in Sect. 2 on both, the high-and the low-fidelity level separately, usually generating significantly more low-than high-fidelity samples. OLH and OIVLH are both applied in the present work, to assess the impact on optimization performance. All highfidelity samples are then calculated. For memory efficiency, a reduced order model is incrementally created during the high-fidelity evaluations as introduced in Sect. 4. Following all initial high-fidelity simulations and as the main adaptation to the originally proposed scheme, the reduced basis is evaluated and the k-nearest neighbor regression model is trained for predictions. The reduced order model is evaluated on the initially created low-fidelity DoE points. From the results, the initial low-fidelity kriging model and subsequently the hierarchical kriging model are fitted. For adaptive improvement, the infill criterion (i.e. VF-EI) is maximized separately on both fidelity levels. Depending on which level yields the better results the next adaptive sample can be either low-fidelity (L=0 in Fig. 3) or high-fidelity (L=1). The objective function is evaluated for the respective new adaptive sample and the kriging model(s) are updated. Then, another infill criterion optimization is started and the iterative improvement continues. Three different criteria determine the termination of the algorithm. First, a minimum allowable value for the maximized infill criterion is specified (here IC th = 10 −5 ). Second, a maximum number of high-fidelity evaluations is specified. Third, a maximum total number of objective evaluations can be specified, which in single-fidelity optimizations is equivalent to the second criterion. Values specified for the criteria are given in Sect. 6 with the respective examples. The first criterion can be interpreted as convergence of the algorithm to an optimal point with little improvement possibilities. The other two criteria are used to represent time restrictions on the optimization runs, i.e. an estimation of the maximum run time.
The proposed scheme is implemented in an in-house Python code. The DoE part of the algorithm is based on Kaps et al. (2021). The implementation of the incremental SVD technique in the present work is adjusted from Baker et al. (2012) and Bach et al. (2019). HK model generation and kernel implementation are based on the scikit-learn library in Python (Pedregosa et al. 2011).

Crashworthiness application problems
In the following, the presented optimization scheme is compared to a baseline multi-fidelity scheme proposed by Zhang et al. (2018a) as well as a single-fidelity scheme based on a high-fidelity model and EI. Additionally, each of the techniques is assessed with a standard OLH sampling as well as the OIVLH method that places samples closer to design space boundaries. These overall six approaches are each evaluated for two crashworthiness problems in terms of result quality and computational requirements. Each optimization run is repeated ten times to ensure reliability of the assessment. An overview of the compared methods and nomenclature is given in Table 1. All analyses are performed on the identical hardware using the explicit finite element software LS-Dyna in its MPP version distributed on eight cores. All objective function values referred to below are based on the high-fidelity model unless explicitly stated otherwise.

Side sill impact problem
The initial problem is a crash model representing the side sill of a car under side pole impact as depicted in Fig. 4. Both ends of the side sill are fixed, and the pole represented by a cylindrical rigid body with radius 35mm has a prescribed initial velocity of 36 km/h and a mass of 86 kg. Further modeling information as well as all material parameters are summarized in "Appendix A". The simulation is terminated when the impactor stops or-as a backup-after 40 ms.
The design variables of the optimization problem are the thicknesses t i of the five horizontal reinforcement ribs in the interior of the side sill (compare detail in Fig. 4). The objective is to minimize the mass of the side sill, i.e. the mass of the horizontal ribs, while keeping the lateral intrusion below u allow = 50 mm. Applying the penalty method with a penalty factor p = 3.75 , the deformation constraint is included into the objective function. The complete optimization problem can then be formulated as Fig. 3 Schematic representation of the proposed optimization scheme to integrate reduced order model (ROM) into a multi-fidelity optimization based on hierarchical kriging and variable-fidelity expected improvement (compare Fig. 2 in Zhang et al. (2018a)). The main difference to the originally proposed scheme is the utilization of initial high-fidelity sample data for ROM creation (marked by the grey box). Not shown is the varying method of the design of experiment performed in the present work The model depicted in Fig. 4 represents the high-fidelity model with an element size of approximately 5mm in each direction. For comparison, the low-fidelity model with the element size doubled to 10mm is depicted in "Appendix A". Thus, the speed-up factor of the low-fidelity model is about five to six. The number of initial high-fidelity samples for the single-fidelity case is set to 50. For all multi-fidelity approaches, the number of initial samples is 20 and 120 for high-and low-fidelity models, respectively. In this example, the only termination criterion specified is a threshold value for the respective infill criterion, i.e. EI for SF methods and VF-EI for MF methods, of IC th = 10 −5 . The reduced order model is constructed with 20 training simulations sampled by the respective method. With the truncation energy of = 0.9999 a subspace containing approximately k = 4 basis vectors is computed. Moreover, the kNN regressor using 5 neighbors is trained for the unknowns in the reduced subspace. Eventually, a single evaluation of the reduced order model has a speed-up factor of 150 in comparison to the high-fidelity simulation.
In Fig. 5, optimization results for the six different approaches are compared with a parallel coordinates plot.
On the x-axes all five design variables are listed while the y-axes represent their normalized ranges. Each curve illustrates an optimized design, whereby the color indicates the respective objective function value. The problem seems to have a clear global optimum with t 1 = t 3 = t 5 = 3.0 mm and t 2 = t 4 = 0.5mm. Distributing three horizontal ribs with the maximum allowable thickness across the whole height of the component and having the two ribs in between vanish results in the best compromise between low weight and sufficiently low impactor intrusion. The respective objective value f (t) = 1.898 is found in the majority of single-fidelity and MF (MOR) repetitions independent of the chosen DoE method. For this set of input variables, the maximum impactor displacement of the high-fidelity model is u max = 50.6 mm. The specified constraint u allow = 50 mm is violated by 1.2%, which is considered acceptable in this exemplary problem. All compared results have a similar level of constraint violation, which is not further investigated to focus on the multi-fidelity optimization itself.
Individual results for single-fidelity, i.e. high-fidelity, EGO show a low variation for different DoE methods and the majority of repetitions terminate at the global optimum. For MF (base) and MF + OIVLH, results vary quite significantly and only three of the overall 20 runs converge to the global optimum. Many of the evaluations terminate at points rather close to the global optimum so that the objective function value is only a few percents of the optimal value. Notably, more MF + OIVLH analyses get close to the global optimum than MF (base). The MF (MOR) approaches yield more consistent results than MF (base) and MF + OIVLH, converging to the optimum for a total of 14 out of 20 assessments across both DoE methods, with two more runs being very close to the optimum. Overall, the OIVLH-based approaches yield a higher consistency within the multi-fidelity schemes. OIVLH places samples significantly closer to design space boundaries. Therefore, the reduced order low-fidelity model being created from the high-fidelity samples is based more on interpolation between samples than on extrapolation.
In a further step, the quality of the two varying types of low-fidelity models is assessed. To that end, the low-fidelity  models are evaluated at the respective determined optimum. The maximum displacement of the impactor d max is considered as a metric to compare high-(HF) and low-fidelity (LF) models. This metric is chosen as it directly represents the constraint and-because the mass of the component can analytically be calculated from the given design variablesit also implies the objective function of the optimization problem. As such it estimates the accuracy of the low-fidelity models in the present example. The error metrics are defined as for the absolute difference e abs and the relative difference e rel between fidelity levels, respectively. Results for all methods are listed in Table 2. Each of the listed values represents the mean value of ten evaluations. For MF (base) and MF + OIVLH the low-fidelity model is identical for all evaluations, while for the MOR based methods, the low-fidelity model depends on the initial high-fidelity samples (compare Fig. 3). Both types of low-fidelity models approximate the (maximum) impactor displacement sufficiently well with an error of 4-6% . The error indicators for the reduced order lowfidelity models are slightly better, i.e. lower, compared to those of the coarse simulation models. For the latter, the use of OIVLH sampling does not significantly impact the lowfidelity result quality. However, for MF (MOR) + OIVLH the low-fidelity model is slightly more accurate than for MF (MOR) which may explain the performance difference between the two methods. It also indicates that the increased interpolation share for OIVLH sampling can in fact increase the ROM quality.
Keeping in mind the quality of the results, the computational requirements for all approaches are evaluated. Average run times along with their respective standard deviations are listed in Fig. 6 for all techniques. The larger standard deviations for SF (base) and SF + OIVLH compared to the other four methods are explained by a single outlier in each of the (17) e abs = |d max,hf − d max,lf |, methods. In both cases, the algorithm finds the global optimum in a reasonable number of iterations but then requires many adaptive samples to reach the specified termination criterion. So both outliers can be explained by the somewhat unrealistic definition of termination criterion here, where no maximum allowed number of iterations was specified. As the standard deviations in all methods are similar apart from that, the following comparisons are focused on the average times. SF (base) is taken as a baseline for all comparisons. Without changing the optimization itself, switching the DoE to OIVLH reduces 14% of computation time in this problem, while maintaining result quality.
Using an MF approach with a coarser low-fidelity model lowers the computational cost about 47% and 35% for OLH and OIVLH, respectively. In fact, this is the only case in the present work, where OIVLH requires more computational time than OLH. The reason here is the same as for the outliers outlined above. Both MF (base) and MF + OIVLH find the respective optimum after a similar number of highfidelity iterations, but the latter variant requires a few more adaptive iterations than the former to reach the termination criterion. The work load to create the low-fidelity model is not considered here, however for complex models this is an additional time intensive step. Especially, when considering that the low-fidelity model may not be needed otherwise. Contrary to that, both MF (MOR) options do not require the manual creation of a low-fidelity model. They reduce computation times by about 51% and 56% for MF (MOR) and MF (MOR) + OIVLH, respectively. The speed-up of MF (MOR) compared to the other MF techniques can be explained by the significantly faster evaluation times of the reduced order low-fidelity models compared to the coarse simulation model. In comparison the MF (MOR) + OIVLH approach yields the best overall improvement for the side sill impact problem. It saves on average more than 50% computation time compared to SF (base) while maintaining an acceptable level of accuracy.

Frontal impact problem
A shape optimization problem of the frontal impact of a crash box, as depicted in Fig. 7 is presented as a second application problem for the suggested multi-fidelity scheme.
For the high-fidelity model, an element length of about 4 − 5 mm is specified, while for the low-fidelity model, the element length is in the region of 10mm. This yields a highfidelity model with 4,928 elements and a low-fidelity model with 1,296 elements. The planar impactor crushing the component from the top has a prescribed mass of 300kg and an initial velocity of 30 km/h. The impactor is modeled as a rigid body, and the crash box is constructed as a tube with a steel-like material and a piecewise linear plasticity model. Exact material parameters are listed in "Appendix B". For the contact formulation, the LS-Dyna '*CONTACT_ AUTOMATIC_SINGLE_SURFACE' is applied. The simulation is terminated when the impactor stops or after 45ms. The crash box of the optimization study and its design variables are depicted in Fig. 7. The first three design variables are the vertical positions of the triggers in the model. The corner elements of the triggers are deleted to increase numerical stability. Three additional design variables are the depths of the triggers, i.e. the respective element rows are shifted in or against the arrow directions indicated in the figure next to the variables x 4 , x 5 and x 6 . The configuration depicted in Fig. 7 represents the setup for x = 0 for all design variables. Model generation features are realized using ANSA preprocessor.
Throughout the optimization analyses, the mass of the crash box remains approximately constant as the effect of design variables is negligibly small. The objective function is chosen as the load uniformity, also called peak amplification, of the force-displacement curve. It is defined as the peak force F max divided by the mean force F mean of the complete force-displacement curve measured at the rigid body impactor. The optimization problem can then be formulated as The processing time of the classical low-fidelity model is one fourth of the high-fidelity computation time. The number of initial samples is accordingly set to 60 in single-fidelity analysis, and 30/90 for high-/low-fidelity samples in the multi-fidelity methods. For this problem, all possible termination criteria presented in Sect. 5 are specified and respective values are listed in Appendix Table 7. To construct the snapshot matrix for the reduced order model, 30 initial high-fidelity samples are utilized. With the truncation energy of = 0.9999 , a subspace with approximately k = 15 bases is computed. Moreover, the kNN regressor using 5 neighbors is trained for the unknowns in the reduced subspace. A single evaluation of the reduced order model has a speed-up factor of 50 to 100 in comparison to the high-fidelity simulation.
To investigate the convergence of the optimization schemes, Fig. 8 shows the best current objective values in each iteration for all (high-fidelity) evaluations of SF (base) (grey lines) and MF (MOR) + OIVLH (black lines). These two methods are representatively chosen from all investigated methods for the purpose of clarity. The adaptive phase of the algorithm starts after 60 evaluations for the single-fidelity case and 30 evaluations for the multi-fidelity case. Both approaches reach objective values below 3.4, i.e. close to the final optimum, mostly within ten adaptive high-fidelity evaluations. Afterwards only slight improvements are achieved.
In addition, the termination criteria can give valuable insights to the performed optimization studies. Here, the observations are contrary to those of the previous lateral impact example, where all optimization analyses are Frontal impact problem of a crash box impacted by a planar rigid wall: The lower end is fixed, the impactor moving from the top downwards has a prescribed mass and velocity. Design variables are the vertical trigger positions as well as depths of the three triggers in the crashbox. The latter are realized by shifting the respective element rows in or against the direction indicated by the respective arrows next to x 4 , x 5 and x 6 terminated by the threshold infill criterion. A majority of analyses for the crash box example terminate after reaching the maximum number of allowed iterations. Of course, there is an argument to be made that these optimization runs may not be fully converged. However, as seen in Fig. 8 and because some runs do actually terminate due to IC th , the authors believe that this is not the case. A maximum allowable number of iterations also represents a situation more akin to a practical application case. It matters more to find an acceptable result in a given time as opposed to the best result in as much time as necessary.
A comparison of the optimization results for the presented methods is depicted in Fig. 9. Low objective function values of slightly below 3.3 depend on x 1 on its upper limit of 10mm, x 2 around its lower limit of −10 mm and x 3 in a mid region albeit slightly above 0mm. Only the first of the trigger depths x 4 , x 5 and x 6 , which is found close the maximum of 4mm throughout all methods, seems to significantly impact the optimization results. A significant share of the crash energy is absorbed in the first fold, i.e. the one controlled by x 1 and x 4 . It usually also contains the total maximum in the force-displacement curve F max . To that end, it appears reasonable that those two design variables appear to always converge to similar values, while the others, especially x 5 and x 6 , may vary between optimizations.
All compared methods yield rather similar results in terms of objective function values except for a total of three runs in MF (base) and MF (MOR). In these, the best objective value is up to 10% off the best overall result because the algorithm terminates in a local optimum. For the SF approaches, OIVLH appears to not affect result quality or consistency. In the multi-fidelity setups however, OIVLH sampling reduces the number of outliers compared to MF (base) and MF (MOR), respectively. All variants utilizing OIVLH sampling as well as the baseline single-fidelity EGO yield very similar results.
In a further step the quality of low-fidelity models is evaluated as done for the previous example. Here the same error metric is chosen mainly for two reasons. First, as the impactor kinetic energy remains constant, the maximum impactor displacement is directly related to the mean force F mean which is part of the objective function. Second, having the same comparison metric as in the previous example allows for comparisons of low-fidelity model accuracy between the two examples. Table 3 lists mean values for the differences in maximum impactor displacement between low-and high-fidelity models in the different methods [compare Eqs. (17) and (18)]. In this example, the coarse simulation model of MF (base) and MF + OIVLH has a significantly larger error value compared to the ROMs. Due to the highly nonlinear nature of frontal impact simulation, the mesh size has a high impact on the crushing behaviour of the component. As the ROM-based low-fidelity models utilize the high-fidelity model mesh for learning and predictions, they are not affected and show better accuracy compared to the coarse low-fidelity model. The present example confirms that the application of OIVLH for reduced order models increases the prediction accuracy. The results of MF (base) and MF + OIVLH, as listed in Table 3 differ mainly due    (17) and (18) to one single outlier in MF + OIVLH, where the low-fidelity model predicts a significantly earlier stop of the impactor, thus inflating the error measure. Overall, the low-fidelity accuracy of the two methods is similar, which is expected as the choice of DoE should not have an impact on the simulation model. Computational requirements are again compared for all different methods using the mean run time of each optimization. The averaged processing times for all techniques are plotted in Fig. 10. Regarding the varying EGO approaches, the findings here are highly similar to those of the previous problem. Compared to SF (base), MF (base) and MF + OIVLH yield about 30-45% speed-up while MF (MOR) and MF (MOR) + OIVLH result in 50-55% time reduction. Here, the use of OIVLH over OLH shows slight time benefits across all variants. Notably, the time benefit SF + OIVLH yields over SF (base) is significantly smaller than for the previous example (compare Fig. 6). This is due to the difference in applied termination criterion between the examples.
Overall, results for the present application example closely match those for the side sill impact example presented above. The MF (MOR) + OIVLH scheme performs best with timesavings of more than 50% compared to SF (base). Here, all approaches utilizing OIVLH perform better than respective OLH variants.

Conclusions
In the present work, a multi-fidelity efficient global optimization based on recently proposed hierarchical kriging and variable-fidelity expected improvement is applied to crashworthiness examples and its performance is investigated. Additionally, two adaptations to the scheme regarding initial design of experiments and the choice of lowfidelity model are proposed. For the former, a recently developed variant of Latin hypercube sampling is chosen that places samples closer to design space boundaries and thus allows for more interpolation instead of extrapolation in surrogate models. For the latter, a non-intrusive model order reduction scheme is applied as low-fidelity model as it integrates nicely into the existing multi-fidelity optimization scheme. All different optimization schemes are investigated on two crashworthiness application examples: one side sill impact size optimization and one frontal impact shape optimization. Already with the rather small problems presented here, results show that multi-fidelity optimization is capable of reducing computational costs of the optimizations significantly while not compromising result quality. Both proposed adaptations independently and especially combined further reduce computation times and also increase result quality compared to the baseline multi-fidelity optimization. Especially the use of nonintrusive reduced order modeling techniques is promising as it removes the need to (manually) create an additional low-fidelity model. Together, a speed-up factor of two in the optimization with next to no influence on result quality is observed.
The problems shown in the present work represent rather small examples with a low number of design variables. Based on the works introducing OIVLH (Kaps et al. 2021) and previous investigations on the projection of large-scale systems (Bach et al. 2019), it seems reasonable to assume that the advantages of the proposed schemes grow as the model size and number of design variables increase.
Multi-fidelity optimization is a wide topic with a variety of different applications and many imaginable adjustments to be explored. Based on the promising results of the present work we have collected some topics and questions that we believe to be interesting for future work: -In the present work all DoEs are performed separately for different fidelity levels and with no connection between levels. It seems only reasonable to use a multifidelity DoE scheme if the initial samples are combined into a multi-fidelity surrogate. Multi-fidelity sampling approaches proposed so far, require the high-fidelity DoE to be a subset of the low-fidelity DoE. It could be investigated, how these approaches perform for multifidelity optimizations shown here and if methodological improvements can be achieved. -We believe the potential of the proposed multi-fidelity scheme(s) should be confirmed in further studies on more complex larger crashworthiness application problems and other fields of applications. -A big challenge in practical applications is robustness with regards to both the method as well as the objective function. An optimization method should produce consistent results for given inputs, as in practice, repeating runs is often infeasible. Moreover, an optimum highly sensitive to small perturbations of the inputs is also not Fig. 10 Frontal impact problem: Average run times for all compared optimization approaches in seconds. Also shown is the standard deviation of the optimization run time desirable. To that end, an effort has to be made to integrate robustness into the optimization framework seamlessly.

Appendix B
See Tables 6 and 7.