Construction of a reduced-order model based on tensor decomposition and its application to airbag deployment simulations

We present a construction method for reduced-order models (ROMs) to explore alternatives to numerical simulations. The proposed method can efficiently construct ROMs for non-linear problems with contact and impact behaviors by using tensor decomposition for factorizing multidimensional data and Akima-spline interpolation without tuning any parameters. First, we construct learning tensor data of nodal displacements or accelerations using finite element analysis with some representative parameter sets. Second, the data are decomposed into a set of mode matrices and one small core tensor using Tucker decomposition. Third, Akima-spline interpolation is applied to the mode matrices to predict values within the data range. Finally, the time history responses with new parameter sets are generated by multiplying the expanded mode matrices and small core tensor. The performance of the proposed method is studied by constructing ROMs for airbag impact simulations based on limited learning data. The proposed ROMs can accurately predict airbag deployment behavior even for new parameter sets using the Akima-spline interpolation scheme. Furthermore, an extremely high data compression ratio (more than 1000) and efficient predictions of the response surfaces and Pareto frontier (2000 times faster than that of full finite element analyses using all parameter sets) can be realized.


Construction of a reduced-order model based on tensor decomposition and its application to airbag deployment simulations Takashi Sasagawa * & Masato Tanaka
We present a construction method for reduced-order models (ROMs) to explore alternatives to numerical simulations. The proposed method can efficiently construct ROMs for non-linear problems with contact and impact behaviors by using tensor decomposition for factorizing multidimensional data and Akima-spline interpolation without tuning any parameters. First, we construct learning tensor data of nodal displacements or accelerations using finite element analysis with some representative parameter sets. Second, the data are decomposed into a set of mode matrices and one small core tensor using Tucker decomposition. Third, Akima-spline interpolation is applied to the mode matrices to predict values within the data range. Finally, the time history responses with new parameter sets are generated by multiplying the expanded mode matrices and small core tensor. The performance of the proposed method is studied by constructing ROMs for airbag impact simulations based on limited learning data. The proposed ROMs can accurately predict airbag deployment behavior even for new parameter sets using the Akima-spline interpolation scheme. Furthermore, an extremely high data compression ratio (more than 1000) and efficient predictions of the response surfaces and Pareto frontier (2000 times faster than that of full finite element analyses using all parameter sets) can be realized.
Owing to the increasing computational power in recent years, product design is typically realized by taking advantage of computer-aided engineering. Moreover, multi-scale analyses 1,2 and multi-physics simulations 3 for complex nonlinear or large-scale problems may also be performed in the product design process. However, the number of such expensive simulations is limited. Therefore, optimal design or uncertainty prediction based on multi-scale or multi-physics analyses cannot be realized due to the large number of computations required. For instance, simultaneously optimizing the structure design of multiple car models requires hundreds of hours, even though the K computer is used, where 10,000 or more large-scale crash simulations are required 4 .
The optimal solution of multiple design variables can be obtained through the response surface methodology, which derives the relationship between design variables and response values sampled in advance. However, constructing a response surface directly from the results of large-scale simulations, such as the relation between the material properties and deformation behavior of all elements obtained by high-fidelity simulations, is impossible. Hence, construction methods of reduced-order models (ROMs) based on proper orthogonal decomposition (POD) have been proposed in various fields, such as crash simulations 5 , aerodynamic optimization 6 , and parameterization of material microstructures 7 . Furthermore, ROMs based on POD are applied to structurepreserving algorithms of Hamiltonian systems 8 and can provide high efficiency when combined with deep learning 9 . However, the POD is assimilated with principal component analysis or singular value decomposition, which is a low-rank factorization of two dimensional matrix data. Moreover, in most ROMs, POD projection coefficients are interpolated using the radial basis function (RBF), which requires tuning a shape parameter depending on the problem.
Tensor decomposition can factorize multidimensional data, which is a generalization of singular value decomposition. Tensor decomposition is known as the canonical polyadic (CP) model (also called PARAFAC or CandeComp) [10][11][12] , the Tucker model 13 , or the tensor-train decomposition 14 . In particular, a tensor describing multidimensional data is decomposed into a core tensor and a set of matrices in the Tucker decomposition. The core tensor is restricted to be superdiagonal in CP decomposition 15 ; that is, the CP model is a special case of www.nature.com/scientificreports/ the Tucker model. Moreover, tensor decomposition has been applied to interpolate missing data scattered over a tensor [16][17][18] . Matrix multiplication for tensor decomposition can be efficiently computed with reinforcement learning 19 . An ROM based on hierarchical tensor approximation has been proposed for Neo-Hookean problems of a simple cube and plate with a hole 20 . Furthermore, spectral tensor-train decomposition 21 , which is based on functional tensor-train decomposition and polynomial approximations, has been applied to linear elliptic partial differential equations. However, its application to strongly non-linear problems (with contact and impact behaviors) and multi-objective optimization problems has not been demonstrated yet to the best of the authors' knowledge. Furthermore, investigations on suitable interpolation methods for the ROMs of strongly non-linear problems are limited. This study aims to rapidly predict the response surfaces and Pareto frontier to reduce product development cost and time. To achieve this purpose, a construction method of the ROM for multidimensional data is proposed using tensor decomposition instead of POD. Moreover, the response surfaces and Pareto frontier in multiple design spaces are efficiently predicted without tuning parameters using piecewise linear or Akimaspline interpolation 22 instead of the RBF. This paper is organized as follows. Firstly, we present the formulation of the CP and Tucker models and then propose the construction method of the ROM using piecewise linear or Akima-spline interpolation. Secondly, we provide a numerical example of the prediction of the response surfaces and Pareto frontier in multiple design spaces with the proposed ROM. Finally, conclusions are summarized.

Construction method of reduced-order model
As previously mentioned, CP and Tucker decomposition are typical in tensor decomposition. First, the CP and Tucker models are formulated, and the construction method of the ROM is then proposed using tensor decomposition with some interpolation methods.
Tensor decomposition. The formulation of the tensor decomposition for three-dimensional data X 3D ∈ R I×J×K is shown as follows. The CP model for X 3D is defined as where X CP ∈ R I×J×K denotes the low-rank three-dimensional data obtained by the CP model, and A CP ∈ R I×R , B CP ∈ R J×R , C CP ∈ R K×R are the mode matrices corresponding to each dimension i, j, k of X 3D , respectively. In addition, E CP ∈ R I×J×K is the decomposition error caused by reducing the rank of X 3D using the CP model. In this study, X is the learning data, implying full-rank data, and X indicates low-rank data.
The Tucker model for X 3D is defined as where X Tucker ∈ R I×J×K is the low-rank three-dimensional data obtained by the Tucker model, A Tucker ∈ R I×S , B Tucker ∈ R J×U , C Tucker ∈ R K×V are the mode matrices corresponding to each dimension i, j, k of X 3D , respectively, and E CP ∈ R I×J×K is the decomposition error produced by reducing the rank of X 3D using the Tucker model. Furthermore, the core tensor G ∈ R S×U×V describes the interaction among the mode matrices A Tucker , B Tucker , C Tucker . Equation (2) is equal to Eq. (1) when S = U = V , and G is equal to the superdiagonal tensor T ∈ R S×U×V defined as Therefore, the mode matrices of the CP model are independent of each other, and the interaction among the mode matrices is considered in the Tucker model. The above formulation can be easily generalized to tensor decomposition for arbitrary n-dimensional data. In one example, the formulation of the tensor decomposition for four-dimensional data X 4D ∈ R I×J×K×L is shown.
The target of the tensor decomposition explained below returns to three-dimensional data. The mode matrices and core tensor of the CP and Tucker models are obtained by solving the following optimization problems: www.nature.com/scientificreports/ where �·� corresponds to the Frobenius norm. Although these optimization problems are solved using alternating least squares 11,13 , several local minima exist. Therefore, the initial values in the optimization are computed using higher-order singular value decomposition 23,24 .
The core consistency is defined as CP decomposition is appropriate when the core consistency is close to 100% because the core tensor is almost superdiagonal. However, the Tucker model is suitable when the core consistency is negative or close to 0% 15 .
Prediction based on tensor decomposition. Here, two methods using piecewise linear and Akimaspline interpolation are proposed to predict multidimensional data based on the mode matrices and core tensor obtained by tensor decomposition. Although the cubic-spline interpolation is normally used, the Akima-spline interpolation is applied in our method to avoid overshooting when using cubic-spline interpolation. In addition, because the CP model is a special case of the Tucker model, only prediction methods using the Tucker model are presented below. Our method can be applied to ( N + 2)-dimensional data comprising a time axis, space axis, and N parameter axes. The conceptual diagram of the proposed method when N = 1 is depicted in Fig. 1.
First, X 3D discussed in the previous section is depicted as the learning data in Fig. 1a. Second, the mode matri- Figure 1. Conceptual diagram of the proposed method to predict multidimensional data based on tensor decomposition. www.nature.com/scientificreports/ ces and core tensor are obtained using the Tucker decomposition, as shown in Fig. 1b. Third, new components of the mode matrices corresponding to the predicted values of the multidimensional data are interpolated. For instance, to predict the multidimensional data for the new parameter value p k ′ ∈ R p k < p k ′ < p k+1 , the interpolated value C predicted of the mode matrices are computed as follows: Piecewise linear interpolation: Akima-spline interpolation 22 : where k ′ is the index for the new parameter value, and p k is the parameter value of the k-th index in the parameter axis. Moreover, M kv is the slope of the segment from p k , C Tucker kv to p k+1 , C Tucker k+1,v , which is defined as To fit a cubic polynomial around the first and final parameters in the Akima-spline interpolation, the slopes preceding p 1 and beyond p K are extrapolated as follows: Furthermore, Z kv in Eq. (10) is defined as Finally, the multidimensional data X predicted for the new parameter value p k ′ are predicted as Although the prediction methods for the new parameter value are formulated based on three-dimensional data, the methods can be easily generalized to the prediction in the time or space axis and that of arbitrary n-dimensional data.

Numerical example
In this section, the proposed method is applied to airbag deployment simulations for validation purposes. First, the analysis conditions of the airbag deployment simulations are summarized, and the construction procedure for learning the data from the results of the airbag deployment simulations is then presented. Second, ROMs are constructed based on the learning data to predict the airbag deployment behavior for new parameter sets. Moreover, the accuracy of the proposed method is investigated by comparing the prediction results with highfidelity computational results. Finally, the response surfaces and Pareto frontier in multiple design spaces are efficiently predicted with the ROMs.
Model description for airbag deployment simulation. To construct learning data for ROMs, airbag deployment simulations were performed using the sample model 25 produced by JSOL Corporation. The deployment and impact behavior of the airbag and impactor were computed using the finite element method, and the motion of the inflation gas was modeled by the corpuscular particle method 26 . The finite element model is displayed in Fig. 2. The airbag fabric, tether, and vent holes were discretized by membrane elements, and the inflator case and impactor were discretized using shell elements. The model consists of 330 nodes and 3352 elements. The material parameters of each component are listed in Table 1. The fabric with the liner was set orthotropic elastic, and the inflator and impactor were considered rigid. The mass density of the inflator case was 8.0 × 10 −9 ton/mm 3 , and the weight m of the impactor was 5.0 × 10 −3 ton . The elastic moduli and Poisson ratios of the inflator case and impactor were 2.0 × 10 5 MPa and 0.3, respectively. The contact www.nature.com/scientificreports/ of the fabric with itself and the impactor was considered. Furthermore, the inflator case was fixed, and the initial velocity in the x-direction of the impactor was 2.0 × 10 4 mm/s. Subsequently, analysis conditions for the corpuscular particle method are indicated. Air of 293 K and 1.0 atm was present inside and outside the airbag as the initial gas. The gates of the inflation gas were located concentrically at every 45 • on the y-z plane, as shown in Fig. 2, and the distance between the center of the inflator case and each gate was 30 mm. Nitrogen gas of 800 K flowed into the airbag from the gates for T inflation s at a mass flow rate of 2.0 × 10 −3 ton/s . The properties of the air and nitrogen gas are displayed in Table 2. Their behavior was simulated using 20,000 particles (7000 for the initial gas inside the airbag and 13,000 for the inflation gas). This number of particles is high enough to describe a smooth pressure distribution on the fabric elements 27 . Moreover, the vent hole coefficient of Wang-Nefske leakage 28,29 was 1.0.
The length L vent of the vent holes and inflation time T inflation were considered as the design variables. The airbag deployment simulations with the parameters shown in Fig. 3 were performed to obtain the learning and validation data for prediction using our method. The termination time of the simulations was 0.1 s, and the time step size was 9.0 × 10 −7 s . The simulations were performed with four cores of the Intel Xeon E5-2667v4 processors using massively parallel processing in the general-purpose finite element software LS-DYNA (R11.1.0).
Construction procedure for multidimensional learning data. We performed a numerical example to investigate the prediction accuracy of the deformation and energy absorption of the airbag, as well as the motion and maximum acceleration of the impactor. The time history response of each nodal displacement described the deformation behavior of the airbag motion of the impactor. Energy absorption was computed based on the initial and final velocities of the impactor. The maximum acceleration of the impactor was indicated from the time history response of its acceleration. Furthermore, the displacement and velocity of the impactor were computed by integrating its acceleration. Therefore, the multidimensional learning data X 4D were constructed relying on the time history response of each nodal acceleration of the impactor and that of each nodal displacement of the   www.nature.com/scientificreports/ airbag. The overline indicates that the variables are computed in this study using airbag deployment simulations. Note that the nodal acceleration of the impactor was employed because the velocity and acceleration computed by differentiating the displacement included noise and vibration. The learning data are four-dimensional data consisting of a time axis, space axis, first parameter axis, and second parameter axis representing the time steps, node numbers, vent hole size L vent , and inflation time T inflation , respectively. The matrix data X kl 4D on the time and space axes included in X 4D are displayed below: where I is the total number of time steps, N is the total number of nodes, and d n i is the displacement vector of the n-th node at the i-th time step, which is defined as where d x , d y , and d z are the nodal displacements in the x, y, and z directions, respectively. Therefore, the total number J of the components of X 4D on the space axis is three times larger than the total number of nodes, that is, J = 3N . Subsequently, the total numbers K, L of the components of X 4D on the first and the second parameter axes are four, which is the number of levels of L vent and T inflation for the learning data, respectively. The matrix data X predicted k ′ l ′ on the time and space axes included in X predicted are described in the same way as Eq. (18). The total numbers K ′ , L ′ of the components of X predicted k ′ l ′ on the first and second parameter axes are three when predicting for the white circles indicated in Fig. 3.
Validation of prediction accuracy. By following the procedure explained in the previous section, the multidimensional learning data X 4D were constructed from the results obtained by the airbag deployment simulation, as indicated in Fig. 3. Subsequently, the core consistency was computed with X 4D using the N-way toolbox 30 in the programming and numeric computing platform MATLAB (R2016b), as listed in Table 3. The results demonstrated that the Tucker model is suitable for this learning data because the core consistency becomes negative or close to zero when the size of the core tensor is larger than (S, U, V , W) = (2, 2, 2, 2) . In particular, the validation results with the Tucker model are shown below.
The relationship between data compression ratios and errors using Tucker decomposition with different sizes of the core tensor is summarized in Fig. 4. The sizes of the core tensor were determined using the function "hosvd" in the Tensor Toolbox 31 with the following different tolerances: 0.01, 0.02, 0.03, 0.05, 0.1, and 0.2. The data compression ratio was computed by dividing the total number of components of the mode matrices and the core tensor by the total number of components of X 4D , which is defined as www.nature.com/scientificreports/ Furthermore, the decomposition and prediction errors were obtained using the following equations: where X predicted is the airbag deployment simulation results using the parameter sets regarding the white circles indicated in Fig. 3. The data compression ratio was 14,545, which was the maximum when (S, U, V , W) = (1, 1, 1, 1) . As the size of the core tensor increased, the data compression ratio and errors decreased. However, the prediction errors reached their limits before reaching zero, whereas the decomposition error asymptotically approached zero. The optimal size of the core tensor was identical to (17,8,4,4) because the prediction errors did not change substantially when (S, U, V , W) ≤ (17, 8, 4, 4) . Additionally, the data compression ratio was 1614 when (S, U, V , W) = (17, 8, 4, 4). Tables 4 and 5 list the decomposition and prediction errors, respectively, regarding each parameter set for (S, U, V , W) = (17, 8, 4, 4) . Moreover, the decomposition and prediction errors were obtained by fixing the parameter indices (k, l) in Eq. (21) and k ′ , l ′ in Eq. (22), respectively. The airbag deployment simulation results were predicted within 17% and 9% of the error using piecewise linear and Akima-spline interpolation, respectively. Furthermore, the errors tend to be large when the inflation time is small.
The deployment and impact behavior of the airbag and impactor when the prediction error using the Akimaspline approximation is at maximum; that is, at L vent = 82.5 mm and T inflation = 0.01 s, are summarized in Fig. 5. The errors in the contour maps show the difference in the nodal displacements between the finite element analysis and the proposed ROMs. The simulation results with the proposed ROMs agreed with those obtained by finite element analysis until t = 0.025 s. The errors of the impactor displacements at t = 0.07 s were large because of the irregular behavior of the impactor; that is, when the inflation time was short, the impactor hit the inflator. Figure 4. Relationship between the data compression ratio and the Tucker decomposition error. www.nature.com/scientificreports/  Prediction of response surfaces and Pareto frontier. The proposed ROMs can be employed in numerous parametric studies and response surfaces and Pareto frontiers in multiple design spaces. To predict the response surfaces and Pareto frontier, the deployment and impact behavior of the airbag and impactor were predicted using the Tucker models and Akima-spline approximation with the parameter sets indicated in Fig. 6.   www.nature.com/scientificreports/ The maximum acceleration of the impactor and the energy absorption E absorption computed using the following equation for the parameter sets are depicted in Fig. 7. The values are calculated using the following equations: where m, v initial , and v final are the mass, initial velocity, and final velocity of the impactor, respectively. v final was computed by integrating the acceleration of the impactor included in X predicted . The maximum acceleration of the impactor rapidly increased when the inflation time was shorter than approximately 0.01 s. The rapid increase in the maximum acceleration was caused by the collision of the impactor with the inflator when the inflation time was short. Figure 7a also shows that the maximum acceleration does not strongly depend on the vent hole size. Moreover, the energy absorption increased as the vent hole size and inflation time increased, as shown in Fig. 7b. The energy absorption was almost zero when L vent = 0 mm, that is, without a vent hole. In addition, the energy absorption became negative when the inflation time was long without the vent hole because the gas inflation energy was transferred to the impactor. The sensitivity of the vent hole size and the inflation time to the energy absorption was investigated by defining the following index g. Moreover, the relationship between the energy absorption and g is displayed in Fig. 8: where the maximum vent hole size L max vent is 99 mm, and the longest inflation time T max inflation is 0.02 s. g is the index that denotes the summation of the normalized L vent and T inflation . The Pareto frontier regarding the energy  www.nature.com/scientificreports/ absorption and product cost is depicted in Fig. 8 when the product cost of the airbag is proportional to the vent hole size and inflation time. Note that while the Pareto frontier with the learning data is unclear, the frontier with the predicted data using the ROMs is clear. Finally, the computational costs are summarized. The average CPU time for the finite element analysis with a parameter set was approximately 114 min. Therefore, the total CPU time for the 130 parameter sets shown in Fig. 6 can be estimated to be 247 h. In contrast, the computational time with the ROMs was only approximately 7 min; that is, the response surfaces and Pareto frontier could be predicted 2000 times faster than the finite element analyses.

Conclusion
In this study, a method for constructing ROMs for strongly non-linear problems involving contact and impact behaviors was proposed using tensor decomposition without tuning any parameters. The proposed method was verified and validated through airbag deployment simulations. Moreover, application to multi-objective optimization and suitable interpolation methods was investigated. First, finite element analyses of airbag deployment with 16 parameter sets were performed. Subsequently, four-dimensional learning data for the ROMs were constructed based on the time history response of each nodal acceleration of the impactor model and each nodal displacement of the airbag model. Second, the ROMs for the airbag deployment simulations were constructed by tensor decomposition of the learning data. Airbag deployment simulation results with nine other parameter sets were predicted using the ROMs with piecewise linear or Akima-spline interpolation. The error with the Akima-spline interpolation was approximately half of that with piecewise linear interpolation because the new parameter sets using the Akima-spline interpolation scheme could avoid unrealistic overshoots that occurred when using the standard cubic spline interpolation procedure. Finally, the response surfaces and Pareto frontier were efficiently predicted by computing the airbag deployment simulation results with 130 other parameter sets with the ROMs, which is 2000 times faster than full simulations required for finite element analyses. In conclusion, the proposed method can reduce product development cost and time by rapidly predicting the response surfaces and Pareto frontier. Our method is currently limited to inefficient grid sampling for gathering learning data; therefore, an extension of our method to random sampling will be considered in future work.

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request. www.nature.com/scientificreports/