Optimization of nonlinear quarter car suspension-seat-driver model.

In this paper a nonlinear quarter car suspension-seat-driver model was implemented for optimum design. A nonlinear quarter car model comprising of quadratic tyre stiffness and cubic stiffness in suspension spring, frame, and seat cushion with 4 degrees of freedom (DoF) driver model was presented for optimization and analysis. Suspension system was aimed to optimize the comfort and health criterion comprising of Vibration Dose Value (VDV) at head, frequency weighted RMS head acceleration, crest factor, amplitude ratio of head RMS acceleration to seat RMS acceleration and amplitude ratio of upper torso RMS acceleration to seat RMS acceleration along with stability criterion comprising of suspension space deflection and dynamic tyre force. ISO 2631-1 standard was adopted to assess ride and health criterions. Suspension spring stiffness and damping and seat cushion stiffness and damping are the design variables. Non-dominated Sort Genetic Algorithm (NSGA-II) and Multi-Objective Particle Swarm Optimization - Crowding Distance (MOPSO-CD) algorithm are implemented for optimization. Simulation result shows that optimum design improves ride comfort and health criterion over classical design variables.

system matrix

amplitude ratio of head RMS acceleration to seat RMS acceleration

amplitude ratio of upper torso RMS acceleration to seat RMS acceleration

frequency weighted RMS head acceleration (m/s2)

frequency weighted RMS sprung mass acceleration (m/s2)

frequency weighted head acceleration (m/s2)

damping coefficient (N s/m)

lumber spine damping (N s/m)

thoracic spine damping (Ns/m)

cervical spine damping (N s/m)

objective function

stiffness (N/m)

lumber spine stiffness (N/m)

thoracic spine stiffness (N/m)

cervical spine stiffness (N/m)

nonlinear spring stiffness (N/m3)

tyre stiffness (N/m)

nonlinear tyre stiffness (N/m2)

mass (kg)

vibration dose value at head (m/s1.75)

road profile (m)

displacement (m), velocity (m/s) and acceleration (m/s2)

sprung

unsprung

frame

seat cushion

thigh and pelvis

lower torso

upper torso

head

Introduction

Suspension system along with seat has been widely used in vehicles to isolate passengers from shock and vibrations arising due to road unevenness. Seat-suspension system thus provides ride comfort, reduces fatigue during driving conditions and improves health and safety of drivers. Other performance requirements for a suspension system are to adequately support the vehicle weight, to maintain the wheels in the appropriate position so as to have a better handling and to keep the tyre in contact with the ground. The passive suspension systems are the trade-off between ride comfort and handling [1]. Due to the conflicting requirements, the suspension system has been investigated by many researchers to find the optimal trade-off amongst the conflicting requirements.

Gobbi and Mastinu [2] presented Multi-Objective Programming and Monotonicity analysis based optimization method for finding the trade-off for conflicting performance requirements such as discomfort, road holding and working space. A 2 DoF quarter car model running on random road profile was used. The optimal settings of the vehicle suspension parameters such as tyre stiffness, spring stiffness, and damping were derived either symbolically and/or numerically. Verros and Natsiavas [3] presented optimization of suspension stiffness and damping. A quarter car model travelling on a random road profile was used for optimization study. Authors had used and presented a critical comparison of quarter car models with passive linear and dual-rate suspension dampers and semi-active sky-hook damping models. Optimization of a light commercial vehicle to improve vehicle ride and handling was performed by Özcan et al. [4] using a quarter car and the half car models in Matlab/Simulink® environment. The performance criterions considered were RMS body acceleration, tyre forces, and body roll. The performance of the optimized suspension unit was verified using Carmaker model.

Molina-Cristobal et al. [5] had presented multi-objective optimization of a passive suspension system using quarter car model using meta-heuristic optimization with the multi-objective genetic algorithm (MOGA) and bilinear matrix inequalities (BMI) techniques. Ride comfort using RMS body vertical acceleration and road holding criterions were used as objective functions during optimization. Chi et al. [6] had presented optimization of linear quarter car model using three different techniques namely genetic algorithm (GA), pattern search algorithm (PSA) and sequential quadratic program (SQP) subjected to body acceleration, suspension working space, and dynamic tyre load as design criterions. Gomes [7] presented optimization of 2-DoF quarter car model travelling over a random road surface. The particle swarm optimization (PSO) algorithm is used for optimization. Minimization of dynamic vehicle load and minimization of suspension deflection were used as objective functions in two optimization examples. Baumal et al. [8] presented GA-based optimization of half car model with an objective to minimize acceleration of the passenger’s seat, subject to constraints such as road holding and suspension working space.

Aforesaid literature had studied suspension system as the most significant factor and optimized the suspension parameters, spring stiffness and damping, for ride comfort, suspension space, road holding, and dynamic tyre force as objective functions.

Kuznetsov et al. [9] had presented optimization of quarter car model coupled with a driver. A 3-DoF driver–car model, a quarter car having 2 DoF and a driver having 1 DoF, is developed for optimization. Ride comfort criteria as per ISO 2631-1 were used for optimization using the algorithm for global optimization problems (AGOP). Gundogdu [10] presented optimization of quarter car suspension system with seat using GA. A two DoF linear quarter car model was developed including 2-DoF lumped mass driver model. Objective function is formulated using head acceleration, crest factor, suspension deflection and tyre deflection objective functions. Objective functions are converted into uni-objective function using non-dimensional expressions, giving equal importance to each of the objective functions. Badran et al. [11] presented optimization of a human-car suspension system using GA. Quarter car model is used as a vehicle model. Seat acceleration, head acceleration, and suspension working space were used as the optimization criterion. The objective function was converted into an uni-objective function using weighting parameters. Results are compared to step and sinusoidal road profile. Thus, it is observed that a quarter car suspension system along with driver model was optimized using uni-objective function although the optimization problem is of multi-objective nature [10], [11]. Also, the human bio-mechanical model considered for the study is either 1 DoF [9] or 2 DoF [10].

Patil and Palanichamy [12] showed that body parts such as head, torso, and pelvis respond to a much greater extent to road induced vibrations as compared to seat response. Hence, car suspension system should be aimed to optimize considering responses of other body parts such as head, torso, thorax, abdomen, and diaphragm. This study presents the multi-objective optimization of nonlinear quarter car model coupled with a 4-DoF driver model. A nonlinear quarter car model having quadratic tyre stiffness and cubic stiffness in suspension spring is modelled. Multi-objective optimization is presented using two optimization algorithms – Non-dominated Sort Genetic Algorithm II (NSGA-II) and Multi-Objective Particle Swarm Optimization with Crowding Distance (MOPSO-CD). Vibration Dose Value (VDV) at head, Frequency weighted Root Mean Square (hereafter called as RMS) head acceleration, amplitude ratio of head RMS acceleration to seat RMS acceleration, amplitude ratio of upper torso RMS acceleration to seat RMS acceleration, crest factor, suspension space deflection and dynamic tyre force in terms of dynamic tyre deflection are considered as objective functions for optimization. Results are presented in tabular as well as graphical format.

Methodology

Mathematical modelling – nonlinear quarter car suspension–seat–driver model

To study the behaviour of a dynamic system and to optimize the same, mathematical model of the system is required. Modelling of a suspension system is done in the vertical plane. Longitudinal or transverse deflections of the suspension components are considered negligible in comparison with vertical deflections. The model is simple, yet consists of basic elements of the suspension system such as sprung mass ms (representing chassis) and unsprung mass mus (representing wheel assembly and axle). Springs and dampers, representing suspension element, are connected between the sprung and unsprung masses and spring, representing tyre, is connected between unsprung mass and ground respectively.

A suspension system of commercial vehicles generally consists of coil springs. During the mathematical modelling of a suspension system, the elements of a suspension system, springs and dampers, are considered as linear. But the spring exhibits nonlinear nature. Also, the stiffness of pneumatic tyre is nonlinear in nature. Hence, while mathematical modelling of a suspension system these nonlinear elements should be considered. McGee et al. [13] presented a detailed study on nonlinearities in the suspension system. Frequency domain technique for characterizing the nonlinearities is presented and validated by laboratory shaker with road data for validation. It was concluded that suspension system has quadratic and cubic stiffness nonlinearities and Coulomb friction. Zhu and Ishitobi [14] presented chaotic response of 7 DoF model subjected to nonlinear tyre stiffness and nonlinear suspension spring stiffness. Lixia and Wanxiang [15] presented bifurcation and chaotic response of 2 DoF vehicle model having nonlinearities in tyre stiffness and suspension spring stiffness.

In the present analysis, a nonlinear quarter car model having quadratic tyre stiffness and cubic stiffness in suspension spring as nonlinearities is considered along with seat suspension model consisting of a frame and cushion is shown in Fig. 1a.

Fig. 1a

Quarter car–seat–suspension – human model.

Human body is very complex and sophisticated dynamic system. In the literature, many mechanical models have been developed based on lumped-parameter models. Coermann [16] developed 1-DoF model consisting of single second order differential equation. Wei and Griffin [17] developed 1 and 2 DoF linear models with an assumption that human body is seated firmly on the seat. In the early studies, 3 DoF seated human subject was modelled by Suggs et al. [18] as damped spring-mass model. Wan and Schimmels [19] and Boileau and Rakheja [20] modelled a 4 DoF lumped parameter model.

In this study a 4 DoF lumped parameter human model suggested by Boileau and Rakheja [20] was used in optimization study. It consists of head and neck mass (mh), chest and upper torso mass (mut), lower torso mass (mlt) and thigh and pelvis mass (mt). A 4-DoF human bio-mechanical model developed by Boileau and Rakheja [20] considers typical driving conditions such as seated posture with feet support and hands held in driving conditions. During model development, magnitude and phase characteristics of driving point mechanical impedance (DPMI) and seat-to-head transmissibility (STHT) are satisfied. As DMPI signifies the dynamic load at the point of input whereas STHT signifies dynamic behaviour of body parts, DMPI of seated subjects was determined from laboratory experiments whereas STHT values are determined from published data. To analyse whole body vibrations, the bio-dynamic human model parameters are estimated by simultaneously optimizing the magnitude and phase responses of DMPI and STHT, under driving conditions [20].

Arslan [21] conducted experimental study on three different bio-dynamic models subjected to three different road types to provide quantified assessment. To assess the bio-dynamic models, experimental data on seated subject are recorded and compared with simulated data of models. Assessment of models is based on root mean square difference (RMSD) and Pearson correlation coefficient (PCC) between simulated and experimental results. According to results, the model suggested by Boileau and Rakheja [20] showed the best correlation with experimental data.

According to D’Alembert’s principle, the governing equations of motion representing nonlinear quarter car suspension–seat–human model are:

The nonlinear quarter car seat–suspension–driver model parameters are as follows:The nonlinear and linear quarter car–seat–suspension model along with human model is simulated in Matlab/Simulink®. A nonlinear quarter car model, having quadratic tyre stiffness nonlinearity in tyre and cubic stiffness in suspension spring, and linear quarter car model are modelled and simulated. During simulation of linear quarter car model, nonlinear parameters such as ktnl and ksnl, are kept zero. Both models are simulated using input as step (step size of 0.1 unit), bump (bump height of 0.1 unit) and class C road.

From step response, it is observed that VDV at head and RMS head acceleration of nonlinear model is 5.8180 m/s1.75 and 1.5636 m/s2 respectively whereas for linear model it is 3.7402 m/s1.75 and 1.0369 m/s2. The VDV at head and RMS head acceleration of the nonlinear system is greater than linear system due nonlinearities present in tyre and suspension spring. Maximum head acceleration of a nonlinear system is higher than the linear system. The increase in head acceleration of nonlinear system than that of the linear system can also be observed in frequency response plots. The crest factor, amplitude ratio AR_h and AR_ut are higher for a nonlinear model as compared to the linear model. As AR_h and AR_ut ratios are higher for a nonlinear system hence more magnitude of accelerations will be transferred from the seat. RMS Sprung mass acceleration of a nonlinear system is also higher than the linear system due to quadratic nonlinearity in tyre and cubic nonlinearity in suspension spring. Refer Fig. 1b and Table 1.

Fig. 1b

Step input – linear and nonlinear quarter model – time and frequency response.

Table 1

Comparative analysis of linear and nonlinear model.

	Step input		Bump input		Class C road input
	Linear	Non linear	Linear	Non linear	Linear	Non linear
VDV	3.7402	5.8180	1.5784	1.6639	3.6110	5.4373
Awh	1.0369	1.5636	0.7023	0.7350	1.0737	1.5685
Max awh	10.8769	16.4397	2.7950	2.9612	14.8633	18.9152
CF	10.4901	10.5141	3.9024	3.9468	11.5218	12.0593
AR_h	1.1225	1.1487	1.1255	1.1276	1.1133	1.1273
AR_ut	1.1199	1.1457	1.1206	1.1229	1.1113	1.1250
Aw_spr	1.2668	1.6464	2.0467	2.0672	0.5247	0.7470

Fig. 1b also represents frequency response of linear and nonlinear quarter car with a human model. For linear and nonlinear models, the first peak at head acceleration is observed at 3 Hz whereas the second peak at wheel hop occurs at 10 Hz. For upper torso frequency response, the first peak is observed at upper torso acceleration at 3 Hz and the second peak is observed at wheel hop at 10 Hz for both linear and nonlinear models. The frequency response of sprung mass acceleration shows the first peak at sprung mass acceleration at 3 Hz for both linear and nonlinear models whereas the second peak at wheel hop is observed at 10 Hz for the linear system and 14 Hz at the nonlinear system. This is due to quadratic nonlinearity in tyre and cubic nonlinearity in suspension spring.

From frequency response plots, shown in Fig. 1b, it is observed that the magnitude of head acceleration gain, upper torso acceleration gain and sprung mass acceleration gain of nonlinear system is greater than that of linear system. This indicates more acceleration transmission in nonlinear system due to nonlinearities. This is also evident from the time responses of nonlinear system. Refer Fig. 1b. Similar trends are observed for bump response. Refer Fig. 1c for time and frequency responses and Table 1 for results.

Fig. 1c

Bump input – linear and nonlinear quarter model – time and frequency response.

McGee et al. [13] already presented a detailed study on frequency domain using laboratory shaker data. Hence for class C road input, only time domain results are shown in Fig. 1d. It is observed that RMS head acceleration, VDV at head, upper torso acceleration, and sprung mass acceleration are on the higher side for a nonlinear model as compared to the linear model. Also CF, AR_h and AR_ut are greater for the nonlinear system as compared to the linear system. Hence, nonlinearities should be adequately addressed during the ride, control and optimization applications of vehicle models.

Fig. 1d

Class C road input – linear and nonlinear quarter model – time response.

Multi-objective optimization

Researchers have invented several meta-heuristic optimization algorithms to optimize the problems in several fields. These algorithms have implemented on several mathematical problems involving single objective optimization to multi-objective optimization and provided excellent results.

The suspension system has to perform several conflicting objectives such as ride comfort, road holding, and suspension/rattle space requirements. Also, in this study, human model is incorporated to optimize the objective functions considering the human body responses rather than only the seat. Thus, the optimization problem becomes multi-objective in nature (consisting of Head VDV, RMS head acceleration, crest factor, AR_h, AR_ut, suspension space requirement and dynamic tyre force/deflection as objective functions) with conflicts. As compared to a single objective optimization problem, a multi-objective optimization (MOO) problem has to satisfy several objectives simultaneously. Hence, multi-objective optimization using genetic algorithm (GA) and particle swarm optimization (PSO) algorithms are implemented to solve the optimization problem.

In the solution of MOO problems, MOO forms a Pareto optimal front consisting of multiple optimal solutions. Genetic Algorithm (GA) is implemented to optimize in multiple domains as it handles complex optimization problems with discontinuities, non-differentials, noisy functions and functions with multi-modality. GA also supports parallel computations with obtaining Pareto front in a single run. Non-dominated sort GA-II (NSGA-II) is one of the MOEAs using GA strategy. NSGA-II implements non-dominated sort algorithm thus reducing the computational complexities. While sorting the parents and children, elitism is introduced in NSGA-II. In NSGA-II, to preserve the diversity and uniform spread of optimal front, a crowding distance (CD) operator is used. Chromosomes with better fitness are assigned highest ranks, and thus, they determine the domination [22].

Multi-objective PSO – crowding distance (MOPSO-CD) is one of the variants of MOPSO family. It uses PSO algorithm to handle MOO problems. It uses external archive/repository of non-dominated solutions to store the global best solutions, thus maintaining elitism. MOPSO used CD operator to select the global best solution and deletion method of the external archive. Along with CD operator, the mutation operator is used to maintain diversity amongst the solutions of the external archive. The non-dominated solutions stored in the external archive are used to guide the particle search [23].

Hence due to these merits, NSGA-II and MOPSO-CD have been widely implemented to solve MOO problems.

Genetic algorithm

Genetic algorithm, invented by Holland [24], [25], is a meta-heuristic optimization algorithm based on the principle of genetics and natural selection. As random numbers are generated during the operation of genetic algorithms hence GA is stochastic algorithms. These random numbers generated determine the search result [24], [25], [26]. For multi-objective optimization, NGPM code (NSGA-II Program in Matlab) is used [27], [28]. NGPM is the implementation of NSGA-II (Non-dominated Sort Genetic Algorithm) in Matlab.

Firstly, non-dominated sorting is done using NSGA-II by comparing each individual with remaining solutions of a population [22] and thus, all non-dominated solutions and non-dominated fronts are identified and ranked. For rank 1 individuals, fitness value 1 is assigned. For rank 2 individuals, fitness value 2 is assigned and so on [22].

A new parameter, Crowding Distance (CD), is introduced by NSGA-II [22]. CD is the measure of diversity of individuals in the non-dominated population. After completing the sorting, CD is assigned to each individual, front-wise. More the CD more is the diversity in the population. Individuals in the boundary are always selected as they have assigned infinite CD.

From the non-dominated front, parents are selected on the basis of tournament selection and comparing the CD. New offsprings are created using crossover operator and mutation operator. New offsprings and current population (parents) are combined to generate a new population. Selection is carried out for next generation individuals.

The binary tournament selection method is used by NSGA-II to handle constraints. In this method, a solution either feasible or infeasible is decided by comparing with other solution. Here constrained dominate solution between two solutions is identified by using following rule –

A solution i is said to be constrained dominate a solution j if

i is a feasible solution, whereas j is not.

Both solutions, i and j, are infeasible; however, i has overall constraint violation smaller as compared to j.

i and j solutions are feasible, but i dominates j.

Here, the number of generations is used as stoppage criterion. Fig. 2a explains the flow chart of GA algorithm implemented for multi-objective optimization.

Fig. 2a

Flow chart – GA.

MOPSO-CD

PSO algorithm proposed by Eberhart and Kennedy [29] was inspired by the social behaviour of birds’ flock searching randomly for food. While searching for food, instead of knowing the exact location of food, birds know their current location from food, and birds searches the bird which is closest to the food. Thus, PSO is a population-based algorithm where each bird is known as a particle. The particle flies through the solution space (or search space) to search the optimal solution. Each particle flies through the search space with a velocity which is determined according to the flying experience of bird’s own flying and its flock. All the particles have objective function value as per objective function. Each particle updates its position in search space based on its current location and previous best location (also called as pbest) and best location of the whole population (flock) (also called as gbest), current velocity. In PSO algorithm, particles’ initial positions and initial velocities are randomly initialized.

The new velocity and new position of every particle can be determined using following equations:where w is inertia weight, c1 and c2 are acceleration coefficients and r1 and r2 are the uniform random numbers in the range [0, 1]. The optimization is an iterative process and thus converges to the global optimal solution.

In single objective optimization, it is easy to calculate pbest and gbest values, as compared to multi-objective optimization. Due to multi-objective nature of the problem and conflict between the objectives, it is quite difficult to calculate the pbest and gbest values. It is impossible for all objective functions to reach maximum values or minimum value at the same time. Thus, in multi-objective PSO (MOPSO) uses the Pareto ranking scheme to take care the multi-objective problem. In MOPSO non-dominated solutions are stored in the archive where the historical records of best solutions are obtained [30].

MOPSO-CD [23] algorithm includes crowding distance (CD), similar to NSGA-II, computation mechanism in the PSO algorithm to solve multi-objective problems. The CD mechanism is used in selection of gbest and in deletion of an external archive of non-dominated solutions. To maintain diversity in the non-dominated archive, the mutation operator is used along with CD mechanism. The global best, gbest, is selected from those having highest CD values. The non-dominated solutions are moved in an external repository; A. External repository A is having solutions with least crowded objective space. Fig. 2b explains the flow chart of PSO optimization algorithm.

Fig. 2b

Flow chart – PSO.

Constraint handling in MOPSO-CD

In multi-objective constrained optimization, the key issue is the constrained handling technique. Here, penalty function method is used for constrained handling as it is simple yet has good convergence [31].

Let us consider the multi-objective problem as follows:where x = [x1, x2, … x] is a decision variable in a decision space X.

In this paper, only inequality constraints are considered for optimization. Here, for constraint handling in the multi-objective optimization problem, the degree of violation of constraint is used in penalty function method. Each particle is checked for the constraint violation. Upon violation of constraint by particle x, the degree of constraint violation is defined aswhere δ = extent of constraint violation and δ = 0 if a constraint is not violated when ith constraint is satisfied. C = Penalty.

Objective functions

One of the key factors in optimization problems is to choose proper objective functions. As stated earlier, body parts respond to a much greater extent as compared to seat response. Hence, suspension system is aimed to optimize considering responses such as RMS head acceleration, Vibration dose value (VDV) at head, the amplitude ratio of head RMS acceleration to seat RMS acceleration, the amplitude ratio of upper torso RMS acceleration to seat RMS acceleration, and crest factor as these are the most important factors which affect the driver’s health. Along with these objective functions, suspension space deflection and dynamic tyre force are also included as objective functions to form a multi-objective optimization problem.

RMS head acceleration: As per ISO 2631-1 [32], RMS head acceleration is given by

A major portion of the vibration experienced by the occupants of an automobile enters the body through the seat [33], [34]. Whole-body vibrations, mostly affect the human body, are characterized by vertical vibrations transmitted to the buttocks and back of the occupant along the vertebral axis via the base. The health risk increases as the exposure time to vibrations increases. Hence is it necessary to measure the whole body vibrations. As per ISO 2631-1 [32], VDV is one measure for whole body vibrations. VDV is also called as fourth power vibration dose. VDV is the method of assessing the cumulative effect (dose) of the vibration.

Vibration Dose Value (VDV): VDV is the fourth power of acceleration time histories. It is expressed as follows:

Crest Factor (CF): It is defined as the ratio of maximum head acceleration to the RMS head acceleration [32].

Amplitude ratio of head RMS acceleration to seat RMS acceleration (AR_h): It defined as the ratio of head RMS acceleration to seat RMS acceleration.

Amplitude ratio of upper torso RMS acceleration to seat RMS acceleration (AR_ut): It is defined as the ratio of upper torso RMS acceleration to seat RMS acceleration.

Suspension travel: Suspension travel is characterized by the relative travel between the sprung mass and unsprung mass. Due to random input, RMS suspension space travel is taken as one of the objective functions.

Dynamic tyre force: Dynamic tyre force is related to tyre deflection. Due to random nature of input RMS of tyre deflection is next objective function.

According to Baumal and et al. [8], at least, 125 mm of suspension travel is required and maximum seat acceleration should not increase 4.5 m/s2 so as to avoid hitting the suspension stop. To minimize dynamic tyre forces, maximum tyre deflection should not increase 0.058 m. These parameters are included as constraints in the optimization problem.

The formulation of optimization problem is as follows:

fobj1 = Minimize (VDVh)

fobj2 = Minimize (Awh)

fobj3 = Minimize (CF)

fobj4 = Minimize (AR_h)

fobj5 = Minimize (AR_ut)

fobj6 = Minimize (RMS suspension travel)

fobj7 = Minimize (RMS tyre deflection)

Subject to constraints:

Search space:

During optimization, the design parameters are suspension spring stiffness; suspension damping and seat cushion parameters. The search space is:

, [10] Hence, , , [35]

Results and discussion

Test problems

To test the performance of NSGA-II and MOPSO-CD with penalty function, four benchmark problems are selected and simulated. Each benchmark function has two objective functions, two constraint functions and two design variables.

Constr problem:

Subject to

TNK problem

Subject to

BHN problem

Subject to

SRN problem

Subject to

Both algorithms are run using following parameters:

NSGA-II: Population size = 100, number of generations = 100, thus having 10,000 evaluation count.

MOSPO-CD: Number of particles 100 with Maximum number of particles in repository A = 100 [30], and Maximum evaluation = 10,000.

Both algorithms are simulated in Matlab/Simulink environment. Pareto front and results are shown graphically in Fig. 3a, Fig. 3b respectively. Pareto fronts of both algorithms are in close agreement with each other. Results are also tabulated in Table 2. The table shows minimum value, the maximum value and the average value of each benchmark problem/function from the set of solutions of 100. Also, Table 2 shows elapsed time for each problem.

Fig. 3a

Pareto front for benchmark problems.

Fig. 3b

Comparison of benchmark functions.

Table 2

Benchmark function values.

Benchmark problem	Optimization algorithm	Objective function	Minimum value	Maximum value	Average value	Elapsed time (s)
Constr problem	MOPSO-CD	f1(x)	0.391933	1	0.834189	11.11693
		f2(x)	1	8.929097	2.161154
	NGPM (NSGA-II)	f1(x)	0.391348	1	0.616499	36.02496
		f2(x)	1	8.91946	3.798027
TNK problem	MOPSO-CD	f1(x)	0.048698	1.036299	0.647836	10.97732
		f2(x)	0.048177	1.046548	0.623831
	NGPM (NSGA-II)	f1(x)	0.078051	1.04836	0.666922	40.06497
		f2(x)	0.056232	1.05603	0.632637
BHN problem	MOPSO-CD	f1(x)	0	136	70.43662	9.913939
		f2(x)	4	50	18.52125
	NGPM (NSGA-II)	f1(x)	0	136	50.81124	28.55074
		f2(x)	4	50	19.64545
SRN problem	MOPSO-CD	f1(x)	5.729107	195.7479	92.60403	7.154332
		f2(x)	−217.471	3.609806	−108.957
	NGPM (NSGA-II)	f1(x)	5.99457	198.688	94.82053	19.75622
		f2(x)	−216.96	7.87181	−110.502

In each benchmark problem, both objective functions are of minimization nature. From Table 2 and Fig. 3b, for Constr problem it is observed that minimum, maximum, and average function values are less for NSGA-II as compared to MOPSO-CD. However, the average function value of f2(x) function is much less for MOPSO-CD algorithm as compared to NSGA-II. In case of TNK problem, it observed that minimum, maximum and average function values are less for MOPSO-CD as compared to NSGA-II. For BHN problem, minimum and maximum function values are same for both algorithms. However, the average value of f1(x) function is greater for MOPSO-CD algorithm as compared to NSGA-II, whereas the average value of f2(x) function is smaller for MOPSO-CD algorithm as compared to NSGA-II. For SRN problem, it observed that minimum, maximum and average function values are less for MOPSO-CD as compared to NSGA-II. From Table 2, it is clear that MOSPO-CD takes much lesser computation time as compared to NSGA-II algorithm for benchmark problems.

Validation on a nonlinear quarter car

Multi-objective optimization of nonlinear quarter car model is simulated in Matlab/Simulink® environment using NSGA-II and MOPSO-CD algorithms.

The stationary road roughness is effectively described by power spectral density (PSD) [36]. When a car moves at a constant velocity u, the road roughness can be viewed as a stationary process in space domain, and the PSD of the road disturbance input can be expressed bywhere is the road PSD, n is the spatial frequency, no is the reference spatial frequency and no = 0.1(cycles/m). is the coefficient of road roughness, and ω is waviness and is generally ω = 2.

In frequency domain, PSD can be expressed aswhere ωo is the lowest cut-off angular frequency and ωo = 2πfo = 2πvno.

The differential equation of road roughness can be expressed as [36]

Eq. (24) is modelled in Simulink environment to model the road surface. Input road condition is modelled as class C road (average road) with degree of road roughness 512 × 10−6 m2/(cycle/m) [36], [37]. The vehicle is travelling with velocity of 80 kmph. Fig. 4 represents the class C road surface.

Fig. 4

Road surface (Class C, velocity 80 kmph).

For GA, the range of design variable ks is (35,250 − 11,750=) 23,500. Thus, design variable ks needs to be divided into 23,500 equal range of size. Hence 214 = 16,384 < 23,500 < 215 = 32,768 i.e. 15 bits are required to store value of design variable ks in the chromosome. Similarly for other design variables, the bit required is tabulated in Table 3.

Table 3

Design variables range.

Design variable	Range	Size
ks	35,250 − 11,750 = 23,500	214 = 16,384 < 23,500 < 215 = 32,768	15 bits
cs	1050 − 350 = 700	29 = 512 < 700 < 210 = 1024	10 bits
kc	20,000 − 2500 = 17,500	214 = 16,384 < 17,500 < 215 = 32,768	15 bits
cc	1649.03 − 131.59 = 1517.44	210 = 1024 < 1517.44 < 211 = 2048	11 bits

Hence, the total length of chromosome or gene is 15 + 10 + 15 + 11 = 51 bits are required where first 15 bits are required for ks, next 11 bits are required for cs, next 15 bits are required for kc and next 11 bits are required for design variable cc.

A lot of work is done on determining the optimum population size for GA based optimization. Rosenthal and Borschbach [38] in their study concluded that the best performance for NSGA-II is achieved with a population size ranging from 70 to 100. Hernandez-Diaz et al. [39] experimentally found that NSGA-II needs minimum population size of 52 to reach a good performance. Reeves and Rowe [40] derived a formula based on the principle that, in the search space every point should be reachable from the initial population. In GA, the initial population is randomly generated. Therefore, the probability that at least one allele is present at each locus is given by equation:where N = population size and ls = string length

Using exponential function approximation, the population size is given by equation:

According to Eq. (26), population size 17 is sufficient for the string of length 50 to exceed the probability of 99.9%.

According to Alander [41], population size is given by relation:

Thus selecting population size of 100 for optimization and the stoppage criterion is the number of generations and is equal to 100 generations.

In MOPSO-CD, algorithm 100 particles are used with a repository of size 100 to store non-dominated solutions and a maximum number of evaluations 100,000 to stop the algorithm [30].

Multi-objective optimization of nonlinear quarter car suspension driver system is established using GA and PSO to search optimum parameters of the suspension system (i.e. spring stiffness and damping coefficient) and cushion (i.e. cushion material stiffness and damping). The aim was to minimize VDV, RMS sprung mass acceleration, crest factor, AR_h, AR_ut, suspension space deflection and tyre deflection.

Fig. 5a shows the trade-off fronts obtained for NSGA-II and MOPSO-CD optimization algorithm with four design variables i.e. ks, cs, kc and cc. Next step is to invoke the decision procedure to select a set of design variables from given sets of design variables. In this case, 100 sets of design variables are available for NSGA-II as well as MOPSO-CD algorithm approach, as population size is 100 for NSGA-II algorithm and repository size is 100 for the MOPSO-CD algorithm. This gives the designer more flexibility to choose a particular set of design variables as a solution for further analysis.

Fig. 5a

Trade-off Front for NSGA-II and MOPSO-CD.

Hence, from the 100 sets of design objective functions, minimum value of each objective function is tabulated in Table 4 for both algorithms. Table 4 also represents the objective function values for the classical parameter.

Table 4

Objective function values.

Objective functions	NSGA-II	MOPSO-CD	Classical
VDVh	1.3178	1.3171	3.1028
RMS Awh	0.5491	0.5481	1.3504
CF	3.8099	3.7961	3.8418
STHT	0.8804	0.8875	1.1058
AR	0.8721	0.8712	1.1033
RMS suspension space	0.0189	0.0196	0.02248
RMS tyre deflection	0.0108	0.0114	0.0100
Constraints
Maximum ah	2.8195	2.8224	5.5455
Suspension space deflection	0.0593	0.0633	0.06424
Tyre deflection	0.0364	0.0370	0.0344
Optimization time (s)	9345.2454	8858.4179	–

To represents the results more clearly, design variables corresponding to the minimum value of RMS Awh from Table 4 are selected for further analysis. The design variables correspond to the minimum value of RMS Awh are tabulated in Table 5. The comparison of objective functions corresponding to minimum RMS Awh for NSGA-II and MOSPO-CD with classical values is shown in Fig. 5b. From Fig. 5b, it is observed that RMS head acceleration and VDV at the head are very less for optimized variables (i.e. NSGA-II and MOPSO-CD algorithm) as compared to classical values. But there is a slight reduction in crest factor of optimized variables as compared to classical variables. Although the crest factors are nearly the same for optimized variables and the classical one, still RMS head acceleration is very less for optimized variables thus providing good ride over classical one. The AR_h and AR_ut are reduced by 20% each for NSGA-II algorithm. Also, for MOPSO-CD algorithm, the reduction in AR_h and AR_ut is 19% and 20% respectively. Hence, optimized parameters provide better ride comfort and health by reducing RMS head acceleration and VDV at head, AR at head and upper torso.

Table 5

Design variables.

	ks	cs	kc	cc
NSGA-II	11,750	546.4720	2510	1509.73
MOPSO-CD	11,750	542.1716	2500.2333	1500

Fig. 5b

Objective Function Values corresponding to minimum Aw.

The responses are simulated for Class C road at 80 kmph vehicle speed. Fig. 5c shows simulation responses of head acceleration, upper torso acceleration, suspension space deflection and tyre deflection for optimized design variables corresponding to the minimum value of Awh along with classical values. In Fig. 5c, it is observed that head acceleration, upper torso acceleration of driver is considerably improved by optimized suspension–seat system.

Fig. 5c

Time response of quarter car seat–suspension–driver model.

Fig. 5d shows the frequency response of un-optimized parameters and optimized parameters using NSGA-II and MOPSO-CD algorithms. It is observed that head acceleration gain, upper torso acceleration gain, suspension deflection gain and tyre deflection gains are on the lower side for optimized NSGA-II and MOPSO-CD parameters compared to un-optimized parameters.

Fig. 5d

Frequency response of quarter car seat–suspension–driver model.

To check the performance on varying speed, optimized variables (shown in Table 5) are further simulated on class C road surface for speed ranging from 50 to 120 kmph. It is observed that RMS head acceleration increases with increase in speed. Maximum RMS head acceleration is 0.5536 m/s2 and 0.5623 m/s2 at 120 kmph for NSGA-II and MOPSO-CD, respectively. A similar trend is observed for VDV at the head, and it increases with increase in speed. Maximum VDV at the head is 1.3653 m/s1.75 and 1.3631 m/s1.75 at 120 kmph for NSGA-II and MOPSO-CD, respectively. Refer Fig. 6a.

Fig. 6a

Performance of optimized variables at various speeds.

Amplitude ratio AR_h and Amplitude ratio AR_ut also increase with increase in speed from 50 to 90 kmph. But, from 90 kmph slight decrease in both ratios is observed. Maximum AR_h is 0.8851 and 0.8908 at 90 kmph for NSGA-II and MOPSO-CD, respectively. At 120 kmph AR_h is 0.8819 and 0.8870 for NSGA-II and MOPSO-CD, respectively. Similarly, maximum AR_ut is 0.8770 and 0.8828 at 90 kmph for NSGA-II and MOPSO-CD, respectively. At 120 kmph AR_ut is 0.8716 and 0.8768 for NSGA-II and MOPSO-CD, respectively.

The factors associated with the assessment of health risk are the magnitude of vibration and duration of vibration. In this study, the magnitude of vibration is determined according to frequency weighted RMS head acceleration as per ISO 2631-1 [32] standard. Annexure B of ISO 2631-1 provides a graphical representation of Health Guidance Caution Zone (HGCZ). The graph has two criterions for assessment of health prediction – first is based on exposure duration and vibration magnitude as per RMS value of Aw and second based on VDV. Here, the first criterion is used to assess prediction of health risk and shown in Fig 6b. This criterion is divided into three different exposure levels – no health risk (NHR), potential health risk (PHR) and likely health risk (LHR), refer Fig. 6b.

Fig. 6b

Annexure B – ISO 2631-1 (Health Guidance Caution Zone, HGCZ) [32].

Table 6 represents RMS Awh for 10-min duration. For 10 min period, it is observed that the RMS Awh values lie well within the NHR zone for both optimization methods. Fig. 6b shows the RMS head acceleration for 120 kmph vehicle speed, represented by a blue line. From Fig. 6b, for 0.5635 m/s2 RMS head acceleration, vibration exposure limit is up to 6 h for NHR and 6–17 h for PHE and 17–24 h is for LHR. Table 6 shows the time duration of both optimization algorithms according to exposure levels as per ISO 2631-1 criterion. Thus, optimized variables improve ride and health criterions.

Table 6

Time limits as per Health Guidance Caution Zone (HGCZ) as per ISO – 2631-1.

Speed (kmph)	NSGA-II (Class C road)				MOPSO-CD (Class C road)
	Awh (m/s2)	Vibration exposure limit (h)			Awh (m/s2)	Vibration exposure limit (h)
		NHR	PHR	LHR		NHR	PHR	LHR
50	0.5138	7	7–18	18–24	0.5133	7	7–18	18–24
60	0.5300	6	6–17	17–24	0.5293	6	6–17	17–24
70	0.5411	6	6–17	17–24	0.5403	6	6–17	17–24
80	0.5491	6	6–17	17–24	0.5480	6	6–17	17–24
90	0.5547	6	6–17	17–24	0.5537	6	6–17	17–24
100	0.5587	6	6–17	17–24	0.5576	6	6–17	17–24
110	0.5616	6	6–17	17–24	0.5604	6	6–17	17–24
120	0.5635	6	6–17	17–24	0.5623	6	6–17	17–24

Conclusions

This paper presents multi-objective optimization of Nonlinear Quarter Car Seat Suspension system with a driver model. An eight DoF model consisting of 2 DoF nonlinear quarter car model having quadratic nonlinearities in tyre and cubic nonlinearity in suspension spring, 4 DoF driver biomechanical model, seat frame and the seat cushion are developed and optimized.

In optimum design problem comfort and health criterion consisting of VDV at the head, RMS head acceleration, crest factor, AR_h, AR_ut, along with stability criterions consisting of suspension space and tyre deflection are used as objective functions. ISO 2631-1 methodology is adopted to assess the objective functions such as RMS head acceleration, VDV at head and crest factor.

Optimization of the car suspension–seat–driver system is successfully implemented using NSGA-II and MOPSO-CD with penalty function algorithms. The MOPSO-CD algorithm takes less computation time as compared to NSGA-II for optimization.

Numerical simulations are presented for optimum design variables of a quarter car suspension–seat–driver system obtained by implementing NSGA-II and MOPSO-CD algorithms. Results of a quarter car travelling over a Class C road (average road) at speed 80 kmph are presented to show its performance. For class C road, RMS head acceleration and VDV at head increase with increase in speed. Simulation result shows that optimum design variables improve ride comfort and health criterions over classical design variables.

Conflict of interest

The authors have declared no conflict of interest.

Compliance with Ethics Requirements

This article does not contain any studies with human or animal subjects.

mh = 5.31;	ch = 400;	kh = 310,000;
mut = 28.49;	cut = 4750;	kh = 183,000;
mlt = 8.62;	clt = 4585;	klt = 162,800;
mt = 12.78;	ct = 2064;	kt = 90,000;
mc = 1;	cc = 200;	kc = 18,000;
mf = 15;	cf = 830;	kf = 31,000;
ms = 290;	cs = 700;	ks = 23,500;	ksnl = 100ks[15];
mus = 40;	kt = 190,000;	ktnl = 1.5 kt[14].