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

Yves Le Guennec, Jean-Patrick Brunet, Fatima Zohra Daim, Ming Chau, Yves Tourbier

► To cite this version:

Yves Le Guennec, Jean-Patrick Brunet, Fatima Zohra Daim, Ming Chau, Yves Tourbier. A parametric and non-intrusive reduced order model of car crash simulation. Computer Methods in Applied Mechanics and Engineering, 2018, 338. hal-01485276

HAL Id: hal-01485276

https://hal.science/hal-01485276v3

Submitted on 2 Mar 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A parametric and non-intrusive reduced order model of car crash simulation

Y. Le Guennec*ª (yves.leguennec@irt-systemx.fr), J.-P. Brunetª (jean-patrick.brunet@irt-systemx.fr), F.-Z. Daimᵇ (fda@esi-group.com), M. Chauª (ming.chau@irt-systemx.fr), Y. Tourbierᶜ (yves.tourbier@renault.com)

a : IRT SystemX, 8 avenue de la Vauve, 91120 Paris-Saclay, France;

b : ESI Group, 99 rue des Solets, 94513 Rungis, France ;

c : Renault, 1 avenue du Golf, 78084 Guyancourt, France.

## 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.

<footnote>* Corresponding author</footnote>
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 highlights the promising power of
the method to shorten design process using optimisation and long-run
simulations.

Keywords: Reduced order model, Crash simulation, regression analysis, linear programming.

# 1. 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–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–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 post-
processings 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–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.

## 2. CUR low rank approximation

Let $a_{HF}^k(t, x)$ be the space time function representing the simulation results observed on a high-fidelity simulation and for the field $k$. 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 $a_{HF}^k(t, x)$ made by the reduced order model is denoted $a^k$. The developed reduced order model is based on the CUR decomposition [34]. It approximates a field of one simulation using separable functions:
$$a_{HF}^{k}(t, x) \approx a^{k}(t, x) = \sum_{i=1}^{n_c} \sum_{j=1}^{n_r} C_{i}(t) U_{ij}^{k} R_{j}(x) \quad (1)$$

Where $C_i(t)$ are the vectors of the time basis, $R_j(x)$ the vectors of the space basis, and $U_{ij}^k$ are the projections of the field $k$ on these bases. $n_c$ and $n_r$ 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 $U^k$ gathering the coefficients $U_{ij}^k$ in the equation (1) will characterize the efficiency of the reduced order model: a lower number of non-zero elements stored in $U^k$ will improve the compression rate. Finally $n_c$ and $n_r$ 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 $U^k$. 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 $t_i$ and the space basis is composed of the entire history of some selected points $x_i$. Moreover it is possible to choose the basis vectors in other fields than the field $k$ or simulations $s$ 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 $U^k$ 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 $U^k$ are then computed such that the reduced order model minimizes a given error between the approximation of the matrix and the high-fidelity 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 ($L_2$ 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.
<figure><img src="image_2.png" /><figcaption>Figure 1: View of the impacted structure: It is constituted by 9 deformable bodies and 2 rigid bodies at the top left (in brown, number 16) and at the bottom (in dark green, number 12). The deformable bodies are constituted by 16328 nodes. Left side shows the structure not deformed and right structure shows the structure at the end of simulation, which is after 60 ms.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Clustering of the columns: The grey lines denote the trajectories of all the nodes, the nodes selected by the clustering being coloured. Abscissa is the time indices and ordinates are the displacements normalized according to the maximum value. Left frame shows the x displacement, middle frame shows the y displacement, and right frame shows the z displacement.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Clustering of the rows: The grey lines denote the spatial values for all the instants, the instants selected by the clustering being coloured. Abscissa are the nodes indices and ordinates are the normalized displacements. Left frame shows the x displacement, middle frame shows the y displacement, and right frame shows the z displacement.</figcaption></figure>

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 high-fidelity field and approximated field such that:

$$ \min_{U^k} \|a_{HF}^k - a^k\| $$

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:

$$ \min_{U^k} \|a_{HF}^k - a^k\| + \lambda \|U^k\| \quad (2) $$

Where $\lambda$ 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.

### 3. 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 $\mu(t, x)$ 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:

$$a^k(t, x, \mu^s) = \sum_{i=1}^{n_c} \sum_{j=1}^{n_r} \sum_{l=1}^{n_f} C_i(t) R_j(x) f_l(\mu^s(t, x)) U_{ijl}^k \quad (3)$$

In equation (3), $\mu^s(t, x)$ denotes the value of the parameter $\mu$ at a location $x$ and an instant $t$ for the simulation indexed by $s$. $n_f$ is the number of parametric functions $f_l$. 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 $U^k$. 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 $R_j$, one column $C_i$, and one
parametric function $f_l$. 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 $U^k$.

Nevertheless a limitation of the previous reduced order model is that the approximation at a location $x$ at an instant $t$ 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:

$$ a^k(t, x, \mu^s) = \sum_{i,m=1}^{n_c, n_h} \sum_{j,n=1}^{n_r, n_q} \sum_{l=1}^{n_f} \tilde{C}_{i,m}(t, \mu^s(t, x_i)) \tilde{R}_{j,n}(x, \mu^s(t_j, x)) f_l(\mu^s(t, x)) U_{i,m,j,n,l}^k \quad (4) $$

Here, $\tilde{C}_i(t, \mu(t, x_i))$ is the time basis vector that now depends on the value of the parameter at the location $x_i$, that is the location corresponding to the basis vector $C_i$. 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 $n_h$ parametric functions may be polynomials or other functions according to a priori knowledge such that:

$$ \tilde{C}_{i,m}(t, \mu^s(t, x_i)) = C_i(t)h_m(\mu^s(t, x_i)) $$

where $h_m$ is a function depending on the parametric field supplied by the user. Accordingly, the parameterized space basis $\tilde{R}_j$ may be decomposed as:

$$ \tilde{R}_{j,n}(x, \mu^s(t_j, x)) = R_j(x)q_n(\mu^s(t_j, x)) $$

where $q_n$ is a function depending on the parametric field $\mu^s$ supplied by the user. $n_q$ 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 $C_i$ and $R_j$ 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 $n_x$ nodes and a $n_t$ temporal grid and are stored in $n_t \times n_x$ matrices. The approximation at a given location now indexed by $x$ and instant now indexed by $t$ of a high-fidelity matrix $A_{HF}^{s,k} \in \mathbb{R}^{n_t \times n_x}$ is then:

$$A_{tx}^{s,k} = \sum_{i,m=1}^{n_c, n_h} \sum_{j,n=1}^{n_r, n_q} \sum_{l=1}^{n_f} \tilde{C}_{tim}^{s} \tilde{R}_{jxn}^{s} F_{txl}^{s} U_{imjnl}^{k} \quad (5)$$

We denote by $\tilde{C}^s \in \mathbb{R}^{n_t \times n_c \times n_h}$ the 3rd order tensor gathering the column vectors $\tilde{C}_{im}^s \in \mathbb{R}^{n_t}$ and $\tilde{R}^s \in \mathbb{R}^{n_r \times n_x \times n_q}$ the 3rd order tensor gathering the row vectors $\tilde{R}_{j,n}^s \in \mathbb{R}^{n_x}$. Then, $\tilde{C}_{tim}^s$, (resp $\tilde{R}_{jxn}^s$) denotes the i-th (resp j-th) basis vector at the instant (resp location) indexed by $t$ (resp $x$) for the simulation $s$ and the parametric function indexed by $m$ (resp $n$). Similarly, $F_{txl}^s$ is the value of the $l$-th parametric function at location indexed by $x$ and instant indexed by $t$. All these values of the parameter are gathered into an order 3 tensor $F^s \in \mathbb{R}^{n_t \times n_x \times n_f}$ 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:

$$A_{tx}^{s,k} = X_{tx}^{s} \cdot U^{k}, \quad (6)$$

, where the dot stands for the usual inner product, then, the modified CUR method clearly appears as a multilinear regression of unknown $U^k$ and coefficients taken in the vector $X_{tx}^s = \tilde{C}^s e_t \otimes \tilde{R}^s e_x^T \otimes F^s$, where $e_i$ is the i-th element of the standard basis. $X_{ij}^s$ is a five order tensor extracted from the $X^s$ tensor that corresponds to the location indexed by $x$ and the instant indexed by $t$.

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 $U^k$ 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:
$$ \min_{\mathbf{U}^k} \left| \left| \mathbf{A}_{HF}^k - \mathbf{A}^k \right| \right| + \lambda \| \mathbf{U}^k \| \quad (7) $$

In equation (7), the matrices $\mathbf{A}_{HF}^k$ and $\mathbf{A}^k$ are the concatenations of the matrices corresponding to all the computations (respectively high-fidelity and approximated) for the field $k$. 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 $\mathbf{X}$ using any desired parameter's values $\mu$. The reduced order model may be considered as a function of the parameter such that:

$$ A_{tx}^k(\mu) = \mathbf{X}_{tx}(\mu) \cdot \mathbf{U}^k \quad (8) $$

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 $\mathbf{U}^k$. 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.

## 4. Numerical computation of the coefficients of the reduced order model

The core of the method is to compute the coefficients of the regression $\mathbf{U}^k$ solving the optimization problem (7). In this part, the computation method of the coefficients is outlined.

### 4.1 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 $L_1$ 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:

$$ 
\begin{gathered}
\min \epsilon_{max}^k + \lambda \sum_{j=1}^{n_{RM}} |U_j^k| \ , \\
\text{such that, } \forall s \in [1, n_s], x \in [1, n_x], t \in [1, n_t] \\
\epsilon_{max}^k \ge 0 \ , \\
-\epsilon_{max}^k \le A_{HF_{tx}}^{k,s} - \mathbf{X}_{tx}^s \cdot \mathbf{U}^k \le \epsilon_{max}^k
\end{gathered}
$$

Here $n_{RM} = n_c n_r n_f n_h n_q$ corresponds to the number of coefficients in $\mathbf{U}^k$; $\epsilon_{max}^k$ 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 $U_j$ in a positive and a negative part such that $U_j^k = \alpha_j^{k,+} - \alpha_j^{k,-}$. The linear program is then:

$$ 
\begin{gathered}
\min \epsilon_{max}^k + \lambda \sum_{j=1}^{n_{RM}} \alpha_j^{k,+} + \alpha_j^{k,-} \ , \\
\text{such that, } \forall s \in [1, n_s], x \in [1, n_x], t \in [1, n_t] \qquad (9) \\
\epsilon_{max}^k, \alpha_j^{k,+}, \alpha_j^{k,-} \ge 0 \ , \\
-\epsilon_{max}^k \le A_{HF_{tx}}^{k,s} - \mathbf{X}_{tx}^s \cdot \mathbf{U}^k \le \epsilon_{max}^k
\end{gathered}
$$
It is now possible to benefit from the efficiency of linear programming. Here the number of variables is equal to $1 + 2n_{RM}$ for each field $k$, and the number of constraints is $2 * n_x n_t n_s$. The number of variables is much lower than the number of constraints because the compression rate $n_{RM}$ 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 $A_{HF}^{k,s}$ or on the rows of $X^s$. This set of constraints is called *List*. It is then possible to solve the following linear problem:

$$ \begin{gathered} \min \epsilon_{max}^k + \lambda \sum_{j=1}^{n_{RM}} \alpha_j^+ + \alpha_j^- , \\ \text{such that, } \forall s \in [1, n_s], \{x, t\} \in \text{List} \tag{10} \end{gathered} $$

$$ \epsilon_{max}^k, \alpha_j^+, \alpha_j^- \ge 0 , $$

$$ -\epsilon_{max}^k \le A_{HF_{tx}}^{k,s} - X_{tx}^s \cdot U^k \le \epsilon_{max}^k $$

The number of constraints in this problem is now 2 times the length of *List*. 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 $A_{HF_{tx}}^{k,s}$ and the approximations made by the current reduced order model $X_{tx}^s$. $U^k$. 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 *List* may be considered as a validation set of the current reduced order model.

## 4.2 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 $k$,

$$ -\delta \leq (\mathbf{X}_{tx}^s - \mathbf{X}_{tx'}^s) \cdot \mathbf{U}^k \leq \delta , \\ \forall x \in C, x' \in C', t \in [1, n_t], s \in [1, n_s] $$

Where $\delta$ is a small margin in order to guarantee the solvability of the problem and $C$ and $C'$ are two set of nodes (typically the borders of each welded part). These set of constraints will ensure that the nodes in the set $C$ and $C'$ will never be further than the distance $\delta$ 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:

$$ -A_{HF_{tx}}^{k,s} - \delta \leq \mathbf{X}_{tx}^s \cdot \mathbf{U}^k \leq A_{HF_{tx}}^{k,s} + \delta , \\ \forall x \in C, t \in T, s \in [1, n_s] $$

Where $C$ and $T$ 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 $\epsilon_{max}^k$. If so, the coefficients $\mathbf{U}^k$ will ensure this error for the selected constraints in *List* 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 $\epsilon_{max}^k$ on the selected constraints in *List* 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.

<figure><img src="image_5.png" /><figcaption>Figure 4: Chart of the construction of the reduced order model. $N_i$ and $N_{max}$ stand respectively for the number of current iterations (initialized to 0) and the maximum number of iterations of the constraint generation process.</figcaption></figure>
## 4.3 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 high-fidelity 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:

<table><caption>**Table 1:** Thicknesses in mm of the 2 computations in the training set. The part number corresponds to the index of the parts on the Figure 1.</caption><thead><tr><th>Part number</th><th>P. 1</th><th>P. 2</th><th>P. 3</th><th>P. 4</th><th>P. 5</th><th>P. 6</th><th>P. 7</th><th>P. 8</th><th>P. 9</th></tr></thead><tbody><tr><th>Computation 1</th><td>0.88</td><td>0.72</td><td>0.77</td><td>1.15</td><td>1.03</td><td>1.05</td><td>0.63</td><td>0.89</td><td>0.72</td></tr><tr><th>Computation 2</th><td>0.7</td><td>0.7</td><td>0.7</td><td>1.4</td><td>1.4</td><td>1.4</td><td>0.7</td><td>0.7</td><td>0.7</td></tr></tbody></table>

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 $L_2$ 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) is 124 per field. They are the same for every cases and have been selected with the constraints generation approach described in the section 4.1 on the case with 75 column clusters. Values has been imposed to the errors $\epsilon_{max}^k$ of the constraints in the linear program. These values have been set empirically in order to minimize violations on the constraints not taken in the linear program in the case number of column clusters equal to 75. It is about 0.23 for the X normalized displacement, 0.33 for the Y normalized displacement, and 0.06 for the Z normalized displacement. Hence only the number of cluster will vary.

The figures 5 shows the evolution of the errors according to the number of columns in the time bases. The mean of the errors over all the computations seems to decrease even if there are some flat areas between 50 and 125 and between 150 and 175. 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.

<figure><img src="image_6.png" /><figcaption>Figure 5: Evolution of the RMSE according to the number of column clusters in the time base. The number of row cluster is set to 6. Each black cross stands for the error of one computation, the red crosses stand for the computations in the training set, and the circles stand for the means.</figcaption></figure>
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 $n_r > 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 $n_r > 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 $n_r = 5$ row clusters times the number of computations will be adopted.

<figure><img src="image_7.png" /><figcaption>Figure 6: Evolution of the RMSE according to the number of row clusters. The number of column is set to 20 for every cases. Each black cross stands for the error of a computation, the red crosses stand for the computations in the training set, and the circles stand for the means.</figcaption></figure>

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.
<figure><img src="image_8.png" /><figcaption>Figure 7: Evolution of the error according to the number of computation used in the reduced order model. The black crosses indicate the error on each simulation and the empty circles denote the mean of the errors. The red crosses denote the error on the computations in the training set of the reduced order model.</figcaption></figure>

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 $L_2$ error is also satisfying. The refinement technique used here give satisfactory results, but other techniques could be investigated. 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 $L_2$ error is not always the best choice to select
computations to refine the model, and perhaps choosing to add the computation
with the maximal $L_{inf}$ error would lead to a refinement in this area.

<figure><img src="image_9.png" /><figcaption>Figure 8: Final state of the crash simulation. The left figure displays the results of the high-fidelity simulation and on the right the interpolation of the same configuration using the reduced order model built with 8 computations. On the latter, the colour scale shows the difference between the high-fidelity simulation and the interpolation, red is the biggest error and blue the lowest. The displayed configuration is the worst predicted using this reduced order model.</figcaption></figure>

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.
# 5 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.

## 5.1 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.

<figure><img src="image_10.png" /><figcaption>Figure 9: Rear wheel assembly model</figcaption></figure>

Our objective was to reproduce the overall physical deformation (represented by local displacements), localisations of areas of high plastic strain as well as its maxima.
<table><caption>**Table 2:** Parameter ranges for the rear wheel assembly model</caption><thead><tr><th></th><th>Yoke</th><th>Crossbar</th><th>Arm</th></tr></thead><tbody><tr><td>Minimum thickness (mm)</td><td>4.0</td><td>2.0</td><td>2.0</td></tr><tr><td>Nominal thickness (mm)</td><td>6.0</td><td>4.5</td><td>4.0</td></tr><tr><td>Maximum thickness (mm)</td><td>8.0</td><td>6.0</td><td>6.0</td></tr></tbody></table>

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 $C_i$ and $R_i$ 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.

<table><caption>**Table 3:** Experiences properties for the rear wheel assembly case. In bold experiences used for the reduction construction.</caption><thead><tr><th>Experiences (A<sub>HF</sub>)</th><th>Yoke thickness (mm)</th><th>Crossbar thickness (mm)</th><th>Arm thickness (mm)</th></tr></thead><tbody><tr><td>1 3</td><td>4.0</td><td>2.0</td><td>6.0</td></tr><tr><td>2 7</td><td>4.0</td><td>6.0</td><td>2.0</td></tr><tr><td>3 15</td><td>6.0</td><td>4.0</td><td>6.0</td></tr><tr><td>4 20</td><td>8.0</td><td>2.0</td><td>4.0</td></tr><tr><td>5 25</td><td>8.0</td><td>6.0</td><td>2.0</td></tr><tr><td>6 28</td><td>6.0</td><td>4.5</td><td>4.0</td></tr></tbody></table>
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).

<figure><img src="image_11.png" /><figcaption>Figure 10: High-fidelity (VPS) calculation, right, and a reduced order model, left. On the reduced order model, nodal projection of the displacement error (%).</figcaption></figure>

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 high-fidelity (right) and reduced order model (left) are presented. To be usable for
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 high-fidelity 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.
<figure><img src="image_12.png" /><figcaption>Figure 11: Nodal projection of the plastic strain for a high-fidelity (VPS) calculation, right, and a reduced order model, left. Areas of high plastic deformation are indicated.</figcaption></figure>

Finding a representative metric to evaluate the quality of a model is challenging and a subject of research in itself. A try to provide such metric for this test-case is made here by calculating the average nodal error of normalized values of a field:

$$e = \frac{\sum_{i=1}^{n_x} \left| \left( \frac{A_i}{\max(|A|)} \right) - \left( \frac{A_{HF_i}}{\max(|A_{HF}|)} \right) \right|}{n_x}$$
Where $n_x$ is the total number of nodes and $A_{HF}$ 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.

<table><caption>**Table 4:** Mean normalized nodal plastic strain error at the last instant for the rear wheel assembly model. In bold, the experiments used for the construction of the reduced order model.</caption><thead><tr><td>Experiences</td><td>e<sub>plstr</sub> (%)</td></tr></thead><tbody><tr><td>1 3</td><td>0.35</td></tr><tr><td>2 7</td><td>2.58</td></tr><tr><td>3 15</td><td>0.77</td></tr><tr><td>4 20</td><td>0.76</td></tr><tr><td>5 25</td><td>3.40</td></tr><tr><td>6 28</td><td>1.08</td></tr></tbody></table>

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.

## 5.2 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, table 5).

<table><caption>**Table 5:** Experiences properties for the vehicle case. CO stands for Carry Over and ROM for Reduced Order Model.</caption><thead><tr><th rowspan="2">Experience (A<sub>HF</sub>)</th><th rowspan="2">Mass (kg)</th><th rowspan="2">CO (%)</th><th colspan="3">Intrusion (mm)</th><th colspan="3">Compression (mm)</th><th colspan="2">Error (%)</th></tr><tr><th>HF</th><th>ROM-1</th><th>ROM-2</th><th>HF</th><th>ROM-1</th><th>ROM-2</th><th>ROM-1</th><th>ROM-2</th></tr></thead><tbody><tr><td>Nominal</td><td>77.8</td><td>1</td><td>116.9</td><td>137.9</td><td>128.9</td><td>188.39</td><td>151.3</td><td>157.8</td><td>18.8</td><td>13.3</td></tr><tr><td>1</td><td>72.8</td><td>38</td><td>151.9</td><td>164.4</td><td>130.4</td><td>199.2</td><td>147.2</td><td>165.8</td><td>17.2</td><td>15.5</td></tr><tr><td>2</td><td>70.5</td><td>27</td><td>132.3</td><td>95.5</td><td>131.3</td><td>181.1</td><td>157.7</td><td>154.1</td><td>20.4</td><td>7.8</td></tr><tr><td>3</td><td>71.2</td><td>49</td><td>128.6</td><td>86.7</td><td>127</td><td>209.2</td><td>133</td><td>166.3</td><td>34.5</td><td>10.9</td></tr><tr><td>4</td><td>70.9</td><td>43</td><td>143.3</td><td>84.6</td><td>129.9</td><td>189.8</td><td>129.8</td><td>167.6</td><td>36.3</td><td>10.5</td></tr><tr><td>5</td><td>70.8</td><td>49</td><td>137.4</td><td>129.2</td><td>126.8</td><td>173.3</td><td>154.6</td><td>153.9</td><td>8.4</td><td>9.5</td></tr><tr><td>6</td><td>72.8</td><td>51</td><td>129.8</td><td>128.6</td><td>130.4</td><td>163.6</td><td>157.2</td><td>150.6</td><td>2.4</td><td>4.2</td></tr><tr><td>7</td><td>72.3</td><td>51</td><td>125.3</td><td>128.6</td><td>128.9</td><td>186.2</td><td>154.4</td><td>152.3</td><td>9.9</td><td>10.5</td></tr><tr><td colspan="5"></td><td>Mean</td><td></td><td></td><td></td><td>18.5</td><td>10.3</td></tr></tbody></table>
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 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.
<figure><img src="image_13.png" /><figcaption>Figure 12: Structural view of the model at final instant made using ROM-2 (left) and VPS (right). On the left model, nodal projection of normalized displacement error (%)</figcaption></figure>

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.
<figure><img src="image_14.png" /><figcaption>Figure 13: Comparisons of reduced order model (left column) and high-fidelity (right column) with the nodal projection of displacement error on the reduced order model. Only the parts used as parameters are shown.</figcaption></figure>

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 rendering 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.

# 6 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 refine 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 $L_1$ or $L_2$ 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 non-intrusive 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
method may also be useful to help numerical model calibration according to
experimental data. It may also be useful to exploit this reduced order model for
sensitivity analysis or control problem.

# 7 Acknowledgements

This research work has been carried out in the framework of the ROM project at IRT SystemX, Paris-Saclay, France, and therefore granted with public funds within the scope of the French Program “Investissements d’Avenir”.

# 8 References

[1] G. Berkooz, P. Holmes, J.L. Lumley, The Proper Orthogonal Decomposition in the Analysis of Turbulent Flows, Annu. Rev. Fluid Mech. 25 (1993) 539-575. doi:10.1146/annurev.fl.25.010193.002543.

[2] I.T. Jolliffe, Principal Component Analysis, Second Edition, Encycl. Stat. Behav. Sci. 30 (2002) 487. doi:10.2307/1270093.

[3] S. Chaturantabut, D.C. Sorensen, Nonlinear Model Reduction via Discrete Empirical Interpolation, SIAM J. Sci. Comput. 32 (2010) 2737-2764. doi:10.1137/090766498.

[4] M. Barrault, Y. Maday, N.C. Nguyen, A.T. Patera, An "empirical interpolation" method: application to efficient reduced-basis discretization of partial differential equations, Comptes Rendus Math. 339 (2004) 667-672. doi:Doi 10.1016/J.Crma.2004.08.006.

[5] D. Ryckelynck, F. Chinesta, E. Cueto, A. Ammar, On the a priori model reduction: Overview and recent developments, Arch. Comput. Methods Eng. 13 (2006) 91-128. doi:10.1007/BF02905932.

[6] R. Everson, L. Sirovich, Karhunen--Loève procedure for gappy data, J. Opt. Soc. Am. A. 12 (1995) 1657-1664. doi:10.1364/JOSAA.12.001657.

[7] J. Fehr, D. Grunert, Model reduction and clustering techniques for crash simulations, PAMM Proc. Appl. Math. Mech. 15 (2015) 125-126. doi:10.1002/pamm.201510053.

[8] B. Bohn, J. Garcke, R. Iza-Teran, A. Paprotny, B. Peherstorfer, U. Schepsmeier, C.A. Thole, Analysis of car crash simulation data with nonlinear machine learning methods, Procedia Comput. Sci. 18 (2013) 621-630. doi:10.1016/j.procs.2013.05.226.

[9] Y. Saad, Krylov subspace methods for solving large unsymmetric linear systems, Math. Comput. 37 (1981) 105-126.
[10] F.Z. Daim, D. Ryckelynck, A. Kamoulakos, Hyper reduction for crash simulation, in: World Congr. Comput. Mech., Seoul, 2016.

[11] K. Willcox, Unsteady flow sensing and estimation via the gappy proper orthogonal decomposition, Comput. Fluids. 35 (2006) 208–226.
doi:http://dx.doi.org/10.1016/j.compfluid.2004.11.006.

[12] D. Ryckelynck, Hyper-reduction of mechanical models involving internal variables, Int. J. Numer. Methods Eng. 77 (2009) 75–89.
doi:10.1002/nme.2406.

[13] K. Carlberg, C. Bou-Mosleh, C. Farhat, Efficient non-linear model reduction via a least-squares Petrov-Galerkin projection and compressive tensor approximations, Int. J. Numer. Methods Eng. 2010 (n.d.). doi:10.1002/nme.3050.

[14] N. Verdon, C. Allery, C. Béghein, A. Hamdouni, D. Ryckelynck, Reduced-order modelling for solving linear and non-linear equations, Int. J. Numer. Method. Biomed. Eng. 27 (2011) 43–58.
doi:10.1002/cnm.1286.

[15] F. Chinesta, A. Ammar, E. Cueto, Proper generalized decomposition of multiscale models, Int. J. Numer. Methods Eng. 89 (2010) 1114–1132.
doi:10.1002/nme.

[16] M. Vitse, D. Néron, P.A. Boucard, Virtual charts of solutions for parametrized nonlinear equations, Comput. Mech. 54 (2014) 1529–1539.
doi:10.1007/s00466-014-1073-6.

[17] L. Boucinha, A. Ammar, A. Gravouil, A. Nouy, Ideal minimal residual-based proper generalized decomposition for non-symmetric multi-field models - Application to transient elastodynamics in space-time domain, Comput. Methods Appl. Mech. Eng. 273 (2014) 56–76.
doi:10.1016/j.cma.2014.01.019.

[18] R.H. Myers, D.C. Montgomery, C.M. Anderson-Cook, Response surface methodology: process and product optimization using designed experiments, John Wiley & Sons, 2016.

[19] N. Wiener, The homogeneous chaos, Am. J. Math. 60 (1938) 897–936.

[20] N. Dyn, D. Levin, S. Rippa, Numerical Procedures for Surface Fitting of Scattered Data by Radial Functions, SIAM J. Sci. Stat. Comput. 7 (1986) 639–659.
doi:10.1137/0907043.

[21] D.E. Rumelhart, G.E. Hinton, R.J. Williams, Learning Internal Representations By Error Propagation (original), 1985.
[22] S. Poles, A. Lovison, A polynomial chaos approach to robust multiobjective optimization, in: Dagstuhl Semin. Proc., 2009.

[23] P. Feliot, Y. Le Guennec, J. Bect, E. Vazquez, Design of a commercial aircraft environment control system using Bayesian optimization techniques, in: EngOpt 2016 - 5th Int. Conf. Eng. Optim., Iguassu Falls, Brazil, 2016: pp. 19–23.

[24] M.-W. Han, P. Kopacek, Optimization using Artificial Neural Networks, IFAC Proc. Vol. 28 (1995) 357–360.
doi:http://dx.doi.org/10.1016/S1474-6670(17)46576-1.

[25] R.F. Gunst, Response Surface Methodology: Process and Product Optimization Using Designed Experiments, Technometrics. 38 (1996) 284–286.
doi:10.1080/00401706.1996.10484509.

[26] R. Benassi, J. Bect, E. Vazquez, Bayesian Optimization Using Sequential Monte Carlo, in: Y. Hamadi, M. Schoenauer (Eds.), Learn. Intell. Optim. 6th Int. Conf. LION 6, Paris, Fr. January 16-20, 2012, Revis. Sel. Pap., Springer Berlin Heidelberg, Berlin, Heidelberg, 2012: pp. 339–342.
doi:10.1007/978-3-642-34413-8_24.

[27] M. Moustapha, B. Sudret, J. Bourinet, B. Guillaume, Metamodeling for Crashworthiness Design: Comparative Study of Kriging and Support Vector Regression, in: Uncertainties 2014 - Proc. 2nd Int. Symp. Uncertain. Quantif. Stoch. Model., Rouen, France, 2014.
doi:http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.723.5147&rep=rep1&type=pdf.

[28] M. Billaud-Friess, A. Nouy, O. Zahm, A tensor approximation method based on ideal minimal residual formulations for the solution of high-dimensional problems, ESAIM Math. Model. Numer. Anal. 48 (2014) 1777–1806.
doi:10.1051/m2an/2014019.

[29] M. Haardt, F. Roemer, G. Del Galdo, Higher-order SVD-based subspace estimation to improve the parameter estimation accuracy in multidimensional harmonic retrieval problems, IEEE Trans. Signal Process. 56 (2008) 3198–3213.

[30] I. Oseledets, E. Tyrtyshnikov, TT-cross approximation for multidimensional arrays, Linear Algebra Appl. 432 (2010) 70–88.
doi:10.1016/j.laa.2009.07.024.

[31] S. Wang, Z. Zhang, Improving CUR Matrix Decomposition and the Nystrom Approximation via Adaptive Sampling, J. Mach. Learn. Res. 14 (2013) 2729–2769.
http://dl.acm.org/citation.cfm?id=2567748&CFID=190252820&CFTOK
EN=50617782.

[32] D. Amsallem, J. Cortial, C. Farhat, Towards Real-Time Computational-Fluid-Dynamics-Based Aeroelastic Computations Using a Database of Reduced-Order Information, AIAA J. 48 (2010) 2029-2037.
doi:10.2514/1.J050233.

[33] P. Benner, S. Gugercin, K. Willcox, A Survey of Projection-Based Model Reduction Methods for Parametric Dynamical Systems, SIAM Rev. 57 (2015) 483-531. doi:10.1137/130932715.

[34] M.W. Mahoney, P. Drineas, CUR matrix decompositions for improved data analysis, Proc. Natl. Acad. Sci. 106 (2009) 697-702.
doi:10.1073/pnas.0803205106.

[35] N. Mitrovic, M.T. Asif, U. Rasheed, J. Dauwels, P. Jaillet, CUR Decomposition for Compression and Compressed Sensing of Large-Scale Traffic Data, in: 16th Int. IEEE Conf. Intell. Transp. Syst., 2013: pp. 1475-1480.

[36] D.C. Sorensen, M. Embree, A DEIM Induced CUR Factorization, SIAM J. Sci. Comput. 38 (2016) A1454-A1482. doi:10.1137/140978430.

[37] J. Macqueen, Some methods for classification and analysis of multivariate observations, in: Proc. 5th Berkeley Symp. Math. Stat., 1967: pp. 281-297.

[38] H.-P. Kriegel, P. Kr{"o}ger, A. Zimek, Clustering high-dimensional data: A survey on subspace clustering, pattern-based clustering, and correlation clustering, ACM Trans. Knowl. Discov. Data. 3 (2009) 1-58.

[39] A. Nouy, Low-Rank Tensor Methods for Model Order Reduction, in: R. Ghanem, D. Higdon, H. Owhadi (Eds.), Handb. Uncertain. Quantif., Springer International Publishing, Cham, 2016: pp. 1-26.
doi:10.1007/978-3-319-11259-6_21-1.

[40] E.J. Candes, T. Tao, The Dantzig selector: Statistical estimation when p is much larger than n, Ann. Stat. 35 (2007) 2313-2351.
doi:10.1214/009053606000001523.

[41] J.S. Hesthaven, S. Ubbiali, Non-intrusive reduced order modeling of nonlinear problems using neural networks, J. Comput. Phys. (2017).

[42] D. Xiao, F. Fang, C.C. Pain, I.M. Navon, A parameterized non-intrusive reduced order model and error analysis for general time-dependent nonlinear partial differential equations and its applications, Comput.
Methods Appl. Mech. Eng. 317 (2017) 868-889.
doi:10.1016/J.CMA.2016.12.033.

[43] B. Iooss, P. Lemaître, A Review on Global Sensitivity Analysis Methods, in: Uncertain. Manag. Simulation-Optimization Complex Syst., Springer, 2015: pp. 101-122. doi:10.1007/978-1-4899-7547-8_5.

[44] M. Charrier, Y. Tourbier, L. Jézéquel, O. Dessombz, Strategic decision support through combinatorial optimization with costly evaluation function, in: Nafems World Congr. Procedings, 2017.

Acknowledgements: This research work has been carried out in the framework of IRT SystemX, Paris-Saclay, France, and therefore granted with public funds within the scope of the French Program "Investissements d'Avenir". It has been part of the ROM project.

## Vitae :

<figure><img src="image_15.png" /><figcaption>Yves Le Guennec:</figcaption></figure>

Yves Le Guennec holds a PhD in numerical and structural mechanics. He is a research engineer at the Technological Research Institute SystemX. Involved in the project ROM (Multiphysic's reduction and optimization), he focuses his research on the optimization under uncertainties and non-intrusive reduced order models.

<caption>Jean-Patrick Brunet:</caption>
<figure><img src="image_16.png" /></figure>

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:

<figure><img src="image_17.png" /></figure>

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:

<figure><img src="image_18.png" /></figure>
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.

<figure><img src="image_19.png" /><figcaption>Yves Tourbier:</figcaption></figure>

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.