Mathematical Modeling and Parameters Estimation of Car Crash Using Eigensystem Realization Algorithm and Curve-Fitting Approaches

An eigensystem realization algorithm (ERA) approach for estimating the structural system matrices is proposed in this paper using the measurements of acceleration data available from the real crash test. A mathematical model that represents the real vehicle frontal crash scenario is presented. The model’s structure is a double-spring-mass-damper system, whereby the front mass represents the vehicle-chassis and the rear mass represents the passenger compartment. The physical parameters of the model are estimated using curve-fitting approach, and the estimated state system matrices are estimated by using the ERA approach. The model is validated by comparing the results from the model with those from the real crash test.


Introduction
Car crash test is usually performed in order to ensure safe design standards in crashworthiness (the ability of a vehicle to be plastically deformed and yet maintain a sufficient survival space for its occupants during crash scenario).Nowadays, due to advanced research in computer simulation software, simulated crash tests can be performed beforehand the fullscale crash test.Therefore, cost associated with real crash test can be reduced.Vehicle crashworthiness can be evaluated in four distinct modes: frontal, side, rear, and rollover crashes.Several researches have been carried out in this field, which resulted in several novel computational models of vehicle collisions in literature.In [1], a mathematical model is proposed to estimate the maximum occupant deceleration, which is one of the main tasks in the area of crashworthiness study by a Kelvin model which contains a mass together with spring and damper connected in parallel.An application of physical models composed of springs, dampers, and masses joined together in various arrangements for simulating a real car collision with a rigid pole was presented in [2].
In [3], the authors presented an overview of the kinematic and dynamic relationships of a vehicle in a collision, whereby the work was to identify the parameters of the vehicle crash model using experimental dataset.In [4,5], a lumped parameter modeling in frontal crash was investigated and analyzed in five degrees of freedom and has been used to analyze the response of occupant during the impact.In [6,7], an optimization procedure to assist multibody vehicle model development and validation was proposed.The authors first devised the topological structure of the multibody system representing the structural vehicle components and described the most relevant mechanisms of deformation.In the work of [8], the authors proposed an approach to control the seat-belt restraint system force during a frontal crash to reduce thoracic injury.
The main challenge in accident reconstruction is the system identification described as the process of constructing mathematical models of dynamical systems using measured input-output data.In case of vehicle crash, system identification algorithm consists in retrieving the unknown parameters such as the spring stiffness and damping coefficient.A possible approach is to identify these parameters directly from experimental dynamical data.
From literature, system identification algorithms (SIA) have been developed for different applications.Among

Vehicle Crash Experimental Test
The real vehicle crash experiment was conducted on a typical mid-speed vehicle for pole collision.Its elaboration was the initiative of Robbersmyr [19].A test vehicle was subjected to impact with a vertical, rigid cylinder.The acceleration field was 100 meters long and had two anchored parallel pipelines.The vehicle was steered using those pipelines that were bolted to the concrete runaway.Setup scheme is shown in Figure 1.
During the test, the acceleration was measured in three directions (-longitudinal, -lateral, and -vertical) together with the yaw rate from the center of gravity of the car.Using normal-speed and high-speed video cameras, the behavior of the safety barrier and the test vehicle during the collision was recorded.The initial velocity of the car was 35 km/h, and the mass of the vehicle (together with the measuring equipment and dummy) was 873 kg.The obstruction was constructed with two steel componentsa pipe filled with concrete and a baseplate mounted with bolts on a foundation.The car undergoing the deformation  is shown in Figure 2. The accelerometer is located at the mass center of gravity of the vehicle in the passenger compartment.Since we are interested in the frontal crash, only the measured acceleration in the longitudinal direction is considered in this study.The acceleration data is imported and processed in MATLAB for analysis.The deformation of the vehicle is obtained by integrating twice the acceleration signal.

Mathematical Modeling Theoretical Background
Mathematical models describe the dynamic behavior of a system as a function of time.During frontal crash, the vehicle is subjected to an impulsive force caused by the obstacle.The model for vehicle crash simulates a rigid barrier impact of a vehicle, where  1 and  2 represent the frame rail (chassis) and passenger compartment masses, respectively.Parameters to be estimated are springs  1 and  2 and dampers  1 and  2 , as shown in Figure 3.When the vehicle impacts on a rigid barrier, the two masses will experience an impulsive force during collision.The method for solving the impact responses of the two masses is adapted from the method used in the free vibration analysis of a two-degreesof-freedom damped system [20].
The dynamic equations of the double-mass-springdamper model are shown as follows: or The solution for   ,  = 1, 2, can be represented as with  = 1, . . ., 4.
Where   and   may be complex numbers.Substituting (3) into (1), we get where After substituting  1 ,  2 ,  1 , and  2 into (4), we get the follwing: Now, for a nontrivial response, that is, for nonzero values of  1 and  2 , the determinant of their coefficient matrix must vanish.That is, Expansion of (10) leads to a characteristic equation of the system obtained as follows: where Equation ( 11) is a fourth-order polynomial in s and is to be solved to get four roots.All of the coefficients of this polynomial are physical parameters of the system shown in Figure 3 and are all positive.For that reason, such a polynomial cannot have positive roots.Three allowable configurations of roots are given as follows [20]: (1) two pairs of complex conjugates, (2) one pair of complex conjugates and two real and negative roots, (3) four real and negative roots.
Case 1 (two pairs of complex conjugates).The system in this case has moderate damping.The rate of decay is defined by  1 , the real part of the root, and the frequency of vibration is specified by  1 , the imaginary part.The two pairs of complex conjugates are the following: (1) where  1 ,

Mathematical Problems in Engineering
Case 2 (one pair of complex conjugates and two real and negative roots).The general displacement solutions are shown as follows: where  = 1, 2.
Case 3 (four real and negative roots).The system has a large damping.When it is disturbed, the system will settle to its equilibrium configuration without oscillation.The solutions of the 4th-order polynomial yield four real and negative roots.
In [2], the authors focused just on the first case.In this paper, we will also focus on the third case.The third case is for the system which has a large damping.When it is disturbed, the system will settle to its equilibrium configuration without oscillation.The displacement signal of the real crash is similar to that of a case of an overdamped vibrating system.Hence, the third case would represent the vehicle frontal crash reconstruction.The solution for Case 3 is the following: with  = 1,2, where   are the roots of the characteristic equation (11).

Estimation of Model
Parameters by Curve Fitting.When ( 16) is curve fitted into displacement experimental data, the constants   and   ( = 1, 2, and  = 1, . . ., 4) can be easily found resulting in the following system of equations that can be solved for  2 and  2 : For  = 1, . . ., 4, (17) can be written in a matrix form as Let Then, ( 19) can be represented as and, using the pseudoinverse, we can obtain the following: The spring stiffness  1 and damping coefficient  1 are calculated from (25) as follows: Remark 1. From ( 6), it was shown that Then, we have that Therefore, the estimated parameters are obtained from ( 22) and (25).

Eigensystem Realization Algorithm.
In this section, a state-space representation (SSR) of the model is derived from the dynamic equation of a vehicle subjected to a frontal crash.The same representation is also retrieved from the system Markov parameters.

Formulation of the SSR from the Model Dynamic
Equations.Considering an input force u 1 acting on the front mass  1 , (1) can be rewritten as The dynamic equation ( 26) can be rewritten in a matrix compact form as with In general, the equation of motion for  degrees of freedom is expressed in a matrix form as or, equivalently, where M ∈ R × , L ∈ R × , and K ∈ R × are the mass, damping, and stiffness matrices, respectively, while B  ∈ R × is the input matrix.ẍ , ẋ , and x are vectors of generalized acceleration, velocity, and displacements, respectively, and the vector u of dimension  × 1 is the input force containing  external excitations acting on the systems.
Let us define the new state variables x 1 =  1 , x 2 = ẋ 1 , x 3 =  2 , and x 4 = ẋ 2 .Substituting these variables into (26) and combining in state equations, we get the following: Using the original state variables and interchanging rows 2 and 3, columns 2 and 3 of (31) and we get the following: The output equation or the measurement vector y(), which may contain any combination of modal displacements, velocities, and/or accelerations, is given by where C  , C V , and C  are the output influence matrices for position, velocity, and acceleration, respectively.In our experiment, the acceleration is measured.Therefore, the output equation is the acceleration measurement From (32), Therefore, the continuous-time state-space model of the dynamic system is written as with where A  is the state matrix, B  is the input matrix or the state influence matrix, C  is the output matrix or the measurement influence matrix, and D  is the feedforward matrix or the direct transmission matrix.
Once the A  , B  , C  , and D  matrices are known, it is easy to find the transfer function (TF) and impulse response function (IRF) of the system.By using these matrices, the system's response to any input can be found in time domain or frequency domain.The state-space representation is useful for constructing the mathematical model in MATLAB environment.
The discrete-time state-space representation of a MIMO system is given by the following: where  is the integer discrete-time index at time instant  = Δ, x() is the state vector at the discrete-time , u() is the force vector, y() is the output vector, Â is the discrete state system matrix, and B is the discrete input influence matrix for the state vector x().The output matrix Ĉ = C  and the direct transition matrix D = D  during the zero-order-hold operations.Because experimental data are discrete in nature, equations (38) form the basis for the system identification of linear time-invariant, dynamical systems.The state matrix Â and the influence matrix B of the discrete-time model are related to the matrices A  , B  of the continuous-time model by the following expression [14]: The continuous-time model is calculated from the discretetime model by Δ , where Δ is a constant interval.
The dimensions of the discrete-time system are equal to those of the continuous system.

ERA from the System Markov Parameters.
ERA is a minimum-order realization technique that uses singular value decomposition technique.A flowchart for the ERA is shown in Figure 4.
If the excitations of the dynamic system are measured by the  input quantities in the vector u, the equations of motions and the set of output equations can both be, respectively, rewritten in terms of the state vector.
ERA begins with the definition of the Markov parameter of a state-space model.The method for deriving the expression for the system matrices is adapted from [14].Consider a discrete-time state-space model (38).
The state-space model (38) has an impulse response The discrete-time Markov parameters can be defined in the same way as (41).The term of ĈÂ −1 B is called the Markov parameter of the system.By using these parameters, one can define the impulse response of the system as follows: Consequently, the identification problem is the following: given the values of []  s, construct the constant matrices to identify the system.
The algorithm begins by constructing an  ×  generalized Hankel matrix.
Given a number of input and output measurements   and   generated by a system of unknown parameters, it is requested to identify the order of the system as well as the discrete state matrices ( Â, B, Ĉ, D) of the system.Then, the identified continuous state matrices (A  , B  , C  , D  ) that have the same size as in the physical model can be estimated from (40), and, finally, extract the respective parameters.
All minimum realizations have the same set of eigenvalues and eigenvectors, which are the modal parameters of the system itself.Assume that the state matrix Â of order  has a complete set of linearly independent eigenvectors {Ψ 1 , Ψ 2 , . . ., Ψ  } with corresponding eigenvalues { 1 ,  2 , . . .,   }; then where Λ is a diagonal matrix of the eigenvalues and Ψ is the matrix of the eigenvectors.The realization { Â, B, Ĉ} can be transformed into the realization {Λ, Ψ −1 B, ĈΨ} by using the eigenvalues and eigenvectors matrices.The diagonal matrix Λ contains the information of modal damping rates and damped natural frequencies.The matrix Ψ −1 B defines the initial modal amplitudes, and the matrix Ψ denotes the mode shapes at the sensor points.All of the modal parameters of a dynamic system can, thus, be identified by the triplet {Λ, Ψ −1 B, ĈΨ}.
The real part of Λ, into the continuous-time model, gives the modal damping rates, while the imaginary part gives the damped natural frequencies.
After identifying the combined system and observer gain Markov parameters, the following step consists in forming the generalized Hankel matrix  (−1) : (45) When  = 1, we get the following: In order to compute a minimum-order realization of the system { Â, B, Ĉ}, it is necessary to construct a shifted Hankel matrix  (1) as follows: Substituting the Markov parameters from (43) into (45) and decomposing ( − 1) into three matrices yield the following: Hankel matrix H (1) Figure 4: Flowchart for the ERA by Juang [14].
where O  and C  are given as The block matrix O  is the observability matrix, whereas the block matrix C  is the controllability matrix.
Denote column submatrices of B by B and row submatrices of Ĉ by Ĉ .The ERA data block matrix can be expressed by Assume that there exists a matrix  † satisfying the relation where   is an identity matrix of order .The matrix  † which is the pseudoinverse of (0) plays a major role in deriving the ERA.It is observed that The ERA process starts with the factorization of the block data matrix (46) using the singular value decomposition where the columns of matrices  and  are orthonormal and ∑ is a rectangular matrix given as with ∑  = diag{ 1 ,  2 , . . .,   ,  +1 , . . .  }.
Let   and ∑    be matrices formed by the first  columns of R and S, respectively.Hence, the matrices (0) and  † become where      =   =      .Examining the singular value ∑  of the Hankel matrix  (0) , it is possible to determine the order of the system.
Comparing ( 56) and ( 51) with  = 0, we get that O  =   ∑ 1/2  and C  = ∑ 1/2     .From (51), the first  columns of the observability matrix O  form the input matrix B, whereas the first  rows of the controllability matrix C  form the output matrix Ĉ.

Mathematical Problems in Engineering
With  = 1 in (51), we get that One obvious solution for the state matrix Â becomes Let   be a null matrix of order , let   be an identity matrix of order , and let the matrices E   and E   be defined as where  is the number of outputs and  is the number of inputs.
Finally, using (50), ( 51), ( 52), (56), and (57), the basic formulation of the minimum-order realization for the ERA/OKID (OKID means Observer/Kalman filter identificationis) is given as follows: Hence, a minimal realization is given as follows: The continuous state matrix A  and the input influence matrix B  are obtained from (40), and the physical parameters-matrices M, L, and K are embedded in the state matrix A  .The result from the curve fitting is shown in Figure 5.

Vehicle Crash Experimental Data Analysis.
It is observed from Figure 6 that the dynamic crash from the real vehicle crash test is 53.17 cm and occurs at time   = 0.078 ms, when the unfiltered data are used in the analysis.The filtered data result in a dynamic crash of 51.11 cm at time   = 74.5 ms as shown in Figure 7.
The initial velocity for both filtered and unfiltered data is closer to 35 km/h (i.e., 34.99 km/h for the unfiltered data and 35.28 km/h for the filtered data).

4.3.
Results from the Model.Four different cases are considered in this section as a sample of results.Let  1 be the mass of the chassis,  2 the mass of passenger compartment, and   = 873 kg the total mass of the vehicle.
Case 1 ( 1 <  2 ).One has that  1 = (1/3)  , and  2 = (2/3)  .From Figure 8, the dynamic crash of  2 is 80 cm which is the displacement of the passenger compartment.Therefore, this model cannot represent the vehicle crash scenario.It is observed that the time for dynamic crash is longer than that for the real crash (i.e., 0.17 s instead of 0.078 s).The dynamic crash of  1 is 42.5 cm and occurs after 0.17 s.
Case 2 ( 1 >  2 ).One has that  1 = (2/3)  , and  2 = (1/3)  .From Figure 9, the dynamic crash of the passenger compartment  2 is 66.2 cm.The time for dynamic crash increases further up to 0.15 s.The dynamic crash of  1 is  45.5 cm and occurs after 0.13 s.Therefore, for this case, the model cannot represent the vehicle crash scenario.
Case 3 ( 1 <  2 ).One has that  1 = (1/4)  , and  2 = (3/4)  .From Figure 10, the dynamic crash of the passenger compartment  2 is 69 cm.The time for dynamic crash is 0.15 s.The dynamic crash of  1 is 30.9 cm and occurs after 0.14 s.This also cannot represent the real vehicle crash test.
Case 4 ( 1 >  2 ).One has that  1 = (3/4)  , and  2 = (1/4)  .From Figure 11, the dynamic crash of the passenger compartment  2 is 49.8 cm, and the time for dynamic crash is 0.11 s.The dynamic crash of  1 is 35.5 cm and occurs after 0.1 s.Therefore, this case can represent the vehicle crash senario because the dynamic crash is much closer to that from the real vehicle crash and the time is relatively small as compared to other cases.A summary of main results is shown in Table 1.The values for  2 and  2 are the first and second entries of vector V 2 in (22).The values for  2 and  2 depend on the value of mass  2 taken into consideration.Values for  1 and  1 are obtained from vector V 1 in (25).The stiffness coefficients which result in a closer vehicle crash reconstruction are found to be  1 = 74681 N/m, and  2 = 45821 N/m, and the damping coefficients are:  1 = 18176 Ns/m and  2 = 11196 Ns/m when the mass of the chassis is 3/4 the total mass of the vehicle (  ), where the dynamic crash of the passenger compartment is equal to 49.8 cm and occurs after 0.11 s (see Subsection 4.3 , Case 4, Figure 11).
It is observed that the passenger compartment  2 is the one that reconstructs the vehicle crash.When the mass of the From the state matrix A  , the eigenvectors Ψ and eigenvalues Λ were found to be and Λ = diag{−18.04± 533.15}.
The natural frequency and the damping ratio of a single degree of freedom are 533 rad/s and 0.0338, respectively.The Bode diagram for the system is shown in Figure 12.
Considering a 4th-order system (i.e, = 4) for a twodegree-of-freedom system and the number of samples to assemble the Hankel matrix equal to  = 60, it is observed Frequency (rad/s)  that the continuous state-space realizations from the ERA are found to be the following:   One can see the stability of the model from pole-zero map as shown in Figure 13 and Figure 14.It is worth noted that for small sample size the system becomes unstable as evidenced by Figure 15.
The natural frequencies and damping ratios are given as   1 = 617 rad/s,   2 = 575 rad/s and  1 = 0.00712,  2 = 0.006680.Remark 3. Based on the state-space realization obtained by ERA, we should be able to transfer the state matrix A  to the original obtained from (36).Hence, extract the mass, stiffness, and damping matrices.

Conclusion and Future Work
In this paper, a method to estimate the parameters of a double-spring-mass-damper model of a vehicle frontal crash is presented.It was observed that the model results were much closer to the experimental test results.The overall behavior of the model matches the real vehicle's crash.Two of the main parameters characterizing the collision are the maximum dynamic crash-which describes the highest car's deformation-and the time at which it occurs-  .They are pertinent to the occupant crashworthiness since they help to assess the maximum intrusion to the passenger's compartment.
It can be concluded from this study that a doublespring-mass-damper model can represent the real vehicle crash scenario when the mass  1 representing the chassis of the vehicle is less than  2 , the mass representing the passenger compartment.In all cases, the front part of the vehicle undergoes smaller deformation than the passenger compartment.The time at which the maximum chassis displacement occurs is slightly shorter than the time for the passenger compartment because of its additional compression by the rest of the car.In our future work, (1) more investigation on data analysis would be interesting to consider effect of sensor delays in measurement [21][22][23]; (2) extraction of mass, stiffness, and damping matrices from the identified state-space representation will be investigated; (3) we would like to extend our study to a three-mass-springdamper model taking into consideration the nonlinearity of the spring and damper: the three masses will be representing the engine, the suspension, and the passenger compartment mass, respectively, interconnected by springs and dampers; (4) injury mechanisms of the occupant such as Head Injury Criterion (HIC) would be considered in the future work.

10 Figure 7 :
Figure 7: Filtered data plot from the real vehicle crash test.
Figure 11: Comparative analysis between vehicle crash test and model results for  2 = (1/4)  .chassis is greater than that of the passenger compartment, the results from the model are closer to the expected values.For example, when  1 = 582 kg (i.e., 3/4 of the total mass of the vehicle) and  2 = 291 kg (i.e., 1/4 of the mass of the vehicle), the dynamic crash of the passenger compartment is 49.8 cm which is closer to 51.11 cm (the dynamic crash from the real vehicle crash).