Comprehensive Analysis and Optimization of Dynamic Vibration-Absorbing Structures for Electric Vehicles Driven by In-Wheel Motors

Distributed-drive electric vehicles (EVs) replace internal combustion engine with multiple motors, and the novel configuration results in new dynamic-related issues. This paper studies the coupling effects between the parameters and responses of dynamic vibration-absorbing structures (DVAS) for EVs driven by in-wheel motors (IWM). Firstly, a DVAS-based quarter suspension model is developed for distributed-drive EVs, from which nine parameters and five responses are selected for the coupling effect analysis. A two-stage global sensitivity analysis is then utilized to investigate the effect of each parameter on the responses. The control of the system is then converted into a multiobjective optimization problem with the defined system parameters being the optimization variables, and three dynamic limitations regarding both motor and suspension subsystems are taken as the constraints. A particle swarm optimization approach is then used to either improve ride comfort or mitigate IWM vibration, and two optimized parameter sets for these two objects are provided at last. Simulation results provide in-depth conclusions for the coupling effects between parameters and responses, as well as a guideline on how to design system parameters for contradictory objectives. It can be concluded that either passenger comfort or motor lifespan can be reduced up to 36% and 15% by properly changing the IWM suspension system parameters.


Introduction
Increasingly strict vehicle emission standards have been the main motivation for extensive research on electric vehicles (EVs), which are likely to replace the traditional internal combustion engine (ICE) in the near future [1][2][3][4]. From the perspective of system configuration, EVs can be categorized into two distinct types: centrally driven type and distributed-drive type. Compared with the centrally driven layout, the latter case incorporates an electric motor in the hub of a wheel for direct drive, and it provides numerous advantages, including more available chassis space, prompt system response and high energy efficiency [5][6][7][8][9]. In this context, distributed-drive EVs are attracting more attention from both academia and automobile industry and are also considered as the direction of current and future EVs.
However, the performance of the in-wheel motors (IWMs) used in distributed-drive EVs is restricted by the limited space inside the wheel and negative effect of the increased unsprung mass on suspension performance. Many noteworthy research projects have focused on these issues.
Current vibration mitigation algorithms used in the IWM suspension system can be categorized into four types.
(1) Motor dynamics optimization Some algorithms have been proposed to improve the noise, vibration and harshness (NVH) performance of the IWM. Takiguchi et al. [10] and Sun et al. [11] proposed current controllers to reduce the vibration caused by the IWM. (2) Suspension dynamics optimization Hredzak et al. [12] used optimization algorithms to improve the system performance by optimizing the suspension component parameters. The results achieved improved ride comfort and road handling. (3) Active suspension control Shao et al. [13] and Wang et al. [14] used controllable suspensions for IWMs, and these can provide better results than passive systems. (4) Novel structure design Other than the three abovementioned methods, the dynamic vibration absorbing system (DVAS) is also a possible candidate to mitigate vibrations caused by road excitation. The DVAS refers to a type of well-tuned spring-mass subsystem with which the main system vibration can be effectively reduced. It shows great potential for active vibration control, and a representative DVAS-based IWM suspension system is shown in Fig. 1 [15,16].
Among these four types, the DVAS-based method is considered as a promising direction, for the reasons that it could achieve 35% improvement of ride comfort and 30% enhancement of road handling in comparison with the traditional ones [15]. However, there exists little research on how to choose DVAS parameters for conflicting objectives, and the coupling effects between system parameters and outputs have yet to be comprehensively discussed.
To understand the relationships between system parameters and responses, local sensitivity analysis (LSA) is widely used. LSA is a one-at-a-time (OAT) technique: The effect of each parameter is evaluated by calculating the partial derivative of the cost function. This method is only valid at the central point where the calculation is carried out, and cannot easily reveal the nonlinear interactions of different factors [17]. In contrast, global sensitivity analysis (GSA) uses a representative set of samples with which the full design space is explored using Monte Carlo techniques. Compared with OAT algorithm, the overall outputs of GSA methods can model the interactions of different changes to dynamic characteristics of the system. This paper adopts a two-stage evaluation strategy to solve the widely recognized computational complexity problem associated with Monte Carlo-based algorithms.
The innovations described in this paper can be summarized as follows: (1) A novel GSA algorithm is proposed to investigate the coupling effects between parameters and responses for the DVAS-IWM system, which can avoid the issues existing in conventional sensitivity analysis methods; (2) PSO is employed to provide two sets of parameters for improving ride comfort or mitigating motor vibration. The combination offers more options for varying driving conditions. The rest of this paper is organized as follows. First, a DVAS-based IWM suspension system and road profile model are described in Sect. 2. Then, the GSA method is introduced and an in-depth analysis of the effect of system parameters is made. The PSO and the performance of the optimized system are presented in Sect. 4. A conclusion and future work are discussed in Sect. 5.

System Modeling
This section introduces the DVAS-based IWM suspension system and road profile generation.

The DVAS Model
The DVAS-based IWM suspension system has been well studied by previous researchers and can be regarded as a potential solution to the downsides of IWMs [15,16]. The structure of the tire-type DVAS-based IWM suspension system is listed in Fig. 2. Its dynamics can be described by Eq. (1), as in [16].
where x * ,̇x * andẍ * represent displacement, velocity and acceleration, respectively. The subscripts b, s 1 , s, r stand for the sprung mass, total unsprung mass, stator and axle mass and rotor mass, respectively. The system parameters are listed in Table 1 [16].

Road Profile Generation
The statistical features of a real-world stochastic road profile can be represented using power spectral density (PSD). Many methods have been proposed to generate random road profiles. This paper uses harmonic superposition algorithm to model the road profile in the time domain [18]: where q(t) is the generated road profile, f mid-K is the kth middle frequency, and G q f mid-K represents the PSD at f mid-K . Variable K is an identically distributed phase over the range of (0, 2π) . f 1 , f 2 stand for the upper and lower time-domain frequency boundaries, accordingly. ISO level C was adopted as the system excitation, and the velocity was set to 60 km/h.

Sensitivity Analysis
This section introduces the two-step GSA and analyzes the effects each parameter has on the system response characteristics.

Introduction to Global Sensitivity Analysis (GSA)
The flowchart of the GSA-based algorithm is illustrated in Fig. 3. After both the input and output parameters are chosen for the IWM suspension system, the elementary effects analysis is carried out. The influential parameters from the first step are then fed into the second step-the GSA. Consideration of the reliance on the sensitivity index, robustness, convergence and credibility assessment is required for the careful use of GSA. Kolmogorov-Smirnov statistic approach is used as a final check that the identified sensitivity factors are comprehensive.

Elementary Effects Analysis
The purpose of the first stage is to find out the appropriate parameter which has the biggest effects on the system responses. Note that the sensitivity index generally refers to a number determined by a predefined algorithm that shows the relative sensitivity of the results to the variables. This part uses elementary effects analysis to offer a good approximation to the sensitivity index. The functionality of the elementary effect (EE) analysis is to reduce the computational cost by finding the inputs that have no effect, which are then used to calculate the sensitivity index in the second step. For any system model y = f (X), X = 1, 2, ⋯ , k , it is assumed that the inputs are bounded in a k-dimensional unit cube Ω k with p levels. The elementary effect of any parameter X i is defined as: where e is a zero vector, except for its ith value, which is 1. This paper uses p = 4 and Δ = 2∕3 [17]. The average value of the ith EE is then denoted as , which is used to evaluate the importance of each parameter. A parameter with a higher value represents a greater effect on a specific output. Based on the above definition, both and the standard deviation of the ith EE can be calculated from Eq. (4), where R is the running time.
Note that for a complex model with interactive effects between parameters, results calculated based on values of may result in Type II errors, leading to the omission of influential parameters. The mean of , denoted as * , is thus considered as the final index for EE analysis.

Global Sensitivity Analysis with the Adoption of FAST Method
GSA has become a hot topic due to its independence from the central assumption of the rest of the parameters. The widely used GSA algorithms include the Fourier amplitude sensitivity test (FAST) [19], the Morris method [20], sampling-based methods [21] and the Sobol method [22]. Among these methods, FAST is applicable to complex nonmonotonic systems and offers superior computational efficiency.
For the model used in Sect. 3.1.1, the parameters in Ω k are given by: are the minimum and maximum values of X i . A search function is then introduced, with which each parameter can periodically oscillate inside the space Ω k at a frequency w i . Numerous search functions have been proposed; we adopted the search function proposed by Saltelli et al. [23]: where F −1 i is an inverse cumulative distribution function and s is a variable. It can be seen from Eq. (7) that the output is a periodic function of s. Such a function can be expanded with a Fourier series Y.
where the coefficients A k and B k are: and k is a finite integer. For an odd sample size N, k is calculated by k ∈ Z = {1, … , N − 1∕ 2} . If we denote the N samples as S = s 1 , s 2 , … , s j , … , s N , the corresponding model outputs are calculated as: where f s j = f X 1j , X 2j , … , X kj . The corresponding Fourier coefficients can therefore be expressed as where K = 1, … , (N − 1)∕2 . The variance of a system output V o is then calculated and decomposed into components at the integer frequencies [24]: is the spectrum and p satisfies pw i ⩽ (N − 1)∕2 . The result V i ∕V denotes the total variance of y.

Sensitivity Analysis Results
In this part, the simulation results of the adopted two-step GSA are discussed.

GSA Setting
According to Fig. 3, the first step in the GSA is to determine the candidate parameters and system responses. For the primary purpose is to investigate the effect of component parameters on an IWM suspension system, all the candidate parameters are listed in the second column of Table 2; the nominal values of these parameters are from the literature [15,16]. For the range of both the stator and rotor masses, which are related to the power requirements and restricted by the wheel hub dimensions, the variation is set within the range from + 20 to − 20%. The range of variation of the tire stiffness, which can be viewed as a function of the tire pressure, is also set within the range from + 20 to − 20%. The parameter variations for the other parameters are set to larger ranges, for their values may significantly change as controllable suspension components are adopted. Apart from the candidate parameter, the system responses also need to be determined; these responses, as well as the corresponding abbreviations, are tabulated in Table 3. In this paper, the first three indexes in Table 3 are widely used for the traditional suspension system evaluation, and the two latter indexes are for the IWM installed inside the wheel [16]. The reasons for these choices can be interpreted as follows: (1) SA corresponds to the vibration of the motor system, and a higher SA will result in premature fatigue failure and component damage [16], and (2) motor eccentricity represents the distance between motor and rotor, and an excessive Gap variation may result in contact with and damage to the stator body, as well as damage to the rotor.

Elementary Effects
As illustrated in Fig. 3, the first step in GSA is to determine the candidate parameters and system responses. Because our primary purpose is to find which parameters among all the candidate parameters have the biggest effect, the normalized EEs of all candidate parameters on all investigated responses are presented graphically in Fig. 4. The number marked in the center of each subsquare represents the EE of the investigated parameter; a smaller number means that a parameter has a bigger effect on corresponding response. Figure 4 presents the importance of each parameter to individual response through the color shade. Darker colors correspond to greater importance. Figure 4 demonstrates that the EE method can provide a qualitative description of the individual parameter effect.
Results shown in Fig. 4 reveal that the parameters of the DVAS structure have little EEs in the three traditional suspension-related indexes, which are SMA, RS and TD, and the parameters that have the biggest effect are k s , c s and k t , respectively. As for the two IWM responses, c s has the largest effect on SA, while the effect of suspension stiffness, DVAS stiffness and DVAS damping is also non-negotiable.  x s − x r m Motor eccentricity Gap

Fig. 4 Elementary effects of the investigated parameters on system responses
The parameter k b has a dominant effect on Gap, but the other parameters affect it much less.

GSA Results
This section discusses the GSA results calculated by FAST (graphically presented in Fig. 5), in which the normalized sensitive indices of all candidate parameters are demonstrated for SMA, RS, TD, SA and Gap. For SMA, it is obvious that k s and c s have larger indices than the other seven candidate parameters, among which k t has the biggest effect with an index of 0.024. This phenomenon can be understood as follows: Both k s and c s are directly connected to the sprung mass, and the vibration induced by an uneven road profile is attenuated by these two suspension components. In Fig. 5b, e, both indices of c s and k b are bigger than 0.7, which means these two parameters are the most important ones for these two responses under the investigation ranges given in Table 2. The situation with TD and SA is more complex as more parameters contribute to the variation in the response variances. In the SA result, c s is the parameter that has the biggest effect, with index values of k d and c d bigger than 0.1, which means that the stator mass vibration is closely related to those components that are connected to itself. The TD results in Fig. 5c indicate that k t plays an important role, while k d , k s and c s also have an impact on the range of the variation. Compared with the sensitivity analysis results for the traditional suspension structure in the study [25], in which k t and c s are the two parameters with the biggest effect on TD, it can be observed that the new parameter c d , which comes from the novel DVAS, can also affect the tire deflection. Based on this observation, tire vibration can be mitigated by varying the damping coefficient of the DVAS, which provides a new perspective for improving vehicle road handling capacity with the novel DVAS. Comparison of the EE analysis results in Fig. 4 and the GSA results in Fig. 5 reveal that the EE analysis can provide a good proxy to the sensitivity analysis results.

Parameter Optimization
This section introduces the DVAS-based system parameter optimization based on particle swarm optimization (PSO) and the effects that the parameter variations have on the investigated responses.

Optimization Algorithm
This section presents the parameter optimization algorithm and the effects that the parameter variations have on the investigated responses.
As the primary goal of suspension optimization is to balance system performance among multiple objectives under constraints, this paper first converts the DVAS-based system parameter design into a multiple-objective optimization problem (MOOP). For the investigated system, because a higher SA will result in premature fatigue issues as well as in component damage, we defined the objectives and the constraints of the MOOP as follows: where P ∈ ℜ is the parameters to be designed and (⋅) stands for the response RMS. In Eq. (13), f 1 (P) represents ride comfort and f 2 (P) corresponds to the vibration effects of the IWM. The three constraints are set to be the same as in [15]; i.e., lim (RS) = 0.1 m, lim (Gap) = 0.001 m, lim (TD) = 0.018 m. The purpose of the following optimization can be summarized as an effort to balance both ride comfort characteristics and stator acceleration by maintaining RS, Gap and TD within the constraints. This paper describes the global sensitivity of the parameters listed in Table 2, in which both the mass and component parameters of the suspension are considered to be variables. In the following parts, we will optimize the system according to the MOOP defined in Eq. (13). Because the masses are determined once the system configuration is fixed, we assume the mass is constant and only the suspension components, that is, k d , c d , k s , c s , c b , k b and the ranges defined in Table 2, are to be optimized by PSO. The PSO approach was proposed in the early 1990s [26]. It originated from animal social behavior, as a stylized representation of the movement of flocks of birds or schools of fish. The candidate solution, which is called a particle, moves inside the search space to obtain the best solution.

Optimization Results
For the MOOP defined in Eq. (13), both SMA and SA are regarded as the responses to be optimized. In this part, we provide two-parameter sets for improving ride comfort (SMA-oriented) or mitigating IWM vibrations (SA-oriented). The PSO results for the six parameters are given in Table 4.
It can be seen from Table 4 that all six parameters for both cases are within the range defined in the last column. Based on the two-parameter sets, the simulation results for a vehicle driving on a level but nonsmooth level C road at 60 km/h are given in Figs. 6, 7 and Table 5. The passive responses in the plots and tables correspond to the IWM suspension system without DVAS, and the parameters are the same as in [15]. The RMS comparisons (given in Table 5) show that the emphasis in the design of the IWM suspension system characteristics can be changed by varying the weights of the objective functions. The variations in the six parameters result in a 15% change in SMA (from 0.544 to 0.627 m/s 2 ) and a 36% change in SA (from 3.67 to 5.02 m/s 2 ). Another observation from Table 5 is that SMA and SA are conflicting objectives: They cannot be improved simultaneously.   The comparison in the frequency domain is shown in Fig. 7. The comparison results for SMA show that the SMAoriented parameter set can improve ride comfort, especially in 4-8 Hz frequency range, which is the most sensitive frequency range for the human body. The SA comparison results show that the performance of the SA-oriented dataset is superior to that of the others, especially in the 8-12 Hz range.
Based on the above simulation results, two conclusions can be drawn for parameter tuning. (1) The most important parameters for SMA are the vehicle suspension parameters, that is, k s and c s . Smaller k s and c s are beneficial to ride comfort improvement. (2) Tuning on an SA-oriented system is more complex than on an SMA-oriented system, as multiple parameters have a non-negligible effect on SA. A bigger c s value and smaller k d and c d values will be helpful for SA reduction.
In summary, it can be concluded that the changes of IWM suspension system parameters can achieve different system characteristic objectives.

Conclusions
This paper presents theoretical and numerical simulation studies of the coupling effects between IWM and suspension system parameters and characterizes responses to parameter changes in IWM-driven EVs. Nine variable system parameters and five dynamic-related responses were selected based on the DVAS-based IWM suspension model, and a two-stage GSA-based algorithm was used to investigate the effect that each parameter had on the responses. In this algorithm, the EE of each individual parameter was first estimated to reduce the calculation complexity. GSA was then used to calculate the sensitivity indices. It was found that the damping between the rotor and stator was effect-free, while the other parameters had prominent influences on system responses. The results of this analysis showed that EE analysis could provide a good proxy for the overall sensitivity analysis result and a guideline on how to change system parameters for desired system responses. Finally, PSO was adopted to optimize the system for the engineer's choice of either improved ride comfort or mitigation of IWM vibrations. Numerical simulation results showed that SMA and SA are conflicting objectives, and parameter changes can improve SMA and SA by up to 15% and 36%, respectively. Future research will focus on: (1) Investigation of the parameter sensitivities for a full vehicle model, considering longitudinal and lateral dynamics. (2) Research on adaptive DVAS-based IWM suspension systems based on the optimization results.