A parametric and non-intrusive reduced order model of car crash simulation

Abstract Industrials have an intensive use of numerical simulations in order to avoid physical testing and to speed up the design stages of their products. The numerical testing is indeed quicker to set-up, less expensive, and supplies a lot of information about the system under study. Moreover, it can be much closer to the physical tests as the computation power increases. Despite the rise of this power, time consuming simulations remain challenging to be used in design process, especially in an optimization study. Crash simulations belong to this category. These rapid dynamic computations are used by RENAULT during the sizing of the vehicle structure in order to ensure that it meets specifications set up to reach safety criteria in case of accidents. They are completed using finite element software such as VPS (Virtual Performance Solver) developed by ESI group that will be used in this study. For car manufacturers, the goal of the optimization study is to minimize the mass of the vehicle (and thus its consumption) by modifying the thicknesses of some parts (from 20 to 100 variables). Industrials such as RENAULT currently perform optimization studies based on numerical design of experiments. The number of computations required by this technique is from 3 to 10 times the number of variables. This is too much in order to be intensively used in a design process. In order to decrease the time-to-market and to explore alternative technical solutions, we explore the potential of using a parametrized reduced order model in the optimization studies. The parametrized reduced order model gives an estimation of the high-fidelity result for a new set of parameters without using the solver, by analysing the existing results of previous computations with various sets of parameters. The developed reduced order model is called ReCUR. It is partly based on a CUR approach embedded in a regression analysis. The regression statistical model uses the data of a few calculations made with the solver. Other tools such as clustering and linear programming are used to get the regression analysis more efficient. It is hoped to drastically reduce the number of required simulations of a standard optimization study. In this paper, the construction of the reduced order model will be presented. Then, the relevancy of using the reduced order model into a design process will be exhibited through the treatment of two industrial test-cases. Some improvements of the method as well as several potential uses will then be outlined. The applications will highlight the promising power of the method to shorten design process using optimization and long-run simulations.


Abstract
Industrials have an intensive use of numerical simulations in order to avoid physical testing and to speed up the design stages of their products. The numerical testing is indeed quicker to set-up, less expensive, and supplies a lot of information about the system under study. Moreover, it can be much closer to the physical tests as the computation power increases. Despite the rise of this power, time consuming simulations remain challenging to be used in design process, especially in an optimization study. Crash simulations belong to this category. These rapid dynamic computations are used by RENAULT during the sizing of the vehicle structure in order to ensure that it meets specifications set up to reach safety criteria in case of accidents. They are completed using finite element software such as VPS (Virtual Performance Solver) developed by ESI group that will be used in this study. For car manufacturers, the goal of the optimization study is to minimize the mass of the vehicle (and thus its consumption) by modifying the thicknesses of some parts (from 20 to 100 variables). Industrials such as RENAULT currently perform optimization studies based on numerical design of experiments. The number of computations required by this technique is from 3 to 10 times the number of variables. This is too much in order to be intensively used in a design process.
In order to decrease the time-to-market and to explore alternative technical solutions, we explore the potential of using a parametrized reduced order model in the optimization studies. The parametrized reduced order model gives an estimation of the high-fidelity result for a new set of parameters without using the solver, by analysing the existing results of previous computations with various sets of parameters. The developed reduced order model is called ReCUR. It is partly based on a CUR approach embedded in a regression analysis. The regression statistical model uses the data of a few calculations made with the solver. Other tools such as clustering and linear programming are used to get the regression analysis more efficient. Introduction Despite the continuous increase in computational power, the speed of the numerical simulations does not follow the same trend on account of the refinement of the numerical models and algorithms complexity. For instance, a crash simulation remains a time-consuming simulation (15h on 72 processors for one simulation). During the design phase, numerous crash simulations have to be run in order to find innovative designs fulfilling the safety requirements as well as improving comfort and performances. Several approaches may be investigated to shorten the design phase, and thus improve the time-to-market. Their goal is to decrease the simulation time building equivalent models and thus to limit the calls to time-consuming models. They may be divided in two categories, the reduced order model and the metamodeling techniques.
The first group gathers the reduced order model techniques also called projection-based method. The simulation time is shorten using appropriate bases. In the Proper Orthogonal Decomposition (POD) [1], the bases are built from the diagonalization of the correlation matrix of snapshots. In this sense, it is strongly related to the Principal Component Analysis [2] or to the Singular Value Decomposition also called discrete Karhunen-Loève transform. This kind of method gave a lot of satisfactory results and extensions. For instance, both the Discrete Empirical Interpolation Method (DEIM) [3,4] and the A Priori Hyper Reduction (APHR) [5] proposed methods based on the Gappy POD method [6] in order to select appropriate snapshots leading to satisfactory reduced basis. These methods allow the use of partial or incomplete data (sample mesh in the DEIM and truncated integration domain in the APHR) to build a reduced order model. It is proposed in [7,8] to identify the parts having a non-linear behaviour and the parts having a linear behaviour using a clustering technique and to solve the non-linear problem with the DEIM and the linear part with Krylov subspaces [9]. As well as the APHR [10], this recent method has been applied to crash simulation but it seems limited to particular parts or small models with a low number of parameters. The Gappy-POD has been developed originally for image processing but has been successfully applied to a number of field [11][12][13]. The A Priori Reduction (APR) [14] is strongly related to the APHR method. It tries to adapt and correct the reduced order model according to the quality of the approximation of the POD method using an iterative process. Moreover, the APR tries to improve the basis by the use of snapshots from alternative simulations. This point is very interesting for industrial applications because reduced order model may benefit from previous studies rather than keeping them underexploited. Among these methods, the Proper Generalized Decomposition [15][16][17] performs a variable separation between the time, the space, and the parameters. The advantage of this method is the possibility to perform quasi real time approximation of simulations with new sets of parameters once the reduced order model is built. Drawback of this method stands in the separability hypothesis. It is rarely verified in highly nonlinear models. Crash simulation is typically this kind of model because it gathers high displacements, contact and material non linearities. Besides, most of these methods require calls to the solver to perform the time integration and sometimes modifications of the governing equations of the simulation model. Unfortunately, such an intrusive method is hardly practical in existing design processes because the modifications of the solver are restricted by property and certification considerations.
Contrarily to the projection based techniques, the metamodeling techniques do not require calls to the solver of the governing equations to perform an interpolation. They make use of the simulations to build statistical response surfaces of some quantities of interest. These quantities are often directly extracted from the design problem specifications. They come from some postprocessings of the high-fidelity simulations. The statistical responses may be interpolation models as the usual kriging or regression models as polynomial chaos expansions [18,19], radial basis function [20], or artificial neural network [21]. These models has been already applied to the resolution of optimisation problems in various physics [22][23][24][25][26] giving satisfactory results. Kriging has also been used to mimic some outputs of crash simulations [27]. Advantages of such method is that they do not require to call a specific solver and thus, they can be easily adapted to alternative kinds of physics. Despite they may give very satisfactory solutions for the quantities modelized by the surface responses but they may lead to unphysical solution for the rest of the system that are not taken into account in the response surfaces. Thus, it may be very favourable to be able to approximate the entire crash simulation rather than just few post-processed quantities in order to show the approximation of the entire simulations to an expert of the emulated physics. To solve the optimization problems, the post-processing of the quantities from the specifications may be done on the approximation of the entire structure.
Moreover, additional specifications may be easily added into the optimisation problems without the need of computing another interpolation function. Finally, such global metamodels may benefit from the set of the simulations done in past studies for its construction, allowing a more efficient management of the data. The link between such kind of method and reduced order models are the low rank approximations such as Tucker decomposition [28], High Order SVD [29], Train Tensor Cross decomposition [30], and CUR or Nystrom decomposition [31]. The interpolation of a new numerical experiment will indeed be performed using such a reduced order model together with a parametric function. In this paper, an approach gathering a low rank approximation for the space and time decomposition with a regression model for the parametric interpolation will be outlines and applied to crash simulation. Once the off-line stage of construction of the reduced order model performed, the interpolation for a new parameter set will be done using this reduced order model without any call to the finite element solver. In this, this work may recall the method used in [32] for solving aerodynamics problem by coupling clustering and interpolation on the variety of the orthogonal matrices. Hereinafter, they will be called parametrized non-intrusive reduced order models. Such parametrized reduced order model are at the crossroads of metamodeling, reduced order model, low-rank approximation, and machine learning [33]. Beyond optimization studies, such models may be useful for sensitivity analysis or control theory that both require a lot of simulation.
In the following, the first part of the paper will be dedicated to the outline of the CUR approximation used in this research. This low rank tensor approximation is usually devoted to the reduction of a tensor solely and thus the second part will show how to introduce parametric interpolation capability into this model. After, the numerical computations of the coefficients of the reduced order model will be outlined and some convergence studies will draw some properties of the reduced order model. Finally, first numerical results will exhibit the potential of the reduced order model in industrial applications. These examples are a car rear wheel assembly side hit and a full car crash model.

CUR low rank approximation
Let ( , ) be the space time function representing the simulation results observed on a high-fidelity simulation and for the field . This field is a physical quantity, for instance in crash simulations, it may be displacements or plastic strain. It must be noticed that displacements along the x, y, and z axes are three different fields. There is no particular assumptions concerning this function. The approximation of the high-fidelity field ( , ) made by the reduced order model is denoted . The developed reduced order model is based on the CUR decomposition [34]. It approximates a field of one simulation using separable functions: Where ( ) are the vectors of the time basis, ( ) the vectors of the space basis, and are the projections of the field on these bases. and denote the size of each basis. If all the rows and columns are selected as basis vectors, then the high-fidelity data will be surely perfectly reproduced. Thus adding rows and columns to given reduced bases will decrease or at least will remain unchanged the gap between the approximations and the references. The gap may be measured in several ways depending on the chosen norm. It will be unchanged if the supplementary vectors are linear combinations of vectors already included in the bases, meaning that the bases are not enriched by these vectors. It may append in mechanical simulations because time steps or locations may be partly correlated because some simulated phenomena are linear and causal. On another side, the size and the shape of the matrix gathering the coefficients in the equation (1) will characterize the efficiency of the reduced order model: a lower number of non-zero elements stored in will improve the compression rate. Finally and may be chosen as big as allowed by the memory of the computer if the selection of the relevant rows and columns ensuring a minimal level of error is done during the computation of the coefficients in . It is interesting to notice that the previous decomposition (1) is similar to Tucker decomposition [28].
In the Principal Component Analysis (PCA) [2] these bases are computed as the eigenvectors of the correlation function. Whereas, in the CUR approach, the bases are built selecting their vectors directly in the high-fidelity results. It means that the time basis is composed of snapshots of the entire system at some selected instants and the space basis is composed of the entire history of some selected points . Moreover it is possible to choose the basis vectors in other fields than the field or simulations to complete the bases. Then each field may be explained by the data from other fields. For instance the displacement along the X axis of some nodes may be explained advantageously by the displacement according to Y or Z of other nodes. Each coefficient in is associated with one vector of the time basis and one vector of the space basis. Hence, the value of this coefficient will give the importance of this combination and then it will supply a skeleton of the crash scenario [35]. It is said by Wang and Zhang [31] that another advantage of CUR compared to SVD concerns the preservation of the structural properties of the original data matrices such that sparsity or non-negativity. One drawback of the CUR technique is that the bases are no more orthonormal. Therefore, more vectors are required in order to get the same precision than an orthonormal basis and numerical issues may arise.
The selection of the columns and rows in Eq. (1) may be done using different approaches. Traditionally in the CUR decomposition, rows and columns are selected randomly according to the distribution of a leverage scores. This scores may be based on energetic or correlation considerations for instance [35]. Enrichment may be performed using iterative greedy or Tabu algorithms. Some of them make a link between the DEIM and the CUR factorization [36]. The coefficients in are then computed such that the reduced order model minimizes a given error between the approximation of the matrix and the highfidelity matrix. The iterative approaches may need numerous iterations to converge that may limit the industrial applications. That is the reason, we choose in the outlined approach to select rows and columns representative of the computations with an unsupervised learning approach. The classic k-means clustering method [37] is used here in a sake of simplicity of application and computational efficiency. The rows or the columns closest to the centers of each cluster are selected as a potential vector of the basis. If the numbers of row and column are sufficient enough then the selection will be a good representation of the high-fidelity data such that the approximation will be close to the reference. Moreover the selected rows and columns will be as different as possible according to the norm used in the k-means algorithm ( norm), limiting stability problems.
For illustrative purpose, the results of a clustering of the displacement of the structure displayed on figure 1 are shown on figure 2 and 3. This structure contains 16878 nodes and 76 instants. The high-fidelity simulation was done using Virtual Performance Solution of ESI Group in about 15 minutes. 15 columns and 5 rows are selected for each field. They will be the points closest to the centers of the clusters. They will constitute the vectors of the space and time bases. It is interesting to notice that even there are only 15 columns and 5 rows, they spread on the high-fidelity data quite well. Moreover they represent the variety of features of the whole columns and rows. Despite the well-known limitation of the k-means algorithm for high-dimensional data [38], it leads to satisfactory results. A study of the impact of the number of clusters on the error of the approximation will be led in the part 4.3 on the same test-case. Nevertheless other clustering and unsupervised learning algorithms of pattern recognition should be tried in future works. And an indicator measuring the quality of the clustering should also be set-up in order to infer the relevance of the bases to build the reduced order model.   As well as the selection of rows and columns, there are several ways to compute the coefficients of the low rank approximation. Each of them corresponds to the minimization of a dedicated norm of the error between highfidelity field and approximated field such that: As mentioned previously, it may be wise to couple this minimization problem with a penalization problem. It will indeed minimize the number of non-zero coefficients and select directly the most important and influent coefficients for the minimization of the error. Moreover this property will improve the stability of the reduced order model that will be useful for interpolation, avoiding hence the well-known overfitting phenomenon. The minimization problem becomes then: Where ' is a scalar defining the trade-off between the minimization of the error and the penalization term. As the minimization of the error is more important, it is chosen much lower than 1. The choice of the norm in the problem (2) and the resolution scheme will be outlined in the section 4.
Before, parameters has to be included in equation (1). In this equation, the approximation is indeed done for one simulation solely and there is no parametric dependency in this formulation. In the next part, the parametric reduced order model will be described.

Introduction of the parametric interpolation capability into the CUR model
In order to consider the parametric dependency of the model, functions depending on the parametric field denoted by )( , ) should be introduced in equation (1). Note that the parameter is considered as a space and time varying field. This modelling will allow to consider the parameters more parsimoniously. For instance, all the thicknesses of the model will be considered as a unique space varying parameter. It may be space and time varying but only scalar parametric field are possible here. For a sake of clarity, only the case of a unique parametric field will be studied here but the method is easily generalized to multiple parametric field case. In this case, each parametric field will correspond to one space time parameter, let say the thickness field and the Young modulus field for instance. Then, it is possible to perform an expansion of the high-fidelity results over a set of functions depending on the parameter. It leads to: In equation (3), ) * ( , ) denotes the value of the parameter µ at a location and an instant for the simulation indexed by . . is the number of parametric functions + , . For instance, they can be polynomials or trigonometric functions. The choice depends on the system's behaviour knowledge. In general polynomial expansion is a clever first guess if there is no specific knowledge about the system and the physic. If frequency simulation instead of crash simulation was considered, then trigonometric functions will be surely more appropriate. As for the number of rows and columns in the bases, the number of parameteric functions must be chosen as large as the memory of the computer allows letting the selection be done during the computation of the coefficients in . Compared to Eq. (1), a supplementary dimension has been added to the tensor gathering the coefficients of the reduced order model. Because, one coefficient corresponds to one row , one column , and one parametric function + , . This formulation are typical Hierarchical Tucker tensor [39]. It may also be identified to TT-cross decomposition [30]. If several parameters are considered, then it will only increase the order of the tensor .
Nevertheless a limitation of the previous reduced order model is that the approximation at a location at an instant depends only on the value of the parameter at the same location and the same instant. In some cases, it may be very restrictive. In crash simulation, the displacement of a part may be influenced by the thickness at another location. This effect may be considered through the introduction of "interactions" into the space and time bases such that: ( , , ) * ) = / ,0 ( , ) * ( , )) 1 , 2 , ) * ( , )3+ , () * ( , )) ,0, , ,, Here, / 2 , )( , )3 is the time basis vector that now depends on the value of the parameter at the location , that is the location corresponding to the basis vector . Then, using this kind of interaction, the approximation at a location and an instant may depend on the value at other locations but at the same instant. The 6 parametric functions may be polynomials or other functions according to a priori knowledge such that: where ℎ 0 is a function depending on the parametric field supplied by the user. Accordingly, the parameterized space basis 1 may be decomposed as: where : is a function depending on the parametric field ) * supplied by the user. ; is the number of time interaction functions. Using this kind of interaction will allow to explain the high-fidelity field at a location and an instant using the value of the parametric field at the same location but a different instant. Obviously if the parametric field is invariant according to the time, then it is useless to consider such dependency. It must be noticed that the basis vectors and can be chosen in any fields or simulations. It is recalled that they are selected using a clustering approach. In order to take into account the interactions efficiently, it may be useful to impose some rows (instants) and columns (nodes) as basis vectors if the parametric field is varying at these nodes and instants between simulations. Automatic selection of such locations and instants may be done by a dedicated clustering approach and will be the subject of ongoing research. Supplementary dimensions had been added to the coefficients tensor according to these new parametric dependencies.
The previous formulation (4) is continuous in space and time. Nevertheless, the finite element simulation results are discretized over a 3D mesh with < nodes and a = temporal grid and are stored in = × < matrices. The approximation at a given location now indexed by and instant now indexed by of a highfidelity matrix ? *, ∈ ℝ B × C is then: We denote by F G * ∈ ℝ B × × 5 the 3 rd order tensor gathering the column vectors F G ,0 * ∈ ℝ B and I G * ∈ ℝ × C × 4 the 3 rd order tensor gathering the row vectors I G , * ∈ ℝ C . Then, / = 0 * , (resp 1 < * ) denotes the i-th (resp j-th) basis vector at the instant (resp location) indexed by (resp ) for the simulation and the parametric function indexed by J (resp n). Similarly, E =<, * is the value of the K-th parametric function at location indexed by and instant indexed by . All these values of the parameter are gathered into an order 3 tensor L * ∈ ℝ B × C ×where the first dimension is the time, the second dimension is the space and the last dimension is the parametric function. If the equation (5) is reformulated as: , where the dot stands for the usual inner product, then, the modified CUR method clearly appears as a multilinear regression of unknown O and coefficients taken in the vector M =< * = F G * P = ⊗ I G * P < R ⊗ L * , where P is the ith element of the standard basis. M * is a five order tensor extracted from the M * tensor that corresponds to the location indexed by and the instant indexed by .
The coefficients of the regression are the same for all the instants and locations of the high-fidelity computations. The regression aspect of the reduced order model obvious in equation (6) leads to the name ReCUR pointing for Regression CUR. The O are computed using the existing high-fidelity simulations. The latter may be considered as the training set of the reduced order model. When considering several computations, the minimization problem (2) becomes, in a matrix form: In equation (7), the matrices ? and ? are the concatenations of the matrices corresponding to all the computations (respectively high-fidelity and approximated) for the field . The minimization is thus done considering all the computations, it means that the norm is calculated over all the nodes, all the instants and all the simulations for each field. Once the coefficients are calculated, the reduced order model is ready to be used.
An estimation for a new parametric field is done simply by building the tensor M using any desired parameter's values ). The reduced order model may be considered as a function of the parameter such that: Although this computation may be done very efficiently, it must be noticed that it can be faster using coarser discretisation grids for the approximation than the ones used for the computation of the unknown vector O . Typically, if the mesh is very fine only to ensure the convergence of the high-fidelity computation, then it may be useless to use all the points of the mesh for the computation of the reduced order model. Similarly it is possible to use the reduced order model to compute only some very specific nodes and instants, for instance those that are useful in the optimization problem. In the next part, a focus on the numerical computation of the coefficients of the reduced order model is done.

Numerical computation of the coefficients of the reduced order model
The core of the method is to compute the coefficients of the regression O solving the optimization problem (7). In this part, the computation method of the coefficients is outlined.

Formulation of the linear programs
The coefficients are computed such that the reduced order model minimizes an error between the approximation done by the reduced order model and the high-fidelity computation. It is recalled that, in order to avoid stability and overfitting problems, it is necessary to have a selection technique of the most influential coefficients. It is done by the use of the penalization term in (7).
It has been chosen to minimize the infinite norm of the error together with the norm of the coefficients. It has been chosen for a sake of performance. The minimization problem may indeed be formulated as a linear and purely continuous problem with a limited number of variable. The first step is to state the optimization problem: Here `a = . 6 ; corresponds to the number of coefficients in O ; S 0T< is the infinite error for the field k and for all the considered computations.
The first term in the objective function stands for the minimization of the error field by field while the second term stands for the minimization of the number of non-zero coefficients. It is very close to the Dantzig selector [40]. This problem may be written easily as a linear problem by decomposing in a positive and a negative part such that = b ,c − b ,d . The linear program is then: It is now possible to benefit from the efficiency of linear programming. Here the number of variables is equal to 1 + 2 `a for each field , and the number of constraints is 2 * < = * . The number of variables is much lower than the number of constraints because the compression rate `a is directly proportional to the number of rows and columns selected in the clustering step.
Nevertheless computational issues may arise in the resolution of the problem (9). The number of constraints is indeed huge enough such that the memory footprint may be very high. To circumvent this issue, a constraint generation process has been implemented. The algorithm is the following: first some constraints in the entire set of constraints are selected by a dedicated process, randomly, uniformly or by a clustering on the values in ? ,* or on the rows of M * . This set of constraints is called g . It is then possible to solve the following linear problem: The number of constraints in this problem is now 2 times the length of g . After solving this linear problem, the violations for every constraints may be evaluated using (8). It corresponds to the gap between the high-fidelity data D BC ,* and the approximations made by the current reduced order model M =< * • O . In the next iteration, the worst predicted values for each field are added to the linear program. Few iterations of this process greatly increase the quality of the approximation decreasing the infinite error. In this approach, the data outside g may be considered as a validation set of the current reduced order model.

Additional constraints
Moreover, an advantage of this formulation is its flexibility. It is indeed possible to introduce other constraints in the linear problem in order to guide the construction of the reduced order model. The only limitation, is that these constraints should be linear according to the coefficients of the reduced order model. It may be possible to integrate specific behaviour into the reduced order model using such constraints.
For instance, one could integrate welding behaviour into the reduced order model. The additional constraints will then be for each displacement field , , ∀ ∈ , l ∈ l , ∈ Z1, = \, ∈ Z1, * \ Where j is a small margin in order to guarantee the solvability of the problem and and ′ are two set of nodes (typically the borders of each welded part). These set of constraints will ensure that the nodes in the set and ′ will never be further than the distance j in each direction for the simulations in the training set. Note that this behaviour will not be ensured for an interpolated simulation that is not in the training set but it will help the reduced order model to be more relevant according to this behaviour.
Such kind of linear constraint may also force the reduced order model to reproduce exactly the quantity of interest of some nodes or instants known as important. The additional constraints will then be for each field involved in the quantity of interest: −D ,* =< − j ≤ M =< * • O ≤ D ,* =< + j , ∀ ∈ , ∈ n, ∈ Z1, * \ Where and n are the sets of the interesting nodes and instants. These constraints may in fact be considered as the inclusion of expert judgements into the construction of the reduced order model. It must be noticed that if the constraints is not directly linear, then a variable change or introduction of a new variable may solve this shortcoming, even if it will add a supplementary parameter in the reduced order model.
A last kind of additional constraints may be helpful for the solving of the optimisation problem (10). It is indeed possible to impose a constraint on the maximal error S 0T< . If so, the coefficients O will ensure this error for the selected constraints in g and the optimization problem will only minimize the number of non-zero coefficients. It may improve the quality of the reduced order model, it has been indeed noticed that a model with a higher error S 0T< on the selected constraints in g has a lower maximal error for the entire set of constraints. In a future study, an iterative process will be set-up in order to find the optimal value of this impose error while limiting the number of constraint generation iteration.
In order to have a synthetic view of the workflow, the figure 4 sums up the entire process of construction of a reduced order model. The inputs of the method are the data coming from the existing simulations and the corresponding parameter fields while the output is the reduced order model ready to be used for interpolation. The data transformations step in the process stands for all the pre-treatments applied to the raw data such as normalization, noise smoothing, or sampling for instance. These data transformations may improve greatly the quality of the interpolation as well as the computation time of the reduced order model. Nevertheless it must be kept in mind that, on account of these transformations, corresponding inverse transformations should be applied to the interpolation supplied by the reduced order model.

First results of parametric studies
To better understand the behavior of the reduced order model according to its inner parameters such that the numbers of clusters or the number of computations in the training set, some convergence studies have been led. The used test-case is the one presented in the section 2. There are 29 computations available on which the quality of the interpolation of the reduced order model will be evaluated. The first 25 simulations have been chosen according to a Latin Hypercube Sampling design of experiments. The last 4 configurations correspond to 4 corners of the design space. It is remind that there are 16878 nodes and 76 instants in the finite element model. For information, the highfidelity computations take about 15 minutes using the fast dynamics software VPS of ESI Group. The studied parameters are the 9 thicknesses of the deformable parts. The following reduced order models are built using Legendre polynomials of order 2 for the parametric dependency. Such polynomials are very interesting according to their advantageous properties for regression models. The reduction had been done on the displacement according to X, Y and Z displacements. Each displacement field has been normalized according to the maximum displacement of all the computations in the training set. By avoiding scale problems between the centers of the clusters of all the fields, it will help the resolution of the linear program. It must be hence kept in mind that the bases will be constituted by vectors extracted by clustering from all the fields (it might be displacements and stresses for instance), and hence may have not at all the same scale. The same numbers of row and column clusters are chosen in each kind of field. The clustering is performed on each field separately but considering all the computations together. The rigid bodies have not be considered in the construction of the reduced order model as well as in the interpolation.
To study the influence of the numbers of clusters, reduced order models using 2 computations were considered. They have the following parameters:   The error between the high-fidelity computations and the interpolations of the reduced order model is expressed through the root mean square error of the magnitude of the displacements at every points and instants. It is normalized by the norm of the high-fidelity computation.
For the convergence study on the size of the time basis, that is the number of column clusters, the number of constraints taken into in the problem (10)  It is interesting to note that the computations in the training set are the ones with the lowest errors. Besides, the variance of the reduced order model seems to be quite steady. In the following the value of 25 column clusters times the number of computations in the training set will be adopted for this test-case. It represents indeed a good tradeoff between the computation time, the bias, and the variance. For the study of the impact of the size of the space basis, that corresponds to the number of the clusters in the row clustering, reduced order models with 396 constraints over which a maximal error of 0.21 for constraints regarding X axis are allowed, 0.11 regarding the Y axis, and 0.08 according to Z axis.
After a great fall, the mean values in the figure 6 seems to be quite steady for > 4. The same behavior is observed for the worst predicted simulation for each reduced order model. Thus, despite the randomness of the k-means algorithm, it shows that selecting a sufficient number of clusters will lead to a robust reduced order model in term of forecast. From > 10, the reduced order models look very similar. It is interesting to note that the model with 2 row clusters has a very low variance but high bias and adding members to the basis is increasing this variance but decrease the bias. For the following, the value = 5 row clusters times the number of computations will be adopted. After, an analysis of the quality of the interpolation capability of the reduced order model according to the number of computations in the training set has been led. In order to refine the model, the worst predicted computation is added in the training set. It leads to the behavior displayed on the figure 7. The figure 7 shows that the means of the errors decrease as computations are added to the model. The standard deviation decreases also. It is interesting to denote that even if the linear program minimizes the infinite error, the behaviour of the error is also satisfying. The refinement technique used here give satisfactory results, but other techniques could be investigate. For instance if the reduced order model is used in an optimisation study, then the added computations could be decided studying the errors on the quantities of interest in the optimisation study. This will be the strategy used in the next section.
For information, the figure 8 shows the deformed body of an interpolation using the reduced order model with 8 computations in the training set. It displays the geometry of the final instant of the worst predicted configuration. The interpolation and the computation of the output files take less than 1 minute. There are 111 coefficients according to X, 59 coefficients according to Y, and 85 according to the Z direction that is very limited compared to the potential number of coefficients in the linear program. The left part of the device is very well predicted while the error is higher on the right part. It could be possible to add relevant computations or constraints in the construction of the reduced order model in order to refine the model in this area. Moreover it shows that perhaps the error is not always the best choice to select computations to refine the model, and perhaps choosing to add the computation with the maximal . error would lead to a refinement in this area. Comparison to existing reduction methods should be very interesting in a future work. Nevertheless there is a lack of such a method that are dedicated to many parameters computations, especially in the field of crash simulation. Comparisons should be done according to the level of error and to the number of computations required to build the reduced order model efficient for its dedicated use. It could be very interesting to have a comparison to recent works coupling POD and neural network [41], or POD and Radial Basis Function [42]. Nevertheless applications of these studies lie in computation fluid dynamics. But the method outlined in the present paper could be easily adapted to any kind of physics according to the full non intrusiveness of the approach.
These preliminary results show relevant features and exhibit the interest of the method. Nevertheless other convergence or sensitivity studies according to the parameters of the method should be led in future works, for instance to have the sensitivity according to the order of the polynomials or other refinement strategies, or the effect of the randomness in the k-means and in the initialization of the linear program. It is indeed necessary to test the robustness of the model. To limit these impacts, the number of iteration of constraint generation and the numbers of clusters are set to a quite high value in the present models. Hence the most penalizing constraints will be most surely selected in the linear program.
The next part will sum-up the results of the method on industrial optimization problems.

Numerical example
The developed method was applied to industrial test-cases using code developed in Julia 0.4. Two examples are discussed here: a rear wheel assembly and a car structure models. In all examples, high-fidelity computations were ran using the VPS solver from ESI group.

Rear wheel assembly
The rear wheel assembly model consists of the left wheel, fixations to the chassis and crossbar. It simulates a hit by a pendulum on the side of the left wheel.
The design objective in this industrial case is to limit the localisation of structural damage to the crossbar while limiting the overall weight of the assembly. Three parameters are considered: thicknesses of the yoke holders, arm and crossbar (cf. figure 9 and table 2). The assembly is considered as fixed to the chassis, as such the simulation have four fixed points located at the two fixtures and on the end of the two choc absorbers.  The model was constituted of 255,675 nodes of which 41,628 were used for the creation of the reduction after removal of the impactor, non-shell parts and non-structural parts. This represents a total of 31 parts and 39,704 shells. 41 instants representing a total simulated time of 200 ms were kept.
A reduction model was constructed using 5 experiences (see table 3). 4 fields were reduced: local displacements (along X, Y and Z) and plastic strain. Shell thickness was used as a parameter field. All shell defined fields (plastic strain and thickness) were projected on the nodes using standard average when a node belongs to several parts. To avoid order of magnitude inconsistency, all the fields were normalized according to the absolute maximum of each field for all the simulations. A total of 250 clusters on the columns (nodes) and 15 clusters for the lines (instants) were extracted using k-means method. These elements constituted the and vectors. 0.005% of the total constraints were kept and additional constraints on all field at first and last instants for 4 points along the crossbar were added. The additional constraints were chosen to improve prediction of crossbar bending. They may be considered as introduction of expert advice in the reduced order model. After a first resolution, a second iteration was made adding worst predictions as constraints. This resulted in the use of 4080 constraints. The linear problem was solved using CPLEX. Imposed εmax was 0.05. The parametric functions are third order Legendre polynomials, normed between ±10% of minimum and maximum of the thickness field. The advantage of these polynomials is that they are orthonormal between [-1,1]. This improves the efficiency of the regression model. The margin allows to benefit of this property also for extrapolation use of the reduced order model.  The reduced order model construction ran in 315 min using slightly above 500 GB of RAM. In the 240,000 possible coefficients for each field, only 52, 50, 64 and 21 are non-null for the X,Y,Z displacements and plastic strain fields respectively. Using this model, generation of a new estimations takes 39s on average using one processor. For sake of comparison calculation of a new high-fidelity experience took 24 hours on 8 processors. It is important to note that this model is in early development stage. As such memory consumption and computation time can be largely improved.
The general quality of the reduction model was evaluated using a sixth case that had been calculated using VPS for comparison. This case was chosen as it was not close to any experiences used in reduction's generation. Criteria used to assess model's quality were both objectives (mean error, localisation of the plastic strain maximum) and subjective (general deformation in accordance between reduced and high-fidelity models, localisation of areas of high plastic strain). The rear wheel assembly seen from the top for the final instant is presented in figure 10. At first glance it appears that model deformation is well predicted. The bend observed on the left side of the crossbeam (red arrows) presents the same shape and localisation in both cases. This observation is confirmed by the displacement error which is limited to less than 15% on the crossbeam and 25% overall.
In figure 11 the plastic strain nodal projection at the final instant on highfidelity (right) and reduced order model (left) are presented. To be usable for Impact Impact industrial developments two aspects need to be respected in this field: localisation of areas of high plastic strain (which would indicate plastic deformation and thus when damage occurs to parts) and localisation of the maximum of the plastic strain field (which shows the localisation of the concentration of forces during the impact). The second criteria is well respected in this case with a distance of 5.3 mm between high-fidelity and reduction. Localisations of areas of high plastic strain (1-5 on figure 11) is generally well predicted. Area 1 and 2 are correctly predicted both on localisations and scales. Area 3 appears to be over-predicted with a larger area and higher values. Area 4 and 5 are predicted while not existing on the highfidelity model. However those areas have a significantly lower value than the main deformation zone (area 1), 0.04 versus 0.09. The reduced order model also predicted negative values for the plastic strain (minimum -0.008) due to interpolation, while those are not physical they can be assimilated to 0 indicating the absence of plastic deformation at those points. Where < is the total number of nodes and ? ./ is the field being evaluated by reduction or high-fidelity models. Using this error for the plastic strain field allows to take into account the error on the location on the maximum value more than its predicted value. The location of the maximum value is in this case the value of interest.  It appears that the model prediction is sufficiently accurate to provide a valuable and usable surrogate to the high-fidelity counterpart during the exploratory phase of design (table 4). The low cost associated with the estimation after the creation of the model (39s for the evaluation of a full model) introduce the possibility for the design team to evaluate considerably more options than with a high-fidelity approach. Moreover, it is possible to use the reduced order model to estimate only few nodal or time data. Estimation time is then much lower. It may be useful in an optimization study. Specifications are evaluated using the reduced order model only for the nodes and instants useful for the computations of these quantities. This allows the use of the reduced order model in an iterative process, at each iteration an approximate reduction model is used to evaluate large amounts of design options with relative accuracy. From these, some good candidates can be selected and evaluated using high-fidelity models. Then those candidates can be used as reference for a new reduction more tailored to the target zone. This approach is expected to considerably reduce the overall number of high-fidelity calculations required, leading to more innovative design phase. The second example consists of a proof of concept of this approach.

Car model
This model consists of a full car structure impacting a barrier type obstacle on the front left side at a speed of 70 km/h. The full vehicle model comprises 2848 parts for a total of 4,250,821 nodes. After removal of non-shell parts and barrier, 1481 parts representing 2,835,289 nodes were left. 19 instants were taken for a total duration of 90 ms.
The goal in this test-case was to use the method with resource constraints to conduct a prototype of carry-over (CO) study. It was conducted using 37 variables corresponding to the thicknesses of 51 parts (including 14 symmetrical ones) representing a total mass of 77.8 kg. The parameters varying between 0.6 and 1.4 mm. The aim being to reduce the overall mass of the vehicle while keeping safety parameters in the same range as the existing vehicle.
The existing computations were the nominal model (commercialized car) and two selected points (1 and 2, table 5) presenting reduction in mass of 5 and 7.3 kg. Those points were chosen based on mass reduction compared to the nominal case and variables spread which would increase learning possibilities for the reduction. The resulting reduction model was then used to calculate specific objectives values on 10,000 experiences (generated in the parametrized domain). Safety objectives were the shortening of the left front girder (compression) and maximum displacement of the pedal floor (intrusion). Profiting from the possibility to evaluate the reduction model at specific nodes and times, calculation of the 10,000 experiences values took 5 hours on 85 processors. From those, two "best candidates" were chosen using both results evaluations and parameters spreads between them (again to maximize learning possibilities). Those candidates were then calculated in high-fidelity which took 15 hours per experiment on 88 processors, compared to 17 min for the evaluation of an entire vehicle using the reduction model on a single processor. Using those two experiences and one from the original model a new reduction was created. Then evaluation of the 10,000 points was repeated leading to the identification of 3 candidates points that, according to the reduction model, would fulfil safety criteria while minimizing mass and carry-over rate. Those three candidates were then calculated in high-fidelity and a final "winner" configuration identified (experience 7, This procedure led in 7 calculations to a configuration with similar safety measures (see table 5) than the nominal case but with a reduction of 4.9kg in mass. While this is a proof of concept study and obviously more constraints are considered in car design than those related to front crash, this clearly shows the potential and efficacy of this reduction method. Using a classical optimisation method would require around 10 high-fidelity calculations per variable this would have amounted to 370 experiences or 488,400 processor hours. Using ReCUR for this problem limited this requirement to 9,240 processor hours for high-fidelity calculations, 25 for the model construction and 860 for the model evaluations, or a total of 10,125 processor hours. This represents 2.1% of a standard optimisation study.
The first model was made using two experiences (1 and 2 in table 5). As such first order Legendre polynomials were used. The second model used three experiences (2, 3 and 4 in table 5) and second order Legendre polynomials were used. 3 fields were reduced: local displacements along X, Y and Z. Shell thickness was used as parameter field. Thicknesses were projected on the nodes. 200 clusters along the columns (nodes) and 20 clusters along the lines (instants) per fields were chosen using k-means method. A proportion of 5E-6 of the total constraints were kept. In total 3060 constraints were used in the second model. Also, hard constraints on 7 points of interest (along known structural position or used for the calculation of objectives functions) were imposed for each field at initial and final instants. The imposed εmax was 0.05. The data were normalized as in the previous test-case.
The second reduction (ROM-2) ran in 12 hours using 300 GB of RAM. The regression found 21, 38 and 32 non-null coefficients for the X,Y and Z respectively among the 291,600 possible coefficients. It is interesting to notice that the all the non-zero coefficients correspond to order 2 polynomials for the column interaction.
The second reduction model appeared more precise than the first (ROM-1) in predicting objective values on points of interest (10.3 vs 18.5 % error). This appears logical considering that with an increased number of experiences to learn, ROM-2 offers a broader coverage of the range of parameters.
The case presented here corresponds to the final "best" configuration found during the study. It was evaluated using ROM-2. In figure 12 the vehicle front structural parts deformation at final instant are presented. The general shape of the deformation appears well predicted by the reduction model. The maximum normalized difference between the two models is below 26% with most of the structural parts having a deformation under 13% ( figure 12). Considering the size of the model, large amount of variable involved and early stage of developments for the reduction model, this is a promising result. On figure 13 the projection of the displacement error between the VPS and reduced order model for the parts used as parameters in the study are shown.
For those parts the error on the norm of the displacement is limited to 10%. This indicates that the model can be used with relatively good accuracy to predict specifications. Errors on the predicted specifications were less than 18% with an average of 10.5%. This indicates that, while this study was made as a proof of concept, using a non-optimized early development model, results could be used in industrial development. The model still suffers from infancy issues, such as a large memory footprint rending the use of large amounts of high-fidelity cases for the generation of reduction impractical. Another issue is the long computation time required to generate the model (12 hours in this case), while this still reduced tremendously the cost of a study compared to the use of high-fidelity model this prevents the use of the model in a truly interactive setting. Finally, in the current version equal weights are given to all constraints in the model, meaning that most of the power of the model is spent in predicting the rigid body motion rather than relevant deformation. This will be solved in the future by filtering rigid body motion before the creation of the model.

Conclusion
In this article, construction of a reduced order model called ReCUR has been outlined as well as its use to estimate new technical solutions. It consists in a regression model coupled with a low rank tensor approximation method. The two application examples show the potential of the method in the context of the fast dynamic simulations. It allows, indeed, to make evaluation in a few minutes for industrial test-cases, with still large potential for code optimization and parallelization to improve efficiency. The domain of application is crash simulation, however the method is not limited to this physic. The method may be useful as soon as time consuming high-fidelity computations are required in a design process. Potential domain can then be combustion simulations or fluid mechanics among others. These preliminary results look very promising for use in optimization study knowing that the ReCUR reduce order model is newly developed and a lot of improvement tracks exist. The number of simulations required to improve the existing system of a full car model is indeed much lower than for the current approach relying on the design of experiment and surface responses.
One asset of the method is the use of the data coming from all the existing computations while keeping the dimension of the problem tractable. It might allow industrials to shorten and improve their conception stages as well as having a better data employment. Another advantage of the method is its flexibility thanks the possibility of introducing expert's judgments directly into the construction of the reduced order model. This possibility enhances greatly the capability of the reduced order model.
Among all the clues of improvement spread in the text, a first lead would be improvement of the resolution by the linear programming. Algorithms may be set up to avoid fixing the value of the desired error. It may find the one minimizing the error for constraints not taken in the linear problem, using a Lagrangian for instance or another iterative process. Moreover, adding a computation in the training set of the reduced order model lead with the actual process to the construction of a new reduced order model from scratch. It may be very useful to start from the previous reduced order model with a lower number of computations and to directly refined it with extra clusters and constraints in the linear programming.
The use of transformed data in the construction of the reduced order model must be also investigated. It may allow to focus on the most relevant phenomena featured in the high-fidelity simulation. An example may be filtering nodes having no deformation and replacing them with a rigid body motion or removing the temporal motions induced by wave propagation. This method may be inspired by the technique used in APHR or DEIM.
Works must be done concerning the estimation of quality of the reduced order models. A systematic pertinent error indicator must be built to compare estimations of reduced order model to high-fidelity data. This point is not straightforward because there are a lot of ways to characterize the relevance of a reduced order model. It may involve a global error as the or error and some particular errors on the values of interest or more likely a combination of the both. Having such error may allow to have a relevant and systematic refinement strategies. Moreover, a benchmark to existing approaches such as PGD, or DEIM may be useful to compare the method according to the state of art, even if there is a lack of such method for many parameters and nonintrusive case. Some recent works combining the POD and machine learning techniques such as neural network [41] or radial basis functions [42] exist and may be use for comparison. Moreover the results from an optimization study using the reduced order model from ReCUR may be compared to alternative approach lying on the use of design and experiment and response surfaces [27].
Last, it lacks a tool to assist in the analysis of rows and columns selected by the clustering as well as the value of the corresponding non-zero coefficients. Such analysis will be indeed very useful to understand which nodes and parts are influential and which instants are critical during crashes. Analysis of the results given by the reduced order model may be aided by the derivation of sensitivity indices. It may help to quantify the influence of the parameters on some objective functions or specifications of the requirements. Derivation of local sensitivity indices may be directly done using the gradient of the reduced order model according to the parameters. Global sensitivity indices such as Sobol indices [43] require more advance computations but may still be done using the reduced order model. As the reduced order model has an analytical form, these derivations may be done very efficiently.
Nevertheless, ReCUR could already be used in several kinds of processes. For instance, the link with the optimization study could be investigated more deeply, especially for optimization with discrete parameters [44]. ReCUR Jean-Patrick is a research engineer at the Technological Research Institute SystemX. Holding two masters in mechanical and petroleum engineering from the Ecole Centrale Nantes and Penn State University he worked for two years on the project ROM (Multiphysic's reduction and optimization) in the IRT SystemX. His research interests are high performance computing and collaborative engineering solutions.
Fatima-Zohra Daim: Fatima-Zohra Daim (PhD in Applied Mathematics) is a research engineer working in ESI Group, a leading simulation software company. She is part of the Technology Center Of Excellence whose role is to pursue innovative multidisciplinary research projects namely in fluid mechanics, structural mechanics, optimization and numerical algorithms. These past few years, her most active research domain focuses on Reduced Order Modeling in Crash simulation.
Ming Chau: Ming Chau is an engineer and researcher in Applied Mathematics and Computer Science. He currently works at IRT SystemX. He received his PhD in 2005 from INP Toulouse, France. His research interests include sparse linear solvers, image processing and numerical simulation, with applications to mechanical engineering and medical imaging.
Yves Tourbier: Yves Tourbier, PhD, assumes the role of technical expert for the field "optimization and decision making" within RENAULT engineering. He has led research projects using numerical simulation and optimization applied to the improvement of mass / performance compromises of vehicle structure (BIW) and engines. For four years Yves Tourbier led the project "Model Reduction & Multiphysics Optimization" at IRT SystemX. His current research focuses on model reduction, deep learning and autonomous driving vehicle validation.