Modeling and parameter identification of seated human body with the reference vector guided evolutionary algorithm

The low frequency vibration of the vehicle in motion has a great influence on the ride comfort of occupants. The research on the vibration response characteristics of human body plays a great role in analyzing and improving ride comfort. The purpose of this study was to investigate the parameter identification of seated human body dynamic model. A seven-degree-of-freedom (DOF) lumped parameter model was established to describe the vibration response characteristics of human body. The derivation processes of apparent mass (AM) and seat to head transmissibility (STHT) were performed. After the theoretical calculation of the human body vibration characteristics, we used several different evolutionary algorithms to identify the 23 parameters of the model, including the mass, stiffness and damping parameters. By comparing the results of the five optimization algorithms and comprehensively analyzing the convergence and distribution of the non-dominated solution set, we found that the reference vector guided evolutionary algorithm (RVEA) shows good competitiveness in solving many-objective optimization problem (MaOP), that is, parameter identification of seated human body model in this paper. The AM and STHT calculated by model identification were compared with their measured by experiment. The result shows that the selected seven-DOF model can well describe the vertical vibration characteristics of seated human body and the identification method used in this paper can accurately identify the parameters of lumped parameter model, which provides convenience for the establishment of a complete “road-vehicle-seat-human body” system dynamic model.


Introduction
Vibrations of vehicle cab caused by road excitation and powertrain operation have bad effects on ride comfort of occupants. In order to improve vehicle ride comfort, damping parts such as suspension, tire and seat, were constantly being studied and optimized. Meanwhile, human body was usually just equivalent to a rigid body with a fixed mass (sandbag, water tank). However, there are great differences in dynamic vibration response characteristics between human body and rigid body due to the complex nonlinear mechanical properties of human body, especially when the excitation frequency is greater than 2 Hz. [1][2][3] Therefore, the research on vertical vibration characteristics of seated human body is an important part of the overall dynamic model of road-vehicle-seat-human body system, and has a positive effect on the ride comfort improvement. [4][5][6] Vehicle vibrations can lead to discomfort and fatigue of occupants, and even cause damage to human body. Studies have shown that seated occupants are particularly sensitive to whole-body vibration caused by low frequency excitation. 6 Therefore, the dynamic response characteristics of human body under low frequency vertical vibration excitation have been paid special attention by scholars. Three biodynamic response functions were usually used to analyze the vibration characteristics of human body -seat to head transmissibility (STHT), apparent mass (AM), and mechanical impedance (MI) respectively. [7][8][9][10][11][12] In order to analyze the vibration of human body and improve ride comfort of vehicle, various models have been built to describe the vibration response characteristics of human body. [13][14][15][16] Lumped parameter model is a typical, simple method to analyses vibration responses of seated human body. The number of degrees of freedom (DOF) of the model gradually increases from 1 to 15 with the development and refinement of research. [17][18][19][20][21][22] Based on extensive experimental studies, ISO 5982-2019 describes the ideal range of human biodynamic responses under whole-body vibration, and recommends a three-DOF model of seated human body. 23 Apparent mass is more convenient to measure accurately and widely used in analysis of whole-body vibration. 24,25 No matter AM, STHT, or MI, they're all frequency response functions used to describe biodynamic characteristics. Therefore, the amplitudefrequency and phase-frequency characteristics of the model should be identified simultaneously.
Many previous research methods took the sum of the errors of different objective functions and test data as the overall goal, and then applied global optimization algorithm to identify the model parameters. 26,27 However, the improvement of these design objectives may conflict with each other, resulting in the tread-off relationship between different objective functions. Obviously, the single objective optimization algorithm may run into locally optimal solution. Therefore, the multi-objective optimization algorithm should be used for model identification of seated human body.
Multi-objective evolutionary algorithms (MOEAs) have many subcategories. 28,29 When the number of objectives exceeds three, the problems are called manyobjective optimization problems (MaOPs). The typical MOEAs, such as the decomposition-based evolutionary algorithm (MOEA/D) and the dominant relation-based evolutionary algorithm (NSGA II), are considered difficult to get a satisfactory result for solving manyobjective optimization problem. Therefore, most current algorithms are based on a combination of the dominant relation and decomposition method (NSGA-III, u-DEA, MOEA/DD, RVEA, et al.). [30][31][32][33][34] Obviously, the parameter identification of human dynamics model is a many-objective optimization problem.
In the following section, apparent mass and seat to head transmissibility calculation method of seven DOF lumped mass parameter model was calculated. Thereafter, in section 3, the algorithm flows of RVEA were introduced. In section 4, different optimization algorithms were used to identify the parameters of human dynamics model, and the calculation results of each algorithm were compared. The final model identification results are given in section 5. The results verify the rationality and accuracy of the analysis and calculation method in this paper. The conclusions were presented in section 6.

Biodynamic model for seated human body vibration
In literature, many dynamic models have been designed to describe the vibration transmission characteristics of human body. The finite element model and multibody model are structurally and computationally complex. For the low degree of freedom lumped parameter model, the structure is too simplified and the fitting precision usually unsatisfactory. In order to establish a dynamic model with both simple structure and good fitting precision, a seven-DOF biodynamic lumpedparameter model with 23 parameters, which is referred to the literature, 1 was selected for dynamics analysis. Figure 1 is the schematic diagram of model. This model is represented by seven rigid bodies coupled with several dampers and springs. m 1 -m 7 are mass of pelvis, abdomen, diaphragm, thorax, torso, back, and head. c i and k i are spring and the damper of corresponding mass, x 0 represents the displacement excitation of the seat input, x 0 -x 7 represent the displacement responses of the corresponding mass units.
According to Newton's second law, the motion equations of human body system were obtained by force analysis of each mass unit, as shown in equation (1).
Equation (2) is the matrix form of equation (1) The frequency domain equation (7) of the vibration system can be obtained by applying the Fourier transform of equation (2), where Q(jw) f g= 1; jw f g½x 0 .
Àc 56 0 Àk 56 0 By rearranging the equation (7), take the matrix A = À w 2 M + jwC + K, then we can get another equation of motion expression Therefore, the transfer function of human body response is The seat to head transmissibility is defined as the ratio of the head displacement to the excitation point displacement. The apparent mass is ratio of force to acceleration of the contact surface between seat and human butt. According to the definitions, the equation of seat to head transmissibility and apparent mass including the transfer function can be obtained It is assumed that 73% of the body weight is applied on the seat for a seated human body, that is, the sum of the masses of each mass block in the lumped mass model should be 73% of the subject weight. 23 On the basis of dynamics analysis, the errors between the calculated values of the indicators (amplitude and phase) and the experiment data were treated as the objectives.
The calculated data of vibration response of human model should be close to the experiment data as far as possible. The objective function is shown in equation (13). Obviously, the parameter identification of the seven-DOF model is a five-objects optimization problem.
In equation (13), AM amplitude calculated represents the calculated amplitude of apparent mass, AM amplitude tested represents the experiment data of AM amplitude, the other three equations follow this pattern; m sum is 73% of tested human weight. n is the number of models DOF.
The experimental data used in this study refer to ISO 5982-2019. 23 The specific test data of the seated human body STHT and AM exposed to vertical vibration are listed in Table 1. These data are applicable to subjects with a total mass of 75 kg. Therefore, the total mass of the lumped parameter model is 75kg Ã 73% = 54:75kg. In this study, we take the requirement of the sum of each mass unit as an object variable.
In order to make the result of parameter identification more consistent with the actual situation, at the same time, to accelerate the convergence of optimization calculation process, it is necessary to constrain model parameters by boundary conditions. The constraint conditions of m 1 -m 7 refer to literature, 1 according to the percentage of each mass units in the total mass, the upper and lower limits of the constraint range are extended to 10% of the actual mass. The stiffness and damping value of each structure of human body ranges from 100 to 300000 N/m 39 and 500 to 4000 Ns/ m. 40 Therefore, the constraint conditions of each parameter are shown in equation (14).
A reference vector guided evolutionary algorithm (RVEA) RVEA algorithm can better balance the convergence and distribution performance in the process of population evolution, and a special concept ''angle-penalized distance'' was proposed. 36 This algorithm mainly includes three parts: offspring creation, reference vector-based environmental selection, and reference vector adaptation. The detailed process of RVEA was introduced in the following.
Step 1: Initialization. According to the problem to be solved, the population size N, the maximum number of iterations gene_ max and the dimension of reference vector num_ V were determined. Then, the initial population and unit reference vectors are generated. The method of reference vector generation refers to Das and Dennis's method. 35 After generating reference points evenly distributed in the objective space, the reference vector is drawn by connecting the origin to the reference point. All the reference vectors are normalized into unit vectors according to the method in reference. 36  Table 2 lists the reference vector number of 3-7 dimensional MOPs. Step 2: Through crossover (e.g. simulated binary crossover) and mutation (e.g. polynomial mutation) operation, the offspring population Q t with the same size as the parent population is produced.
Step 3: Selection operation. After combining the parent population with the offspring population, N individuals that meet the selection conditions are selected from the 2 N size population to enter the next iteration. The combined population is called S t , and the selected population called P_ t + 1 .
Step 4: If the current population meets the calculation termination requirement, the iteration is stopped. If not, go back to Step 2. Figure 3 is the flow chart of RVEA.
The uniqueness of RVEA is reflected in the operation of environment selection. First and foremost, RVEA uses reference vectors to divide the objective space into multiple subspaces and performs selection in each subspace separately. Next, all the objective variables are shifted with the minimum point-z Ã as the reference origin. z Ã is the minimum value of each objective of all individuals.
Then, the angles between the individual P i and all the reference vectors are calculated, as shown in equation (17). P i is associated with the reference vector corresponding to the smallest angle between them. Figure 4 shows the association operation between individual and reference vector schematically.
After the association operation, several situations can occur: a vector associates an individual, a vector associates several individuals, or a vector has no individuals associated with it. Angle-penalized distance (APD) is proposed for environment selection operation. The calculation equation is shown in equation (18), where, P u i, j À Á and f 0 i are used to maintain the distribution and convergence of the population   respectively. The calculation equations are shown in equations (19) and (20). g where, i = 1, 2, ., 2 N; 2 N represents the size of combined population S t ; N v is the size of reference vector set. M represents the number of objective variables; a is a coefficient related to distribution. t is the number of current iterations. According to the parameter setting, the reference vector is adaptively adjusted after a certain iteration interval, the equation can be written as where, v t + 1 is the adapted reference vector, v 0 is the initial uniformly distributed reference vector; and the 8 represents the Hadamard product. By intermittently adaptive adjustment of the reference vector, the distribution of the population can be improved under the condition of maintaining the convergence stability of the algorithm.

Model identification and analysis
According to the experiment data of vertical vibration characteristics of seated human body in section 2, the single-objective optimization algorithm -genetic algorithm (GA) and multi-objective optimization algorithms-MOEA/D, NSGA-II, NSGA-III, RVEA were used to identify the model parameters. When using the single-objective algorithm (GA), the sum of five normalized objective functions was regarded as the objective function.
The parameter settings of different algorithms used in this paper are summarized as follows.
(1) For all the algorithms, the population size is set to 210, and the maximum number of iterations is set to 1000. (2) For NSGA-3 and RVEA involving the use of reference points or reference vectors, objective variables for each dimension were evenly divided into six parts, and the final size of the reference set was 210 3 5. (3) Settings for evolution parameter: We selected simulated binary crossover and polynomial mutation in evolutionary step. The distribution index and crossover probability were set to 30 and 1 for crossover operation, and the distribution index and mutation probability were set to 20 and 1/23 for mutation operation, where 23 is the number of parameters to be identified in the model. (4) Some specific parameter settings for different algorithms were listed as following. For MOEA/D, the number of weight vectors in the neighborhood field was set to T = 20, the Tchebycheff method was selected in polymerization process. The parameter a in RVEA for balancing convergence and distribution is set to 2. (5) The different algorithms are run independently for 20 times; other parameter settings refer to corresponding references.
We normalized the objective variables calculated by each algorithm using the same scale, the normalized parallel coordinate plots are shown in Figure 5. After independently running 20 times for each algorithm, the first pareto front combined by all solutions can be regarded as the true pareto front of the fivedimensional many-objective optimization problem in this paper, as shown in Figure 5(a). Figure 5(b) is an approximate PF obtained by MOEA/D after 20 operations, Figure 5(c) to (f) are approximate PF obtained by NSGA-II, NSGA-III, RVEA, and GA, respectively. Objective values obtained by NSGA-II and GA range from 0 to 1.35 and 0-1.81 respectively. Compared with the true pareto front (range from 0 to 1), the convergences of these two algorithms are not very satisfactory. The approximate PF obtained by NSGA-III fails to maintain good distribution. The calculation results of MOEA/D and RVEA show good convergence and distribution simultaneously. Meanwhile, it is obviously that RVEA can make all five objective variables smaller at the same time.
Inverted generational distance (IGD) is an indicator to synthetically evaluate the distribution and convergence of the solution set, computed as where, PF is the true pareto front of MOPs with uniformly distributed points, PF Ã is the approximate non dominated solution set obtained by different algorithms, and PF j j is the solution number of PF set. In order to provide a uniformly distributed true pareto front, the normalized non dominant solution set obtained from the previous calculation was filtered with reference vectors. The generation process of reference vectors was described in the section 2. For the fivedimensional many-objective problem to be solved in this paper, each objective is divided into 11 parts, and the final number of reference vectors is 1365. Each nondominated solution was bound to its nearest vector. Finally, 1365 approximately uniformly distributed solutions were selected as the true PF for IGD calculation. Table 3 lists IGD value calculated by RVEA, MOEA/D, NSGA-II, NSGA-III, and GA respectively. The IGD mean value calculated by RVEA obviously lower than that calculated by other algorithms. The distribution and convergence of solution set calculated by RVEA are relatively satisfactory among the five algorithms used in this paper. When using the singleobjective algorithm (GA), the IGD value is the largest. This means that the convergence and distribution of the non-dominant solution set are unsatisfactory. This phenomenon can also be seen in Figure 5(f).

Result of identification
In order to evaluate the human model goodness-of-fit using the different algorithms, the parameter e was introduced 37,38 : where, t m À t c is the error between the experimental value and the calculated value, N is the size of the experimental data. Except for the fifth objective variable, the e coefficients of the other four objective variables are listed in Table 4. The fitting effect is accurate when the e close to 1. Based on the non-dominated solution set obtained by running the five algorithms independently for 20 times in the previous section, 2315 non-dominated solutions are obtained. We calculate the goodness of fit e for all the solutions, and average the e of the four objective functions, then sort them in order of the highest to lowest mean value. Finally, we get the best average goodness of fit is 0.8901. At this time, the goodness of fit of the AM amplitude and phase are 0.9456 and 0.9125, and the goodness of fit of the STHT amplitude and phase are 0.9193 and 0.7828 respectively. The parameter identification results of mass, stiffness, and damping corresponding to the best goodness of fit are given in Table 5. At this time, the mass sum of m 1 -m 7 is 56.8704 kg. The mass error between human model and experimental data is 2.1204, meets our modeling requirement.
The identification curves of AM amplitude and phase calculated by this set of parameters are shown in Figure 6. The calculated AM amplitude and phase can fit well with the experimental data. The frequency corresponding to the AM peak value of identified model is consistent with the measured values (4 Hz). The  Figure 7. The STHT amplitude and phase of the identified model are consistent with the change trend of the experimental data. Both the calculated and experimental resonant frequencies are 5 Hz corresponding to the STHT resonance peak. The calculated data of STHT phase and the corresponding test data show good agreements in the range of 0-10 Hz. The STHT phase calculated by the model in the range of 10-20 Hz is slightly larger than the experimental value. The RMS of STHT amplitude error is 0.0788, and the corresponding RMS of AM phase error is 6.3237°.

Conclusions
In this paper, a model describing the vertical vibration of the seated human body was established, and the parameter identification of multi-parameter complex model was studied. The motion equation of a seven-DOF biodynamic lumped-parameter model with 23 parameters was calculated. The seat to head transmissibility (STHT) and apparent mass (AM) of biodynamic responses of seated human body were derived. In order to obtain an accurate dynamic model, a constrained many-objective optimization problem with five objective was established. The parameters of the model were identified by solving the many-objective optimization problem (MaOPs). The target of identification is to make the amplitude and phase of the AM and STHT calculated by the model coincide with the experimental data as much as possible.
Next, we introduced an algorithm called reference vector guided evolutionary algorithm (RVEA) to solve the parameter optimization identification problem of human body model. In order to verify the quality of nondominated solution set calculated by RVEA, the other four optimization algorithms (NSGA-II, NSGA-III, MOEA/D, and GA) were selected for comparative analysis. The simulation results show that the RVEA obtained superior performance in inverted generational distance (IGD) and parallel coordinate chart. RVEA shows highly competitive performance on parameter identification of human body model based on manyobjective optimization problems.
Indicator e was introduced to evaluate the parameter identification accuracy of different non-dominated solution sets. All the non-dominated solutions were  sorted according to the indicator e, and the solution with the best fitting effect was selected. Finally, the identification results of AM and STHT, including amplitude and phase curves were shown and compared with the experimental data supplied by ISO 5982-2019.
The results show that the seven-DOF model is an efficient model to characterize the vibration transmissibility characteristic of the seated human body exposed to whole body vibration under a wide range of excitation conditions (0-20 Hz). Therefore, the model and parameter identification methods used in this paper show satisfactory effectiveness and accuracy on analyze ride comfort of occupant body biodynamics.

Data Availability
No data were used to support this study.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The presented work was funded by the National Key Research and Development Program of China (Grant no. 2018YFB0106200). The authors are grateful for the financial support.