Chapter 1 Introduction
1.3 Motives and objectives
ATVs are usually driven on rough terrains. In order to maintain the stability and improve the driving comfort, an effective method for designing suspension parameters is needed. A frame of optimization design was developed in this paper.
For passive suspension systems shown in Figure 1.2, the parameter such as spring coefficient k and damping coefficient c are untunable during driving. So the presetting of k and c affects the driving comfort significantly. Based on optimization algorithm, the optimized k and c
are designed suitable for certain terrain. Therefore, different terrains have corresponding different settings of k and c to improve driving comfort.
Besides, the rattle-space limit is also concerned in this paper to insure the travel of suspension is within its safety limit.
Chapter 2 Dynamic model of ATV with passive suspension system
2.1 Derivations of equations of motion
The model we consider here has three degree-of-freedom including vertical displacement of the center of gravity, roll and pitch angles. A full-car model of ATV is illustrated in Figure 2.1. A and D are the front two passive suspensions, C and D are the rear two passive suspensions, and O is the center of gravity. ZJKA
, ZJKB , ZJKC
, ZJKD
are the road profile inputs to the four suspensions, respectively. ZJK
is the vertical displacement along the z-axis. When vehicle rotated in the inertial reference, the included angles are roll, pitch, and yaw angles, with notationsα, β, and γ respectively, called the Euler angles. The signs of rotate angles are according to right-hand rule. Here we assume γ is very small and can be neglected when driving alone a straight lane. Therefore, the rotation matrix of the vehicle can be represented as follow
1
cos sin sin cos sin
( , , ) ( , ) ( , ) 0 cos sin
sin sin cos cos cos
v x y z v x y x y
pitch axis:
With the rotation matrix Rv, the rotated vectors from O to each suspension can be represented as follows
'
JJJK JJJK JK JJK JK JJK
A
JJJK JJJK JK JJK JK JJK
A
JJJK JJJK JK JJK JK JJK
A
JJJK JJJK JK JJK JK JJK
A
(2-3)
Here ZK0
is the origin point of local reference. Then we can derive the vertical displacement of each suspension, and for simplicity, we applied linearization by assuming α and β were small.
'
To get the net vertical displacement of each suspension, we subtract the road input form vertical displacement of each suspension.
0 0
0 0 With the net vertical displacement of each suspension and its first derivative, we can derive the equations of motion by applying Lagrangian Dynamics. where T represents the kinetic energy and V represents the potential energy of the system. For simplicity, we assume
A D f
Here k is the spring coefficient and c is the damping coefficient of each suspension.
Then we have and the generalized forces can be derived as follows.
( ) ( ) ( ) ( )
Substituting (2-9) and (2-10) into the Lagrangian equation as follows.
d ( L) L Z
We obtained the three equations of motion of the ATV.
( 2 ) ( 2 )
[ ( 2 ) ( 2 )]
The above equations of motion (2-12), (2-13), (2-14) can be rearrange as follows
and the state-space form of the system is written as follows
1 1 1
Equation (2-16) can be used to construct the dynamic model with Simulink. The scheme of Simulink model is shown in Figure 2.2
2.2 Whole body vibration
Whole body vibration is the mechanical vibration (or shock) transmitted to the body as a whole. It is often due to the vibration of a surface supporting the body [11]. The effects of vibration could be subdivided into three main headings, namely (a) interference with comfort, (b) interference with activities and (c) interference with health. Each criterion has different conditions and limits associated with it.
Frequencies below about 0.5Hz may eventually lead to symptoms of motion sickness. Different parts of the body resonate at different frequencies. In the vertical direction resonance starts at about 2Hz, but the first major resonance occurs at about 5Hz. The transmissibility of vertical vibration to the head is sometimes a maximum at 4Hz. The voice may be caused to warble by vibrations between 10Hz and 20Hz. Vision may be affected at any frequency, but blurring occurs between 15Hz and 60Hz. The dominate vibration transmitted through the seats of vehicles is often at frequency below 20Hz [11]. It was thus decided to limit the investigation to frequencies below 20Hz.
Comfort is more difficult to determine since it is subjective and has a lot to do with the psyche. It is important to note that frequency and not only amplitude contribute to the perception of comfort. An experiment done by Fothegrill and Griffin in 1977 showed that for a seated person excited by a 10Hz sinusoidal vibration for RMS values below 0.4 m/s2 it was noticeable,
but not uncomfortable. A level of 1.1 m/s2 was considered to be mildly uncomfortable, 1.8 m/s2 uncomfortable and RMS values above 2.7 m/s2 as very uncomfortable.
Most researchers agree that seated humans have a vertical vibration natural frequency at about 4 to 6Hz and a horizontal vibration natural frequency at about 1 to 2Hz. In this frequency ranges the seat motion is most easily transmitted to the upper parts of the body and is not just confined to the area of the body close to the source of vibration [12].
2.3 Legislation
WBV exposures are to be determined separately for the three axes in accordance with ISO 2631-1:1997. Gallais found that comfort is time dependent and suggested that time dependency curves should be incorporated in the evaluation of WBV. ISO 2631 of 1997 doesn’t compensate for tome dependency of comfort. It was noted that such dependency curves were included in the ISO of 1974, but was since removed.
Up to date there has been much debate about whether the VDV measurements should be used to evaluate WBV on seated humans. The arguments for and against the VDV measurements are usually that it is better at indicating the presence of peaks which tend to be more harmful to humans, but is also more sensitive to human induced motion such as shuffling [13].
VDV is based on the fourth power of acceleration and is therefore more sensitive to shocks compared to the RMS (root mean square) magnitude (ISO 2631-1:1997). The general formula for VDV is
1/ 4 weighted acceleration data. This parameter is time dependent and gives an objective measure of the amount of vibrations a person had to experience with a certain period.
VDV is used as an index of vertical vibration. Besides, the RMS angular velocity is used in this paper as the index of rotation vibration. The general form for RMS angular velocity is
1/ 2 weighted acceleration data. The reason of frequency weighting is that human body is more sensitive to certain vibration frequency and is described as follows.
The perception and influence of whole body vibration in humans is strongly dependent on the direction, magnitude and frequency of the vibration exposure. From the literature review it was seen that seated humans have a vertical resonance at about 4-6Hz and a horizontal resonance at about 1-2Hz [11]. At these frequencies the motion of the seat
is transferred most easily to the upper parts of the body. Frequency weightings are applied to the acceleration data. The frequency weightings are designed to not affect those frequencies where the body is most sensitive and to attenuate at those frequencies where the response of the body is less sensitive. In principle, weightings do not amplify at any frequency. Although there are problems with the use of the frequency weightings, there are currently no alternative methods of assessment.
A couple of weightings curves (filters) are specified by ISO 2631:1997 depending on the orientation of the person and the direction of vibration.
For a seated person the Wk curve is used to weigh the frequency contribution for vertical vibrations and the Wd curve is used to weigh the frequency contribution for lateral vibrations [12]. The frequency weighting curves are shown in Figure 2.3 [16]. The black line represents the Wk
weighting curve and the blue one represents the Wd weighting curve. The lateral axis is frequency represented in log scale and the vertical axis is magnitude represented in log scale, too. The idea of the frequency weighting curve is to attenuate the magnitude of acceleration (or velocity) whose frequency lies outside the human resonance frequency. For vertical vibration, the frequency weighting curve is Wk which attenuates the magnitude whose frequency doesn’t lie between 4 and 6 Hz. For lateral vibration, the frequency weighting curve is Wd which attenuates the magnitude whose frequency doesn’t lie between 1 and 2 Hz.
Assessments are made independently in each direction. Figure 2.4 shows which weighting curve should be applied to which axis [16]. For lateral vibration, the weighting curve is Wd and for vertical vibration, the weighting curve is Wk. A second order shaped curve of the form
2 has been used in [13, 14] to approximate the ISO weighting curve Wk and shown in Figure 2.5. Another second order shaped curve of the form
2 has been used to approximate the ISO weighting curve Wd and shown in Figure 2.6.
Chapter 3 Suspension optimization
3.1 History of Genetic algorithm
Computer simulations of evolution started in 1954 with the work of Nils Aall Barricelli. From this beginning, computer simulation of evolution by biologist became more common in the early 1960s. Artificial evolution became a wildly recognized optimization method as a result of the work of Ingo Rechenberg in the early 1970 – his group was able to solve complex engineering problems through evolution strategies. GA in particular became popular through the work of John Holland, especially his book in 1975. Research in GA remained largely theoretical until the mid-1980s, when the first International Conference on Genetic Algorithm was held at the University of IIIinois.
Genetic Algorithm (GA) or Evolutionary Algorithm (EA) is applied to the optimization of the suspension parameters in this paper. GA is computer-based techniques that mirror natural genetic evolution, and they have been found to be successful in application to a wide range of problems that are difficult to solve analytically. GA is categorized as global search heuristics, and use techniques inspired by evolutionary biology such us inheritance, selection, crossover, and mutation. Comparing with the gradient-based optimization methods, GA needs no derivative information and can deal with a large number of parameters. Besides, GA searches from a wide sampling of the cost surface simultaneously. Figure 3.1 shows the gradient-based optimization algorithm and Figure 3.2 shows the GA-based
The evolution search usually starts from a population of randomly generated individuals and happens in generations. In each generation, the fitness of every individual in the population is evaluated, multiple individuals are stochastically selected as parents from the current population (based on their fitness) to proceed crossover, and some of the offspring are randomly mutated. Then the selected parents and offspring construct the new population used in the next iteration of algorithm. The scheme of the procedure of GA is shown in Figure 3.3. [8].
In this paper, the programming of GA is written using MATLAB and combined with the ATV model build with simulink. Each part of the program will be introduced in next paragraph briefly.
3.2 GA procedures
The genetic algorithm program contain following parts: initial population, cost evaluation, mate selection, crossover and mutation [16].
Each part of above will be stated as follows.
3.2.1 Initial population
The initial population is generated randomly between the upper and lower bounds of parameters.
( ) { pop, par}
IPOP= hi−lo ×random N N +lo (3-1) where random N{ ipop,Npar} is a function that generates an Npop
x N
par matrix of uniform random numbers between zero and one, hi is the highest value in the parameter range and lo is the lowest value in the parameter range.In our problem, each individual in the population contains two parameters, spring coefficient k and damping coefficient c, which values are chosen according the constraint functions. The size of initial population
N
pop and numbers of generations will both affect the time cost of computation. The recommend size of initial population and numbers of generations are 40 and 100.3.2.2 Cost evaluation
Mathematical optimization is the process of the formulation and then the solution of a constrained optimization problem of the general mathematical form: called the design variables, f(x) the objective function or cost function and
g
j(x) and h
j(x) are respectively referred to as the inequality and equality
constraint functions.The formulation of the optimization problem requires identifying the cost function, the variables and the constraints. In our optimization problem, the variables are the spring coefficient k and damping coefficient c. The constraints are inequality functions which limit the upper and lower bounds of k and c.
chapter, and the other is the rattle-space limit. A kind of fitness measure, ϕ, that penalize controls that pass too close the rattle-space limits [15]. The penalty function is shown in Figure 3.4 and written of the form
( 1) parameters which are used to penalize a suspension that is too close to its rattle space limits. At a distance of m1+m2 it has reached its safety limit;
within a distance of m1 it is inside the safe travel limit. The summation of ϕ, equation (3-4), is used to indicate the extent to which the system stays within the limits of the rattle space during a simulation.
1 Thus, the two fitness measure used in the optimization was a weighted sum of the measures in equation (2-20) and (3-3), as shown in equation (3-5) below. The change in outcomes can be compared as the weighting, λ, is varied from zero to one.
1/ 4 1/ 2 1/ 2
Fitness of each individual of the initial population is calculated using the cost function, equation (3-5) and is sorted in descending order. That is,
3.2.3 Mate selection
The selection of mate is according to the fitness of each individual. The better the “gene” is, the higher the possibility of be preserved. The mate selection method used in this paper is the roulette wheel selection. Here we defined the xfitness which equals fitness of each individual divided by the sum of all fitness of the population. Therefore, the comparison table can be created by the cumulative probability. The MATLAB code is as follows
xfitness = fitness/sum(finess); % Roulette wheel selection
cum_probability = cumsum(xfitness); % Comparison table created
Selection of mate is to create a random value varied form one to zero and choose the first individual whose xfitness is larger than the random value. Repeating the selection procedure, the chosen parents will proceed with the crossover.3.2.4 Crossover and Mutation
Parents selected above reproduce offspring through the crossover procedure. Kinds of method for crossover were developed and the single point crossover is used here because of the numbers of parameters. The single point crossover means the parameters before crossover point exchanged and the parameters after crossover point mixed using following equation. where μ is varied form zero to one. The crossover procedure is illustrated
and 8 % often works well. Number of mutate parameters is decided through equation (3-7)
_ * *
mut pop par
N =mutate rate N N (3-7) where Nmut
is the number of mutate parameters. The mutation is, for
example, if Nmut equals two then two parameters of the offspring are replaced by new values which are generated randomly between the parameters constraints.By producing offspring using the above methods of crossover and mutation, a new solution is created which typically shares many of the characteristics of its "parents".New parents are selected for each offspring, and the process continues until a new population of solutions of appropriate size is generated. These processes ultimately result in the next generation population of parameters that is different from the initial generation.
Generally the average fitness will have increased by this procedure for the population, since only the best individuals from the first generation are selected for mating. This generational process is repeated until a termination condition has been reached. Common terminating condition is when a fixed number of generations reached.
3.3 Optimization of passive suspension
The passive suspension systems consist of the passive components spring and oil damper. When the vehicle bumped along the rough road, the work of suspension system is to absorb energy and reduce the vibration.
Since the spring coefficient k and damping coefficient c of passive
c is important for the comfortable and vehicle handling. In order to find the appropriate value of k and c for certain road condition, our scheme is showed in Figure 3.6. For passive suspension system, we first collect certain road profiles as the inputs ZJKA
, ZJKB , ZJKC
, ZJKD
. Then we collect the acceleration data of ATV and through the optimization procedure introduced above, the optimized k and c will be suit to this road condition. Different road profiles result in corresponding different setting of k and c. We call this procedure an offline design method.
Chapter 4 Simulation result
4.1 Road profile data
The road profile data is collected from an ATV driving simulation game made by IMON corp. The game contains kinds of different road conditions including loess, timberland, highway, sand beach and jouncing areas. When player driving ATV on these areas, the setting of suspension parameters k and c will affect the feeling of the whole body vibration transferred from the road. In order to improve the comfort of driver for different road conditions, we use the road profile data as inputs ZJKA
, ZJKB , ZC
JK , ZJKD
, then proceed the optimization procedure to find the optimized values of k and c.
Two kinds of road profiles were used here and the left and right road profiles which pass through A, B and C, D are similar but not the same.
Road profile one is jouncing road with a 35 centimeters high-low variation and higher bump frequency. The number of data points is 3180, and sampling period is 0.005 seconds. The graph is shown in Figure 4.1 and Figure 4.2 which represent the road profiles pass through left and right of ATV, respectively. Road profile two is mild road with a 10 centimeters high-low variation and lower bumping frequency. The number of data points is 2900, and sampling period is 0.005 seconds. The graph is shown in Figure 4.3 and Figure 4.4 which similarly represent left and right road profiles, respectively.
4.2 Optimization results
The ATV model used here is HONDA TRX250EX Sport ATV with specifications listed in Table 4.1. The suspension system is made by DNM Suspension with adjustable spring and damper design for ATV using. The default setting of suspension parameters is k=22000N m/ and
1200 /
c= N s m. Before optimization, the setting of genetic algorithm is listed in Table 4.2 and the limit constraints of k and c are
10000< <k 25000 and 1000< <c 2500.
Simulation one is using road profile one shown in Figure 4.1 and Figure 4.2 as the inputs from ground to wheels A, B and C, D. The optimized k=14480.6587N m/ and c=1002.6731N m/ . The optimized weighted vertical acceleration z is shown in Figure 4.5 with VDVz attenuated from 16.8403 to 13.6141m s/ 1.75. The blue line is the response of vertical acceleration with optimized k and c, the red line is the response with k and c before optimized. The zoom in of Figure 4.5 is shown in Figure 4.6. In Figure 4.6, the attenuation of weighted vertical acceleration is obviously.
The optimized weighted roll velocity α is shown in Figure 4.7 with Jα attenuated from 0.5442 to 0.5341rad s/ 0.5. The zoom in of Figure 4.7 is shown in Figure 4.8. The maximum value of weighted roll velocity shown in Figure 4.8 attenuate from 0.55 to 0.5 rad s/ 0.5.
The optimized weighted pitch velocity β is shown in Figure 4.9 with Jβ attenuated from 0.3838 to 0.3765rad s/ 0.5. The zoom in of Figure 4.9 is shown in Figure 4.10. The maximum value of weighted pitch velocity shown in Figure 4.10 attenuate from 0.31 to 0.28 rad s/ 0.5.
The fitness of genetic algorithm is shown in Figure 4.11. The lateral axis is generation and the vertical .the best generation occurred at number 90.
Simulation two is using road profile two shown in Figure 4.3 and Figure 4.4 as the inputs from ground to wheels A, B and C, D. The optimized k =12503.0626N m/ and c=1009.9273N m/ . The optimized weighted vertical acceleration z is shown in Figure 4.12 with VDVz attenuated from 4.28 to 3.4457m s/ 1.75. The blue line is the response of
Simulation two is using road profile two shown in Figure 4.3 and Figure 4.4 as the inputs from ground to wheels A, B and C, D. The optimized k =12503.0626N m/ and c=1009.9273N m/ . The optimized weighted vertical acceleration z is shown in Figure 4.12 with VDVz attenuated from 4.28 to 3.4457m s/ 1.75. The blue line is the response of