# Pareto Optimization of a Half Car Passive Suspension Model Using a Novel Multiobjective Heat Transfer Search Algorithm

Authors: Vimal Savsani; Vivek Patel; Bhargav Gadhvi; Mohamed Tawhid
Journal: Modelling and Simulation in Engineering
 (2017)

---

## Abstract

Most of the modern multiobjective optimization algorithms are based on the search technique of genetic algorithms; however the search techniques of other recently developed metaheuristics are emerging topics among researchers. This paper proposes a novel multiobjective optimization algorithm named multiobjective heat transfer search (MOHTS) algorithm, which is based on the search technique of heat transfer search (HTS) algorithm. MOHTS employs the elitist nondominated sorting and crowding distance approach of an elitist based nondominated sorting genetic algorithm-II (NSGA-II) for obtaining different nondomination levels and to preserve the diversity among the optimal set of solutions, respectively. The capability in yielding a Pareto front as close as possible to the true Pareto front of MOHTS has been tested on the multiobjective optimization problem of the vehicle suspension design, which has a set of five second-order linear ordinary differential equations. Half car passive ride model with two different sets of five objectives is employed for optimizing the suspension parameters using MOHTS and NSGA-II. The optimization studies demonstrate that MOHTS achieves the better nondominated Pareto front with the widespread (diveresed) set of optimal solutions as compared to NSGA-II, and further the comparison of the extreme points of the obtained Pareto front reveals the dominance of MOHTS over NSGA-II, multiobjective uniform diversity genetic algorithm (MUGA), and combined PSO-GA based MOEA.


---

## Body

## 1. Introduction

In recent time, the multiobjective evolutionary algorithms (MOEAs) have gained enormous attention in solving the engineering optimization problems with more than one objective. The multi/many objective optimization problems (MOOPs) differ from their single objective counterparts (SOOPs) in terms of both their problem definitions/statements and methods to solve such problems. MOOPs have the conflicting objectives to be optimized simultaneously and solving these yields a set of optimal or trade-off or Pareto solutions (Pareto front points), whereas SOOPs have a single objective at a given time and the solution is usually a single optimal point. The classical methods such as calculus based methods, gradient based methods, and elimination and interpolation methods and the like for solving SOOPs and methods, for instance, weighted sum method andε constraint method for solving MOOPs are limited to the problems with simple objective functions and constraints. This is due to their nature to get stuck to a suboptimal solution, their convergence dependency on initial guess, and their unsuitability in solving a large variety of optimization problems [1]. Alternatively evolutionary algorithms (EAs), or metaheuristics, have had a remarkable success in finding the global optimum of complex problems nearly in all the disciplines of the knowledge. This success of EAs and their nature of using a population of solutions had led the researchers to employing the search techniques of EAs to optimize the MOOPs. Such algorithms to solve MOOPs are referred to as multiobjective evolutionary algorithms (MOEAs) when they use EAs as their basic search techniques and are referred to as simply multiobjective optimization algorithms (MOOAs) when they use any metaheuristics in general. Primarily all the methods to solve MOOPs have two goals to attain. The first is to find the Pareto front solutions as close as possible to the optimal Pareto front and the second is to maintain diversity among the optimal set of solutions [1].One of the oldest attempts to have employed the EA to form a population based multiobjective evolutionary algorithm was by Schaffer [2] and the method is known as vector evaluated genetic algorithm (VEGA). VEGA does the selection for each objective separately [3], but its incapability in finding all the Pareto front solutions had limited its applications to very few real world optimization problems. These demerits of VEGA were overcome by algorithms which were inspired by the work of Goldberg and Holland [4], such as multiobjective genetic algorithm (MOGA) by Knowles and Corne [5], nondominated sorting genetic algorithm (NSGA) by Srinivas and Deb [6], and niched Pareto genetic algorithm (NPGA) by Horn et al. [7]. These new algorithms introduced the concept of Pareto optimality into the selection process which awarded them a great success in application to various disciplines of engineering as thoroughly described in [8]. Consequently the researchers had started experimenting various ways in assigning fitness values to the populations and in maintaining the diversity of optimal points which led to the development of new state-of-the-art MOEAs, for example, strength Pareto evolutionary algorithm (SPEA2 and SPEA) by Zhou et al. [9] and by Zitzler and Thiele [3], Pareto archived evolution strategy (PAES) by Knowles and Corne [5], Pareto envelope-based selection algorithm (PESA and PESA-II) by Corne et al. [10] and Horn et al. [7], an elitist based nondominated sorting genetic algorithm-II (NSGA-II) by Deb et al. [11], and a multiobjective evolutionary algorithm based on decomposition (MOEA/D) by Zhang and Li [12] and its improved versions described in [9]. In addition, swarm intelligence based search algorithms for multiobjective optimization [13–15] have also been developed and applied to a variety of MOOPs. Several more recently developed multiobjective algorithms [16–19] based on the search techniques other than the ones used in genetic algorithm [4] or EAs and particle swarm optimization [20] have also shown a great potential in obtaining the Pareto front closest to the true Pareto front.In this study a nonlinear MOOP of optimizing the vehicle suspension system design having conflicting objectives is optimized. The goal is to optimize the design of passive suspension system of the passenger road vehicle by using a 2-dimensional pitch-heave ride model, also known as half car model. This half car ride model had been used by many researchers as an MOOP to test the performance of their proposed MOOAs, for instance, multiobjective uniform diversity genetic algorithm (MUGA) by Jamali et al. [21] and hybrid PSO and GA algorithm for MOOPs by Mahmoodabadi et al. [22]. The following paragraph illustrates the MOOP of half car ride model.Vehicle suspension system has essentially three objectives, good ride comfort, good road holding or handling, and limited suspension deflection or suspension rattle space. These objectives are conflicting with one another and the designers developing such products are always seeking an optimized passive suspension system fulfilling largely all three conflicting objectives simultaneously. In the present study the multiobjective optimization problem is solved by employing a lumped mass-spring-damper vibration model for the half car vehicle suspension system and then applying MOOAs to get the desired set of optimal design solutions, that is, trade-off or Pareto optimum design points. The employed 2D half car model has five degrees of freedom with one lumped sprung and two lumped unsprung masses (front and rear tires) and thus has a set of five second-order, linear, ordinary differential equations of motion.In the present paper a multiobjective heat transfer search (MOHTS) algorithm for solving MOOPs is proposed and has been applied to the MOOP of optimizing the suspension parameters of the half car model with driver’s seat having total five degrees of freedom as used in [21–23] to evaluate the potential of the proposed MOHTS algorithm. Moreover, optimization results of MOHTS are compared with the ones obtained from the popular MOEA named NSGA-II and with the other algorithms which were developed and applied to solve the multiobjective half car optimization problem. The proposed algorithm (MOHTS) employs a very recently developed metaheuristic algorithm named heat transfer search algorithm [25] as a base search technic. Additionally, for sorting of the population an elitist nondominated sorting approach and for maintaining the diversity among multiple optimal solutions of each nondominated front the crowding distance approach of NSGA-II are adopted.The second section briefly introduces the basic concepts and definitions pertaining to multiobjective optimization. The third section describes the overview of the search technic of heat transfer search (HTS) algorithm, elitist based nondominated sorting and crowding distance approach of NSGA-II, and then explains the architecture and working of the proposed MOHTS algorithm. The fourth section of the paper illustrates two different half car multiobjective optimization problems which are solved numerically, elucidating the potential of the proposed MOHTS algorithm. Firstly the half car MOOP with five conflicting objective functions as used in [21–23] is employed in numerical study 1 to test and compare the performance of MOHTS with NSGA-II, MUGA [21], and with the combined PSO and GA based MOEA [22]; secondly another half car MOOP with a different set of five more realistic conflicting objective functions from the ones used in numerical study 1 is formulated and tested on the platforms of MOHTS and NSGA-II.


## 2. Terminologies in Multiobjective Optimization

A multiobjective optimization problem can be stated as follows [1, 26].Find the design vectorX=x1x2⋯xnT, which minimizes(1)FX=f1Xf2X⋯fmXTsubject top inequality and q equality constraints, respectively, and bounds on design vector members (2)gjX≤0,forj=1,2,…,p(3)hkX=0,fork=1,2,…,q(4)xiL≤xi≤xiU,fori=1,2,…,n.Here,X is the vector of design variables known as design vector with its bounds as shown in (4); X forms an n-dimensional Cartesian space called design space D created by representing each variable with one coordinate axis. F is the vector of objectives or decisions known as objective vector which is to be minimized meaning all the objectives are to be minimized simultaneously. Unlike the single objective optimization problem (SOOP), MOOP creates a multidimensional objective space in the present case it is an m-dimensional objective space with Cartesian each axis as one objective function. Equations (2) and (3) show the inequality and equality constraints, respectively, which must be satisfied by the design vector in order to create the feasible solutions. The set of all feasible solutions is known as the feasible region S. In general, there exists no solution vector X which can minimize all the m objective functions simultaneously; hence a concept called Pareto optimum solution which is used in solving MOOPs is now introduced.


### 2.1. Pareto Dominance and Optimality

For the multiobjective minimization problem as stated above, the set of optimum solutions (Pareto optimal design vectors) includes those vectors which cannot minimize any member of the corresponding objective vector without deteriorating another member [3]. To illustrate the Pareto optimality mathematically, consider two design vectors u, v∈S; then vector u is said to be dominant over vector v (denoted by u≺v) iff ∀l∈1,2,…,m, Fl(u)≤Fl(v)∧∃l∈1,2,…,mFl(u)<Fl(v). In a given set of solutions, the design vectors that are nondominated by any other design vector of the same set are known as nondominated regarding that set. Furthermore, the design vectors which are nondominated over the entire feasible region S are known as Pareto optimal solutions and these solutions collectively form a set called Pareto optimal set (PS) or Pareto optimal front (PF).


### 2.2. Solving MOOPs by Extending EAs and Metaheuristics

Evolutionary or population based algorithms for solving SOOPs yield a set of solutions known as populations instead of a single solution after one iteration. This nature of metaheuristics in obtaining multiple solutions in each iteration makes it suitable to solve MOOPs. As discussed in the first section, MOOPs have primarily two goals: the first is to obtain the Pareto front as close as possible to the optimal Pareto front and the second is to preserve the diversity between the Pareto optimal solutions. The first aim can be achieved by choosing an appropriate fitness assignment scheme that prefers nondominated solutions, and the second can be accomplished by using an appropriate strategy that preserves the diversity among the solutions of each nondominated front. For instance, NSGA [6] used nondominated sorting as a fitness assignment procedure in which each individual from the population is compared with the others to find its nondominancy and thus to obtain the first Pareto front this will be followed by sorting all the nondominated individuals from the population and this will be repeated until the entire population is sorted in different fronts. In addition, for maintaining the diversity among the solutions in a front, a front wise sharing function method was used in NSGA which calculates the Euclidean distance for each solution in the Pareto front from another solution in the same front. The performance of NSGA was further improved in terms of converging the solution more close to the Pareto front by employing an elitist based nondominated sorting, in maintaining a widespread of solutions by applying crowding distance approach instead of sharing function method, and in reducing the computation time. This improved NSGA was named NSGA-II [11]. The proposed algorithm (MOHTS) adopts the fitness assignment and diversity preserving mechanisms of NSGA-II over the search method of HTS for solving MOOPs.


## 2.1. Pareto Dominance and Optimality

For the multiobjective minimization problem as stated above, the set of optimum solutions (Pareto optimal design vectors) includes those vectors which cannot minimize any member of the corresponding objective vector without deteriorating another member [3]. To illustrate the Pareto optimality mathematically, consider two design vectors u, v∈S; then vector u is said to be dominant over vector v (denoted by u≺v) iff ∀l∈1,2,…,m, Fl(u)≤Fl(v)∧∃l∈1,2,…,mFl(u)<Fl(v). In a given set of solutions, the design vectors that are nondominated by any other design vector of the same set are known as nondominated regarding that set. Furthermore, the design vectors which are nondominated over the entire feasible region S are known as Pareto optimal solutions and these solutions collectively form a set called Pareto optimal set (PS) or Pareto optimal front (PF).


## 2.2. Solving MOOPs by Extending EAs and Metaheuristics

Evolutionary or population based algorithms for solving SOOPs yield a set of solutions known as populations instead of a single solution after one iteration. This nature of metaheuristics in obtaining multiple solutions in each iteration makes it suitable to solve MOOPs. As discussed in the first section, MOOPs have primarily two goals: the first is to obtain the Pareto front as close as possible to the optimal Pareto front and the second is to preserve the diversity between the Pareto optimal solutions. The first aim can be achieved by choosing an appropriate fitness assignment scheme that prefers nondominated solutions, and the second can be accomplished by using an appropriate strategy that preserves the diversity among the solutions of each nondominated front. For instance, NSGA [6] used nondominated sorting as a fitness assignment procedure in which each individual from the population is compared with the others to find its nondominancy and thus to obtain the first Pareto front this will be followed by sorting all the nondominated individuals from the population and this will be repeated until the entire population is sorted in different fronts. In addition, for maintaining the diversity among the solutions in a front, a front wise sharing function method was used in NSGA which calculates the Euclidean distance for each solution in the Pareto front from another solution in the same front. The performance of NSGA was further improved in terms of converging the solution more close to the Pareto front by employing an elitist based nondominated sorting, in maintaining a widespread of solutions by applying crowding distance approach instead of sharing function method, and in reducing the computation time. This improved NSGA was named NSGA-II [11]. The proposed algorithm (MOHTS) adopts the fitness assignment and diversity preserving mechanisms of NSGA-II over the search method of HTS for solving MOOPs.


## 3. Multiobjective Heat Transfer Search Algorithm

This section explains the working of the proposed MOHTS algorithm. The different elements of the MOHTS algorithm are basic search technic which is HTS and nondominated sorting method and diversity preserving crowding distance approach of NSGA-II. These elements which constitute the proposed algorithm are briefly explained first followed by the MOHTS procedure section described by the pseudocode for better understanding and software implementation for the readers and users of this algorithm.


### 3.1. Basic Heat Transfer Search Algorithm

A very recently developed metaheuristic algorithm called heat transfer search (HTS) algorithm [25] which is based on the thermodynamic mechanism of the heat transfer has been employed as a basic search strategy in the proposed MOHTS algorithm. The search agents in the algorithm are the molecules which interact with one another and with the surrounding for gaining the thermal equilibrium. Unlike other metaheuristic algorithms, the HTS has been developed by putting more efforts in setting a well-balanced trade-off between the intensification (exploitation) and diversification (exploration). The algorithm is divided into three phases as conduction, convection, and radiation; and in the course of entire search process each phase is performed with an equal probability. All the three phases intensify the search space in initial iterations and later diversify during remaining iterations. The following sections briefly explain all three phases.


#### 3.1.1. Conduction Phase

As stated above in HTS algorithm the population of solutions is represented as number of molecules, design variables are represented as molecules’ temperature level, and fitness values are represented as energy level of molecules. In this phase energy from the molecules with high energy is transferred to the ones with the lower energy. The population is updated during this phase and can be described as follows:If number of iterationsn ≤ maximum number of iterations/CDF (initial populations) andiffXm<fXp then(5)Xp,lnew=1-R2Xm,lold.ElsefXp<fXm; then(6)Xm,lnew=1-R2Xp,lold.Else, for the remaining populations (number of iterationsn> maximum number of iterations/CDF),iffXm<fXp then(7)Xp,lnew=1-rlXm,lold.ElsefXp<fXm; then(8)Xm,lnew=1-rlXp,lold.Here,m and p indicate randomly picked solutions and vary from 1 to u (population size); and l indicates the design variable index ranging from 1 to v (total number of design variables), which is also selected randomly, and CDF is the conduction factor. Moreover, R and rl are probability value for conduction phase (this value is between 0 and 0.3333) and random number between 0 and 1, respectively. For the first half of the total generations (CDF should be taken as 2), a function evaluation comparison is made between two randomly selected solutions; the inferior of the two is replaced by (5) and (6). For the latter half the inferior solution is replaced by (7) and (8). Thus the switching from (5) and (6) to (7) and (8) obtains both the intensification and diversification for maintaining the algorithm’s potential to seek an optimum solution.


#### 3.1.2. Convection Phase

In this phase, a thermal equilibrium is attained as system’s mean temperature (Xmean) interacts with the surrounding temperature (Xs). The latter temperature is assumed to be the best solution and the solution is updated as follows:If current iteration numbern≤ maximum number of iterations/CVF (initial populations) then(9)Xm,lnew=Xm,lold+RXs-Xmean·absR-rl.Else, for the remaining population,(10)Xm,lnew=Xm,lold+RXs-Xmean·round1+rl.It should be noted that, along with each solutionXm, each design variable is also updated during this phase. Here R is taken in the range from 0.6666 to 1, and CVF (convection factor) is taken as 10. That is to say, for the first one-tenth of the iterations, the algorithm explores the feasible design space and for the remaining iterations it exploits the same, but in actual practice this condition does not take place as exploration (diversification) and exploitation (intensification) are taken care of by the last terms of (9) and (10).


#### 3.1.3. Radiation Phase

In this phase, system (solution) communicates either with the surroundings which is the best solution or within the system which is some other solution to attain the thermal equilibrium. The solutions are improved as described below.If number of iterationsn≤ maximum number of iterations/RF (initial populations) andiffXm<fXp then(11)Xm,lnew=1+RXm,lold-R·Xp,lold.ElsefXp<fXm; then(12)Xm,lnew=1-RXm,lold+R·Xp,lold.Else, for the remaining populations,iffXm<fXp then(13)Xm,lnew=1+rlXm,lold-rl·Xp,lold.ElsefXp<fXm; then(14)Xm,lnew=1-rlXm,lold+rl·Xp,lold.Herem and p are randomly selected solutions for which their fitness values are compared with each other and that determines how each solution is updated using (11), (12), (13), and (14). The range of R in this phase is from 0.3333 to 0.6666, and radiation factor is taken to be 2.All three phases of the HTS algorithm are employed as basic search technic in the proposed multiobjective heat transfer search (MOHTS) algorithm. Notice that the proposed MOHTS algorithm uses the values of the constants CDF, CVF, and RF as used in [25]. Readers are advised to refer to [25] for clear justification behind using these particular values.


### 3.2. Elitist Nondominated Sorting and Crowding Distance Approach

Elitist nondominated sorting method and diversity preserving crowding distance approach of NSGA-II are introduced in the proposed MOHTS algorithm for sorting of the population in different nondomination levels with computed crowded distance. Firstly for each solution obtained from the basic search method (HTS) or from initially generated random populationP0, all the objectives from the objective vector F are evaluated. In addition, a domination count np defined as number of solutions dominating the solution p and Sp which is a set of solutions dominated by solution p are calculated. Secondly, all the solutions p are assigned a domination count zero and are put in first nondominated level also known as Pareto front (PF) and their nondomination rank (NDRp) is set to 1. Thirdly, for each solution p with np=0, each member q of the set Sp is visited and its domination count nq is reduced by one. While reducing nq count if it falls to zero the corresponding solution q is put in second nondomination level and NDRq is set to 2. The procedure is repeated for each member of second nondomination level to obtain the third nondomination level, and subsequently the procedure should be repeated until the whole population is sorted into different nondomination levels.In crowding distance approach for maintaining diversity among the obtained solutions firstly the population is sorted according to value of each objective function in ascending order. An infinite crowding distance is then assigned to the boundary solutions,i=1 and i=l, of each objective. Here l is the total number of solutions in a particular nondominated set. The boundary solutions are the minimum (i=1) and maximum (i=l) function values. Except the boundary solutions, all the other solutions of the sorted population (i=2tol-1) for each objective j (j=1,2,…,m) are assigned the crowding distance (dmi) as(15)dji=fji+1-fji-1fjmax-fjmin.In (15) the right hand side term is the difference in values of objective function j for two neighbouring solutions (i+1 and i-1) of solution i. Now each solution i is assigned two entities, nondomination rank NDRi and crowding distance CDi. A crowded comparison operator (≺n) is used as follows to compare two solutions (i and j) i≺nj, if (NDRi<NDRj) or ((NDRi=NDRj) and (CDi>CDj)). That is to say, between two solutions, the one with the lower nondomination rank is preferred and if both the solutions have the same nondomination rank then the one with the higher crowding distance is preferred.


### 3.3. MOHTS Procedure

The procedure of the proposed MOHTS algorithm has been shown in Pseudocode1. Firstly the parameters such as population size (Npop), termination criteria, here being the maximum number of generations (Ngen), conduction, convection, and radiation factors (CDF, CVF, and RF, resp.) are initialized. Secondly a random parent population Po in feasible region S is generated and each objective function of the objective vector F for Po is evaluated. Next, elitist based nondominated sorting and crowding distance computation as explained in earlier section is applied on Po. Thirdly HTS algorithm is employed to create the offspring population Pj, which is then merged with Po to form the merged parent-offspring population Pi. This Pi is sorted based on elitism nondomination, and based on the computed values of NDRp and CDp the best Npop solutions are updated to form a new parent population. This process is repeated till the maximum number of generations (iterations) is reached. It should be noted that the same algorithm can also be used with the termination criteria set on the basis of number of function evaluations. Since nondominated sorting and crowding distance assignment of MOHTS area is adopted from NSGA-II, the computational complexity of MOHTS is also the same as NSGA-II which is OmNpop2, where m is the total number of objective functions and Npop is the population size.Pseudocode 1:Pseudocode of MOHTS algorithm.
Initialize
N
p
o
p
,
N
g
e
n
, CDF, CVF, RF and Rand
Generate
P
0
randomly in S and evaluate F for the generated
P
0
Sort the
P
0
based on elitist non dominated sort method and find the NDR and fronts
Compute the CD for each front
For
i
=
1
to
N
g
e
n
// MOHTS main loop starts
For
j
=
1
to
N
p
o
p
If Rand ≤ 0.3333           // Conduction phase
Get
P
j
from equations (5)–(8)
Else if 0.3333 < Rand ≤ 0.6666       // Radiation phase
Get
P
j
from equations (11)–(14)
Else                 // Convection phase
Get
P
j
from equations (9)-(10)
End if
End for
Merge
P
0
and
P
j
to create
P
i,i.e.Pi  =  P0  ∪  Pj
For
P
i
apply the elitist non dominated sort method and find the NDR of the fronts
Apply the crowding distance for the sorted set of solutions for each objective
Based on NDR and CD sort
P
i
in descending order
Replace
P
0
with
P
i
for first
N
p
o
p
members of
P
i
End for                  // MOHTS main loop ends


## 3.1. Basic Heat Transfer Search Algorithm

A very recently developed metaheuristic algorithm called heat transfer search (HTS) algorithm [25] which is based on the thermodynamic mechanism of the heat transfer has been employed as a basic search strategy in the proposed MOHTS algorithm. The search agents in the algorithm are the molecules which interact with one another and with the surrounding for gaining the thermal equilibrium. Unlike other metaheuristic algorithms, the HTS has been developed by putting more efforts in setting a well-balanced trade-off between the intensification (exploitation) and diversification (exploration). The algorithm is divided into three phases as conduction, convection, and radiation; and in the course of entire search process each phase is performed with an equal probability. All the three phases intensify the search space in initial iterations and later diversify during remaining iterations. The following sections briefly explain all three phases.


### 3.1.1. Conduction Phase

As stated above in HTS algorithm the population of solutions is represented as number of molecules, design variables are represented as molecules’ temperature level, and fitness values are represented as energy level of molecules. In this phase energy from the molecules with high energy is transferred to the ones with the lower energy. The population is updated during this phase and can be described as follows:If number of iterationsn ≤ maximum number of iterations/CDF (initial populations) andiffXm<fXp then(5)Xp,lnew=1-R2Xm,lold.ElsefXp<fXm; then(6)Xm,lnew=1-R2Xp,lold.Else, for the remaining populations (number of iterationsn> maximum number of iterations/CDF),iffXm<fXp then(7)Xp,lnew=1-rlXm,lold.ElsefXp<fXm; then(8)Xm,lnew=1-rlXp,lold.Here,m and p indicate randomly picked solutions and vary from 1 to u (population size); and l indicates the design variable index ranging from 1 to v (total number of design variables), which is also selected randomly, and CDF is the conduction factor. Moreover, R and rl are probability value for conduction phase (this value is between 0 and 0.3333) and random number between 0 and 1, respectively. For the first half of the total generations (CDF should be taken as 2), a function evaluation comparison is made between two randomly selected solutions; the inferior of the two is replaced by (5) and (6). For the latter half the inferior solution is replaced by (7) and (8). Thus the switching from (5) and (6) to (7) and (8) obtains both the intensification and diversification for maintaining the algorithm’s potential to seek an optimum solution.


### 3.1.2. Convection Phase

In this phase, a thermal equilibrium is attained as system’s mean temperature (Xmean) interacts with the surrounding temperature (Xs). The latter temperature is assumed to be the best solution and the solution is updated as follows:If current iteration numbern≤ maximum number of iterations/CVF (initial populations) then(9)Xm,lnew=Xm,lold+RXs-Xmean·absR-rl.Else, for the remaining population,(10)Xm,lnew=Xm,lold+RXs-Xmean·round1+rl.It should be noted that, along with each solutionXm, each design variable is also updated during this phase. Here R is taken in the range from 0.6666 to 1, and CVF (convection factor) is taken as 10. That is to say, for the first one-tenth of the iterations, the algorithm explores the feasible design space and for the remaining iterations it exploits the same, but in actual practice this condition does not take place as exploration (diversification) and exploitation (intensification) are taken care of by the last terms of (9) and (10).


### 3.1.3. Radiation Phase

In this phase, system (solution) communicates either with the surroundings which is the best solution or within the system which is some other solution to attain the thermal equilibrium. The solutions are improved as described below.If number of iterationsn≤ maximum number of iterations/RF (initial populations) andiffXm<fXp then(11)Xm,lnew=1+RXm,lold-R·Xp,lold.ElsefXp<fXm; then(12)Xm,lnew=1-RXm,lold+R·Xp,lold.Else, for the remaining populations,iffXm<fXp then(13)Xm,lnew=1+rlXm,lold-rl·Xp,lold.ElsefXp<fXm; then(14)Xm,lnew=1-rlXm,lold+rl·Xp,lold.Herem and p are randomly selected solutions for which their fitness values are compared with each other and that determines how each solution is updated using (11), (12), (13), and (14). The range of R in this phase is from 0.3333 to 0.6666, and radiation factor is taken to be 2.All three phases of the HTS algorithm are employed as basic search technic in the proposed multiobjective heat transfer search (MOHTS) algorithm. Notice that the proposed MOHTS algorithm uses the values of the constants CDF, CVF, and RF as used in [25]. Readers are advised to refer to [25] for clear justification behind using these particular values.


## 3.1.1. Conduction Phase

As stated above in HTS algorithm the population of solutions is represented as number of molecules, design variables are represented as molecules’ temperature level, and fitness values are represented as energy level of molecules. In this phase energy from the molecules with high energy is transferred to the ones with the lower energy. The population is updated during this phase and can be described as follows:If number of iterationsn ≤ maximum number of iterations/CDF (initial populations) andiffXm<fXp then(5)Xp,lnew=1-R2Xm,lold.ElsefXp<fXm; then(6)Xm,lnew=1-R2Xp,lold.Else, for the remaining populations (number of iterationsn> maximum number of iterations/CDF),iffXm<fXp then(7)Xp,lnew=1-rlXm,lold.ElsefXp<fXm; then(8)Xm,lnew=1-rlXp,lold.Here,m and p indicate randomly picked solutions and vary from 1 to u (population size); and l indicates the design variable index ranging from 1 to v (total number of design variables), which is also selected randomly, and CDF is the conduction factor. Moreover, R and rl are probability value for conduction phase (this value is between 0 and 0.3333) and random number between 0 and 1, respectively. For the first half of the total generations (CDF should be taken as 2), a function evaluation comparison is made between two randomly selected solutions; the inferior of the two is replaced by (5) and (6). For the latter half the inferior solution is replaced by (7) and (8). Thus the switching from (5) and (6) to (7) and (8) obtains both the intensification and diversification for maintaining the algorithm’s potential to seek an optimum solution.


## 3.1.2. Convection Phase

In this phase, a thermal equilibrium is attained as system’s mean temperature (Xmean) interacts with the surrounding temperature (Xs). The latter temperature is assumed to be the best solution and the solution is updated as follows:If current iteration numbern≤ maximum number of iterations/CVF (initial populations) then(9)Xm,lnew=Xm,lold+RXs-Xmean·absR-rl.Else, for the remaining population,(10)Xm,lnew=Xm,lold+RXs-Xmean·round1+rl.It should be noted that, along with each solutionXm, each design variable is also updated during this phase. Here R is taken in the range from 0.6666 to 1, and CVF (convection factor) is taken as 10. That is to say, for the first one-tenth of the iterations, the algorithm explores the feasible design space and for the remaining iterations it exploits the same, but in actual practice this condition does not take place as exploration (diversification) and exploitation (intensification) are taken care of by the last terms of (9) and (10).


## 3.1.3. Radiation Phase

In this phase, system (solution) communicates either with the surroundings which is the best solution or within the system which is some other solution to attain the thermal equilibrium. The solutions are improved as described below.If number of iterationsn≤ maximum number of iterations/RF (initial populations) andiffXm<fXp then(11)Xm,lnew=1+RXm,lold-R·Xp,lold.ElsefXp<fXm; then(12)Xm,lnew=1-RXm,lold+R·Xp,lold.Else, for the remaining populations,iffXm<fXp then(13)Xm,lnew=1+rlXm,lold-rl·Xp,lold.ElsefXp<fXm; then(14)Xm,lnew=1-rlXm,lold+rl·Xp,lold.Herem and p are randomly selected solutions for which their fitness values are compared with each other and that determines how each solution is updated using (11), (12), (13), and (14). The range of R in this phase is from 0.3333 to 0.6666, and radiation factor is taken to be 2.All three phases of the HTS algorithm are employed as basic search technic in the proposed multiobjective heat transfer search (MOHTS) algorithm. Notice that the proposed MOHTS algorithm uses the values of the constants CDF, CVF, and RF as used in [25]. Readers are advised to refer to [25] for clear justification behind using these particular values.


## 3.2. Elitist Nondominated Sorting and Crowding Distance Approach

Elitist nondominated sorting method and diversity preserving crowding distance approach of NSGA-II are introduced in the proposed MOHTS algorithm for sorting of the population in different nondomination levels with computed crowded distance. Firstly for each solution obtained from the basic search method (HTS) or from initially generated random populationP0, all the objectives from the objective vector F are evaluated. In addition, a domination count np defined as number of solutions dominating the solution p and Sp which is a set of solutions dominated by solution p are calculated. Secondly, all the solutions p are assigned a domination count zero and are put in first nondominated level also known as Pareto front (PF) and their nondomination rank (NDRp) is set to 1. Thirdly, for each solution p with np=0, each member q of the set Sp is visited and its domination count nq is reduced by one. While reducing nq count if it falls to zero the corresponding solution q is put in second nondomination level and NDRq is set to 2. The procedure is repeated for each member of second nondomination level to obtain the third nondomination level, and subsequently the procedure should be repeated until the whole population is sorted into different nondomination levels.In crowding distance approach for maintaining diversity among the obtained solutions firstly the population is sorted according to value of each objective function in ascending order. An infinite crowding distance is then assigned to the boundary solutions,i=1 and i=l, of each objective. Here l is the total number of solutions in a particular nondominated set. The boundary solutions are the minimum (i=1) and maximum (i=l) function values. Except the boundary solutions, all the other solutions of the sorted population (i=2tol-1) for each objective j (j=1,2,…,m) are assigned the crowding distance (dmi) as(15)dji=fji+1-fji-1fjmax-fjmin.In (15) the right hand side term is the difference in values of objective function j for two neighbouring solutions (i+1 and i-1) of solution i. Now each solution i is assigned two entities, nondomination rank NDRi and crowding distance CDi. A crowded comparison operator (≺n) is used as follows to compare two solutions (i and j) i≺nj, if (NDRi<NDRj) or ((NDRi=NDRj) and (CDi>CDj)). That is to say, between two solutions, the one with the lower nondomination rank is preferred and if both the solutions have the same nondomination rank then the one with the higher crowding distance is preferred.


## 3.3. MOHTS Procedure

The procedure of the proposed MOHTS algorithm has been shown in Pseudocode1. Firstly the parameters such as population size (Npop), termination criteria, here being the maximum number of generations (Ngen), conduction, convection, and radiation factors (CDF, CVF, and RF, resp.) are initialized. Secondly a random parent population Po in feasible region S is generated and each objective function of the objective vector F for Po is evaluated. Next, elitist based nondominated sorting and crowding distance computation as explained in earlier section is applied on Po. Thirdly HTS algorithm is employed to create the offspring population Pj, which is then merged with Po to form the merged parent-offspring population Pi. This Pi is sorted based on elitism nondomination, and based on the computed values of NDRp and CDp the best Npop solutions are updated to form a new parent population. This process is repeated till the maximum number of generations (iterations) is reached. It should be noted that the same algorithm can also be used with the termination criteria set on the basis of number of function evaluations. Since nondominated sorting and crowding distance assignment of MOHTS area is adopted from NSGA-II, the computational complexity of MOHTS is also the same as NSGA-II which is OmNpop2, where m is the total number of objective functions and Npop is the population size.Pseudocode 1:Pseudocode of MOHTS algorithm.
Initialize
N
p
o
p
,
N
g
e
n
, CDF, CVF, RF and Rand
Generate
P
0
randomly in S and evaluate F for the generated
P
0
Sort the
P
0
based on elitist non dominated sort method and find the NDR and fronts
Compute the CD for each front
For
i
=
1
to
N
g
e
n
// MOHTS main loop starts
For
j
=
1
to
N
p
o
p
If Rand ≤ 0.3333           // Conduction phase
Get
P
j
from equations (5)–(8)
Else if 0.3333 < Rand ≤ 0.6666       // Radiation phase
Get
P
j
from equations (11)–(14)
Else                 // Convection phase
Get
P
j
from equations (9)-(10)
End if
End for
Merge
P
0
and
P
j
to create
P
i,i.e.Pi  =  P0  ∪  Pj
For
P
i
apply the elitist non dominated sort method and find the NDR of the fronts
Apply the crowding distance for the sorted set of solutions for each objective
Based on NDR and CD sort
P
i
in descending order
Replace
P
0
with
P
i
for first
N
p
o
p
members of
P
i
End for                  // MOHTS main loop ends


## 4. Multiobjective Optimization of Half Car Ride Model

A half car model of the vehicle’s passive suspension system with five degrees of freedom has been adopted from the literature by Bouazara [23] and is shown in Figure 1. The lumped mass model consists of driver’s seat mass mc, a sprung mass ms, and two unsprung masses as front and rear tires m1 and m2, respectively. The elastic and damping properties of the driver’s seat are represented as stiffness and damping coefficients kss and css, respectively; similarly elastic properties of the front and rear tires are depicted as stiffness coefficients kp1 and kp2, respectively. The damping properties of the tires are assumed to be minimal and thus are neglected. Sprung mass ms along with mc rests on the front and rear suspensions. The stiffness and damping coefficients of front suspension are ks1 and cs1, respectively, and of rear suspension are ks2 and cs2, respectively. In addition to the four heave (vertical linear) degrees of freedom of four lumped masses, a rotational degree of freedom (pitch motion) has been considered as a fifth degree of freedom of the system. In Figure 1, Is, r, l1, and l2 are mass moment of inertia of sprung mass, position of the driver’s seat, and positions of front and rear tires with respect to the centre of mass of the vehicle, respectively. In addition, the heave or vertical motion of the vehicle is assigned by the vertical displacements z with subscripts showing displacement of the associate lumped masses. The rotational motion (pitch) with respect to the axis perpendicular to the vehicle’s plane and passing through the centre of mass of the vehicle is shown by θ. In Bouazara [23] model, the variables g1 and g2 are the damping coefficients for active suspension system. Since the present study has adopted this model as a passive half car model, these coefficients are neglected.Figure 1
Half car passive ride model with five degrees of freedom, adopted from [23].Applying the Newton-Euler equations to the multibody spring-mass-damper model leads to the equations of motion of the half car model which are given by(16)zps=zs-θzs1=zs-l1θzs2=zs+l2θFss=ksszc-zps+cssz˙c-z˙psFs1=ks1zs1-z1+cs1z˙s1-z˙1Fs2=ks2zs2-z2+cs2z˙s2-z˙2mcz¨c=-Fssmcz¨s=-Fs1-Fs2+FssIsθ¨s=l1Fs1-l2Fs2+rFssm1z¨1=Fs1-kp1z1-zp1m2z¨2=Fs2-kp2z2-zp2.The equations of motion given by (16) are a set of five, second-order ordinary differential equations. Here, zc, zs, z1, z2, and θ are the vertical linear displacements of the driver’s seat, the gravitational centre of the sprung mass, front tire, and rear tire, and the angular displacement of the sprung mass with respect to its centre of mass, respectively. In addition, zs1, zs2, zp1, and zp2 depict the vertical displacements of front and rear end of the sprung mass and the road disturbances on the front and rear tires, respectively.Practically, the road disturbances are random in nature and thus for simulating the mathematical models they are treated as Gaussian isotropic random field [27], stationary and ergodic stochastic process, or are modeled as single obstacles such as speed bumps and potholes [24, 28, 29]. Single obstacle is a simple way to model the road excitation; the double bump excitation as used by Bouazara [23] has been used in the present study and is as shown in Figure 2. The equations of motion of the half car were modeled in Simulink/MATLAB and were solved using fourth-order Runge-Kutta numerical method (RK4) using MATLAB.Figure 2
Double bump road excitation.The multiobjective optimization of the half car model using the proposed MOHTS and NSGA-II is presented in two numerical studies. The first study adopts the optimization problem from [21, 22] with a total of five objectives and seven design variables, whereas the second numerical study considers the more realistic approach in choosing the objective functions for optimizing the passive half car model and the same design variables as used in study 1. The seven design variables taken for the optimization studies are stiffness and damping coefficients of driver’s seat (kss, css), front suspension (ks1, cs1), and rear suspension (ks2, cs2) and the position of the driver’s seat from the centre of mass of the vehicle (r). The rest of all the parameters of the half car model as shown in Figure 1 are considered as vehicle constants and are the same for both the numerical studies and are given in Table 1. Furthermore, the bounds on the design variables are also kept the same for both studies and are as specified in Table 2. The overall procedure followed for optimizing and simulating the half car ride model is shown in Figure 3.Table 1
Vehicle constants.
Parameter
Value
Unit
Parameter
Value
Unit
m
c
75
kg
l
1
1.011
m
m
s
730
kg
l
2
1.803
m
I
s
1230
kg m2
k
p
1
175500
N/m
m
1
40
kg
k
p
2
175500
N/m
m
2
35.5
kg
g
1, g2
0
N s/mTable 2
Design variable bounds.
Design variable
Lower bound
Upper bound
k
s
s (N/m)
50000
150000
c
s
s (N s/m)
1000
4000
k
s
1, ks2 (N/m)
10000
20000
C
s
1, Cs2 (N s/m)
500
2000
r (m)
0
0.5Figure 3
Optimization and simulation procedure followed in optimization of the half car ride model.


### 4.1. Numerical Study 1

The multiobjective optimization problem’s statement for study 1 is given as follows.Find the design vector,X=ksscssks1cs1ks2cs2rT, which minimizes(17)FX=z¨cz˙1z˙2zs-z1zs-z2Twith no constraints, and the bounds on the design variables are as presented in Table 2. The primary objectives of the suspension system are explained in Section 1 and are developed as measured performance criteria of the suspension system as given in (17). The first objective in (17) is that maximum ride comfort can be obtained by minimizing the absolute vertical linear acceleration of the driver’s seat. The fourth and fifth objectives in (17) ensure minimum suspension rattle space and can be achieved by minimizing the front and rear suspension deflections, respectively. Notice that second and third objectives in (17) do not fulfill any primary objective of the suspension system. These are only chosen for study 1 merely to be consistent in the comparison of the proposed MOHTS algorithm with other algorithms which published the results based on the objective function vector as given by (17). More importantly all five selected objectives take the absolute values of their respective quantities and are conflicting with each other; hence they suite the multiobjective optimization quite effectively.Firstly, four 2-objective optimization studies, as taken by Jamali et al. [21], are executed with the pairs of objectives as z¨c,z˙1, z¨c,z˙2, z¨c,zs-z1, and z¨c,zs-z2 using MOHTS and NSGA-II. For both MOHTS and NSGA-II, the population size is selected as 80, and the total number of generations is set to be 240. Furthermore for NSGA-II the crossover and mutation probabilities were taken as 0.9 and 0.1, respectively. The obtained Pareto front for the 2-objective optimization studies is presented in Figures 4–7. In the figures, points Ai and Bi are the extreme left and bottom trade-off points, respectively, of the Pareto front given by MOHTS and points Ci and Di are the extreme left and bottom trade-off points, respectively, of the Pareto front given by NSGA-II. In addition, points Ei and Fi refer to the intermediate trade-off points given by MOHTS and NSGA-II, respectively. These trade-off points are different from the extreme points but show significant design information which would be discussed individually later in this section.Figure 4
Pareto front of seat acceleration and front tire velocity in 2-objective optimization using MOHTS and NSGA-II.Figure 5
Pareto front of seat acceleration and rear tire velocity in 2-objective optimization using MOHTS and NSGA-II.Figure 6
Pareto front of seat acceleration and front suspension deflection in 2-objective optimization using MOHTS and NSGA-II.Figure 7
Pareto front of seat acceleration and rear suspension deflection in 2-objective optimization using MOHTS and NSGA-II.As can be seen from Figure4, overall MOHTS has produced a Pareto front with the better diversity among the optimal solutions and with the best nondominated set of solutions compared to that of NSGA-II. Extreme left trade-off point of MOHTS A1 has lower seat acceleration than its equivalent point C1 of NSGA-II but at the cost of higher front tire velocity. Similarly extreme bottom trade-off point B1 of MOHTS has lower front tire velocity than its match point D1 of NSGA-II, though in expense of higher seat acceleration. In addition to the extreme points on the Pareto front, the intermediate trade-off points are equally important. For instance, the optimal point of MOHTS E1 has lower seat acceleration than any optimal point of NSGA-II and lower front tire velocities than many optimal points of NSGA-II. Even more E1 registers 48.5% lower seat acceleration than B1 and 24.65% lower front tire velocity than A1. Likewise, F1 of NSGA-II has 44.34% lower seat acceleration than D1 and 4.12% lower front tire velocity than C1. The points E1 and F1 would not have been obtained without using Pareto based multiobjective optimization algorithms; thus such algorithms are able to help the decision maker by giving alternative design trade-off points. Similar to points E1 and F1, improvements in other intermediate trade-off points in comparison with their respective Pareto fronts’ extreme trade-off points are demonstrated in Table 3. It should be noted that the chosen intermediate trade-off points have the highest average percentage improvement (as compared to the rest of the intermediate trade-off points) in both the objectives with respect to the extreme worst points of the same objective. The dashes in Table 3 show that the chosen intermediate trade-off point has no relevance to the selected objective function and hence no improvement or deterioration can be measured.Table 3
Design information of intermediate trade-off points.
Trade-off point
Percentage improvement (positive) or deterioration (negative) of trade-off point in
In comparison with
z
¨
c
z
˙
1
z
˙
2
z
s
-
z
1
z
s
-
z
2
E
1
−77.08
24.65
—
—
—
A
1
48.5
−0.44
—
—
—
B
1
F
1
−0.17
4.12
—
—
—
C
1
44.34
−4.39
—
—
—
D
1
E
2
−17.24
—
1.53
—
—
A
2
65.82
—
−1.23
—
—
B
2
F
2
−2.66
—
7.58
—
—
C
2
6.59
—
−0.18
—
—
D
2
E
3
−12.37
—
—
9.28
—
A
3
5.73
—
—
−6.23
—
B
3
F
3
−15.48
—
—
10.49
—
C
3
3.25
—
—
−4.61
—
D
3
E
4
−20.27
—
—
—
39.81
A
4
0.26
—
—
—
−0.39
B
4
F
4
−15.53
—
—
—
31.99
C
4
0.77
—
—
—
−0.8
D
4Furthermore, as can be observed in Figures5, 6, and 7, MOHTS also yields better nondominated set of Pareto optimal solutions with better and wider spread of the optimal solutions than NSGA-II in the rest of the three 2-objective optimization studies.The comparison of the extreme trade-off points given by the proposed algorithm MOHTS, NSGA-II, MUGA [21], and a combined PSO-GA based MOEA [22] is presented in Tables 4–7. The points Ai′ and Bi′ are the extreme left and bottom points given by MUGA, and points Ai′′ and Bi′′ are the extreme left and bottom points given by combined PSO-GA based MOEA. As can be observed from the last two columns of the tables, MOHTS either outperforms or performs similar to NSGA-II but clearly outperforms the previous work by Bouazara [23], Jamali et al. [21], and Mahmoodabadi et al. [22]. Since Bouazara [23] decomposed the multiobjective optimization into single objective optimization by using the weighted coefficients, his study yielded a single optimum solution.Table 4
The comparison of the extreme trade-off points (as shown in Figure3) in terms of the values of objective functions and their corresponding design variables between present work (MOHTS), NSGA-II, and previous work by Bouazara [23] Jamali et al. [21], and Mahmoodabadi et al. [22].
Trade-off point
k
s
s (N/m)
k
s
1 (N/m)
k
s
2 (N/m)
c
s
s (Ns/m)
c
s
1 (Ns/m)
c
s
2 (Ns/m)
r (m)
z
¨
c (m/s2)
z
˙
1 (m/s)
A
1  (MOHTS)
77448.364
10000
19940.983
3991.464
500
1997.748
0.5
1.831
0.495
B
1  (MOHTS)
61311.317
20000
20000
4000
2000
1904.71
0.047
4.06
0.358
C
1  (NSGA-II)
146432.050
10000
10000
1369.147
1739.995
1798.133
0.415
2.129
0.388
D
1  (NSGA-II)
81907.529
19999.348
17542.899
3517.244
2000
1767.371
0.496
3.832
0.358
A
1
′  from [21]
146863
10039.2
10313.7
3941.18
1917.65
1958.82
0.5
2.807
0.417
B
1
′  from [21]
57451
10039.2
18196.1
2517.65
964.706
1476.47
0.3412
3.297
0.407
A
1
′
′  from [22]
146728.6
10006.86
10020.65
3964.849
1806.955
1947.486
0.498
2.794
0.415
B
1
′
′  from [22]
149700.2
10003.23
19907.76
1796.857
940.0627
548.968
0.499
3.331
0.405
P  from [24]
105630
14958
14985
2548
1956
1967
0.288
3.089
0.375Table 5
The comparison of the extreme trade-off points (as shown in Figure4) in terms of the values of objective functions and their corresponding design variables between present work (MOHTS), NSGA-II, and previous work by Bouazara [23], Jamali et al. [21], and Mahmoodabadi et al. [22].
Trade-off point
k
s
s (N/m)
k
s
1 (N/m)
k
s
2 (N/m)
c
s
s (Ns/m)
c
s
1 (Ns/m)
c
s
2 (Ns/m)
r (m)
z
¨
c (m/s2)
z
˙
2 (m/s)
A
2  (MOHTS)
76695.747
10000
19925.858
4000
500
1994.622
0.5
1.831
0.361
B
2  (MOHTS)
149997.760
18093.905
19998.477
1586.362
500
1998.436
0
6.282
0.351
C
2  (NSGA-II)
146158.890
10000.508
10000
2283.425
1776.585
1905.497
0.415
2.102
0.386
D
2  (NSGA-II)
122584.350
10159.024
19999.999
3893.357
1988.289
1999.973
0.270
2.311
0.357
A
2
′  from [21]
146863
10039.2
10313.7
3941.18
1917.65
1958.82
0.5
2.807
0.434
B
2
′  from [21]
52745.1
10000
10000
3635.249
670.588
1658.82
0.058
3.475
0.430
A
2
′
′  from [22]
131574.7
10027.98
10107.79
3292.868
1906.304
1935.256
0.478
2.818
0.433
B
2
′
′  from [22]
50504.37
11852.13
10000.4
3737.814
502.1251
1112.495
0.003
5.1074
0.428
P  from [24]
105630
14958
14985
2548
1956
1967
0.288
3.089
0.377Table 6
The comparison of the extreme trade-off points (as shown in Figure5) in terms of the values of objective functions and their corresponding design variables between present work (MOHTS), NSGA-II, and previous work by Bouazara [23], Jamali et al. [21], and Mahmoodabadi et al. [22].
Trade-off point
k
s
s (N/m)
k
s
1 (N/m)
k
s
2 (N/m)
c
s
s (Ns/m)
c
s
1 (Ns/m)
c
s
2 (Ns/m)
r (m)
z
¨
c (m/s2)
Z
S
-
Z
1 (m/s)
A
3  (MOHTS)
77446.455
10000
19982.171
3996.929
500.853
1998.217
0.5
1.83
0.079
B
3  (MOHTS)
92412.553
10000
20000
4000
2000
2000
0.5
2.182
0.068
C
3  (NSGA-II)
77397.693
10005.791
18820.197
3997.672
500
1914.903
0.5
1.867
0.08
D
3  (NSGA-II)
107565.51
10147.906
19940.144
3604.317
2000
1999.218
0.466
2.229
0.068
A
3
′  from [21]
146863
10039.2
10313.7
3941.18
1917.65
1958.82
0.5
2.807
0.09
B
3
′  from [21]
57058.8
10000
10000
3588.24
2000
1676.47
0.114
3.045
0.08
A
3
′
′  from [22]
147769.6
10012.38
10037.35
3926.99
1997.283
1879.508
0.498
2.796
0.0896
B
3
′
′  from [22]
50006.19
10000.09
10000.09
3999.446
1999.998
785.631
4.28
E
-
05
3.605
0.0866
P  from [24]
105630
14958
14985
2548
1956
1967
0.288
3.089
0.0844Table 7
The comparison of the extreme trade-off points (as shown in Figure6) in terms of the values of objective functions and their corresponding design variables between present work (MOHTS), NSGA-II, and previous work by Bouazara [23], Jamali et al. [21], and Mahmoodabadi et al. [22].
Trade-off point
k
s
s (N/m)
k
s
1 (N/m)
k
s
2 (N/m)
c
s
s (Ns/m)
c
s
1 (Ns/m)
c
s
2 (Ns/m)
r (m)
z
¨
c (m/s2)
Z
S
-
Z
2 (m/s)
A
4  (MOHTS)
78255.336
10000
19984.602
3996.583
502.137
1998.292
0.5
1.832
0.056
B
4  (MOHTS)
104827.01
10000
20000
4000
1999.516
2000
0.5
2.209
0.034
C
4  (NSGA-II)
92998.785
10002.354
19987.524
3995.824
601.103
1780.096
0.478
1.945
0.051
D
4  (NSGA-II)
114107.32
10354.485
19910.39
3637.269
2000
1983.220
0.485
2.264
0.035
A
4
′  from [21]
146863
10039.2
10313.7
3941.18
1917.65
1958.82
0.5
2.807
0.06
B
4
′  from [21]
92745.1
10000
18313.7
2941.18
576.471
1994.12
0.484
3.551
0.052
A
4
′
′  from [22]
147107.3
10001.45
10020.01
3941.79
1843.581
1999.742
0.499
2.791
0.059
B
4
′
′  from [22]
149232.8
10001.75
19998.5
3959.178
500.222
1999.908
0.499
3.655
0.05
P  from [24]
105630
14958
14985
2548
1956
1967
0.288
3.089
0.054The response of the seat acceleration, front suspension deflection, and rear suspension deflection of the Pareto optimal design points obtained from MOHTS and NSGA-II and that of the optimal point found by Bouazara [23] are as shown in Figures 8–10, respectively. Figure 8 shows the seat acceleration response of the half car ride model for the design points A3 of MOHTS, C3 of NSGA-II, and P of the work by Bouazara [23]. Figure 9 depicts the front suspension deflection response of the system for the design points which yield the minimum front suspension deflection, and these points are B3 of MOHTS, D3 of NSGA-II, and P of the work by Bouazara [23]. Similarly rear suspension deflection response of the system for the optimal points B4 of MOHTS, D4 of NSGA-II, and P of the work by Bouazara [23] is demonstrated in Figure 10.Figure 8
Seat acceleration response of the half car model for different optimal design points.Figure 9
Front suspension deflection response of the half car model for different optimal design points.Figure 10
Rear suspension deflection response of the half car model for different optimal design points.Lastly, 5-objective optimization study, which seeks the Pareto front for simultaneous optimization of all five objectives of (17), with the same algorithm parameters of 2-objective optimization study results, is investigated. Since 5-objective optimization study deals with the five-dimensional objective space and hence for better understanding of the obtained Pareto front, we use the idea proposed by [30] of using biobjective cross sections (slices) of the Pareto front. Such slices are known as decision maps and they give a clear picture of trade-offs between all the objectives. Such decision maps help the decision maker to better understand the obtained Pareto front. However this approach does not take into account the instability of the Pareto front. The same four pairs of 2-objective optimization as presented in study 1 are presented in study 2 as well with their original Pareto front being superimposed over the Pareto front of 5-objective optimization study. Figures 11–14 are the various planes of the Pareto front given by MOHTS and Figures 15–18 are the same planes of the Pareto front given by NSGA-II. In comparison with the 2-objective optimization study, 5-objective optimization study offers more choices for the decision maker. Figure 11 demonstrates the plane of z¨c,z˙1 from the 5-objective optimization study obtained by the proposed MOHTS algorithm and as can be seen Pareto front is widely spread as compared to 2-objective optimization Pareto front for the same two objectives. This 5-objective Pareto front of Figure 11 has also been depicted in Figures 12, 13, and 14 with a view from the planes z¨c,z˙2, z¨c,zs-z1, and z¨c,zs-z2, respectively. Similarly Figures 15–18 show the Pareto fronts of 5-objective optimization study with their corresponding 2-objective Pareto front superimposed in the planes of z¨c,z˙1, z¨c,z˙2, z¨c,zs-z1, and z¨c,zs-z2, respectively. In summary, the comparison of the proposed MOHTS algorithm with NSGA-II indicates that MOHTS has produced a widespread and better Pareto front than NSGA-II. We ran both MOHTS and NSGA-II for 10 runs and computed the average time taken to solve the optimization problem in study 1. MOHTS and NSGA-II took 1567 and 1364 seconds on average, respectively, showing NSGA-II edges over MOHTS in computational time required to solve the numerical problem. Apart from the Pareto front and computational time, parameter or performance sensitivity could be selected as a performance comparison criterion between MOHTS and NSGA-II. However tuning the algorithm parameters to obtain the best possible outcomes requires exhaustive studies and thus to be fair in making the comparison we choose the same mutation and crossover probabilities for both the algorithms as indicated earlier in this section.Figure 11
Pareto front of seat acceleration and front tire velocity in 5-objective optimization and 2-objective optimization using MOHTS for numerical study 1.Figure 12
Pareto front of seat acceleration and rear tire velocity in 5-objective optimization and 2-objective optimization using MOHTS for numerical study 1.Figure 13
Pareto front of seat acceleration and front suspension deflection in 5-objective optimization and 2-objective optimization using MOHTS for numerical study 1.Figure 14
Pareto front of seat acceleration and rear suspension deflection in 5-objective optimization and 2-objective optimization using MOHTS for numerical study 1.Figure 15
Pareto front of seat acceleration and front tire velocity in 5-objective optimization and 2-objective optimization using NSGA-II for numerical study 1.Figure 16
Pareto front of seat acceleration and front tire velocity in 5-objective optimization and 2-objective optimization using NSGA-II for numerical study 1.Figure 17
Pareto front of seat acceleration and front tire velocity in 5-objective optimization and 2-objective optimization using NSGA-II for numerical study 1.Figure 18
Pareto front of seat acceleration and front tire velocity in 5-objective optimization and 2-objective optimization using NSGA-II for numerical study 1.


### 4.2. Numerical Study 2

In contrast to the first numerical study, we choose more realistic approach in developing the objective functions for optimizing the design of a half car model which is considered with the same set of design variables, vehicle constants, and bounds on the design variables as chosen in study 1. As stated earlier, vehicle suspension system’s three opposing performance indices or objectives are used to evaluate the suspension system’s dynamic response. A better method to estimate the ride discomfort is to measure the root mean square (RMS) of the vertical acceleration felt by the passenger/driver. In this study RMS of the seat acceleration is chosen as the first objective. Equivalently road holding is determined through the RMS of tire deflection, and suspension rattle space is evaluated by measuring the RMS of suspension deflection. The RMS of front and rear suspension deflection and RMS of front and rear tire deflections are selected as other four objectives to form the 5-objective optimization problem in this study. The optimization problem is stated as follows.Find the design vector,X=ksscssks1cs1ks2cs2rT, which minimizes(18)FX=f1Xf2Xf3Xf4Xf5XTf1X=1T∫0Tz¨c2dtf2X=1T∫0Tzs-z12dt,f3X=1T∫0Tzs-z22dtf4X=1T∫0Tz1-zp12dt,f5X=1T∫0Tzs-zp22dtwith no constraints.Unlike the previous study, the Pareto fronts obtained from 2-objective studies are superimposed over the corresponding Pareto fronts yielded from the 5-objective studies to shorten the length of the paper. We again use the same slicing method for visualizing the Pareto front in two-dimensional space to better demonstrate the effectiveness of MOHTS and NSGA-II. Figures19–22 show the planes (slices) of f1,f2, f1,f3, f1,f4, and f1,f5, respectively, for the 5-objective optimization study using MOHTS and likewise Figures 23–26 depict the same planes in the same order for the 5-objective optimization using NSGA-II. The objective function and design variable values of the extreme points of the 2-objective studies’ Pareto fronts as shown in Figures 19–26 are given in Table 8. It is clear from the table that MOHTS is superior to NSGA-II in minimizing the objectives of 2-objective optimization studies and the minimum values of each objective are shown in bold letters.Table 8
The comparison of the extreme trade-off points (as shown in Figures18–25) in terms of the values of objective functions and their corresponding design variables between present work (MOHTS) and NSGA-II.
Trade-off point
k
s
s (N/m)
k
s
1 (N/m)
k
s
2 (N/m)
c
s
s (Ns/m)
c
s
1 (Ns/m)
c
s
2 (Ns/m)
r (m)
f1(m/s2)
f
2 (m)
f
3 (m)
f
4 (m)
f
5 (m)
A
5  (MOHTS)
62325.038
10000.000
20000.000
4000.000
500.000
1999.952
0.500
0.850
0.035
0.026
0.003
0.003
B
5  (MOHTS)
105868.150
10000.000
20000.000
4000.000
2000.000
2000.000
0.500
1.062
0.033
0.016
0.003
0.002
C
5  (NSGA-II)
111650.910
10004.447
13795.435
3676.384
1411.603
1970.654
0.434
1.017
0.037
0.022
0.003
0.002
D
5  (NSGA-II)
113125.130
10000.024
13996.260
3668.958
1477.650
1997.922
0.454
1.021
0.036
0.021
0.003
0.002
A
6  (MOHTS)
66552.042
10000.000
20000.000
4000.000
500.000
1967.310
0.500
0.850
0.035
0.025
0.003
0.003
B
6  (MOHTS)
104962.350
10000.000
20000.000
4000.000
2000.000
1979.535
0.500
1.063
0.033
0.016
0.003
0.002
C
6  (NSGA-II)
95623.254
10030.399
15187.577
1680.633
1478.532
1995.701
0.500
1.045
0.036
0.021
0.003
0.002
D
6  (NSGA-II)
132369.180
10409.407
19999.716
2728.990
2000.000
1730.412
0.450
1.138
0.034
0.016
0.003
0.003
A
7  (MOHTS)
63426.812
10000.002
19976.225
3999.546
500.357
1962.724
0.500
0.849
0.035
0.026
0.003
0.003
B
7  (MOHTS)
62028.590
10000.000
20000.000
4000.000
500.000
500.000
0.500
1.584
0.048
0.028
0.002
0.009
C
7  (NSGA-II)
57838.283
10434.762
19752.225
2306.963
736.508
1672.299
0.500
1.019
0.038
0.026
0.003
0.003
D
7  (NSGA-II)
58174.332
10430.336
20000.000
2334.326
838.576
500.000
0.500
1.463
0.045
0.025
0.002
0.007
A
8  (MOHTS)
62635.017
10000.000
20000.000
4000.000
500.000
2000.000
0.500
0.850
0.035
0.026
0.003
0.003
B
8  (MOHTS)
50000.000
20000.000
10000.000
4000.000
2000.000
2000.000
0.486
1.908
0.040
0.026
0.007
0.001
C
8  (NSGA-II)
65778.794
10000.000
18885.969
3710.928
500.052
1932.241
0.500
0.874
0.036
0.026
0.003
0.003
D
8  (NSGA-II)
50308.460
19948.963
10000.000
3521.482
1971.198
1972.075
0.486
1.900
0.040
0.026
0.007
0.001Figure 19
Pareto front of seat acceleration RMS and front suspension deflection RMS in 5-objective optimization and 2-objective optimization using MOHTS for numerical study 2.Figure 20
Pareto front of seat acceleration RMS and rear suspension deflection RMS in 5-objective optimization and 2-objective optimization using MOHTS for numerical study 2.Figure 21
Pareto front of seat acceleration RMS and front tire deflection RMS in 5-objective optimization and 2-objective optimization using MOHTS for numerical study 2.Figure 22
Pareto front of seat acceleration RMS and rear tire deflection RMS in 5-objective optimization and 2-objective optimization using MOHTS for numerical study 2.Figure 23
Pareto front of seat acceleration RMS and front suspension deflection RMS in 5-objective optimization and 2-objective optimization using NSGA-II for numerical study 2.Figure 24
Pareto front of seat acceleration RMS and rear suspension deflection RMS in 5-objective optimization and 2-objective optimization using NSGA-II for numerical study 2.Figure 25
Pareto front of seat acceleration RMS and front tire deflection RMS in 5-objective optimization and 2-objective optimization using NSGA-II for numerical study 2.Figure 26
Pareto front of seat acceleration RMS and rear suspension deflection RMS in 5-objective optimization and 2-objective optimization using NSGA-II for numerical study 2.As can be observed from Figures19 and 23, the Pareto front for the 5-objective optimization is quite similar for both MOHTS and NSGA-II and this is also true for the rest all the planes, thus, showing that the potential of the proposed MOHTS algorithm is as good as NSGA-II in finding the Pareto front solutions of the multiobjective optimization problem. Moreover, when 2-objective optimization Pareto fronts from MOHTS are compared for all the shown pairs of objectives with NSGA-II, it appears that the proposed MOHTS algorithm dominates NSGA-II in obtaining the better nondominated Pareto front with the widespread (diveresed) set of optimal solutions. The reason behind the success of MOHTS in yielding a better Pareto front can be attributed to the search mechanism of HTS, as it takes a good care in maintaining the diversified set of solutions in the search space with its better exploration capabilities, which in turn allows the nondominated sorting and crowding distance approach mechanisms to produce a widespread Pareto front. Furthermore, as exploitation capability has also been wisely implemented in HTS algorithm it helps the dominance sorting and diversity preserving mechanisms in producing a better nondominated set. In contrast to the first numerical study, the average computational time for MOHTS and NSGA-II in this study is 1545 and 1633 seconds, respectively. The reason behind this conflict can be attributed to the uncertainties in the numerical solution of the half car model’s dynamic equations of motion.


## 4.1. Numerical Study 1

The multiobjective optimization problem’s statement for study 1 is given as follows.Find the design vector,X=ksscssks1cs1ks2cs2rT, which minimizes(17)FX=z¨cz˙1z˙2zs-z1zs-z2Twith no constraints, and the bounds on the design variables are as presented in Table 2. The primary objectives of the suspension system are explained in Section 1 and are developed as measured performance criteria of the suspension system as given in (17). The first objective in (17) is that maximum ride comfort can be obtained by minimizing the absolute vertical linear acceleration of the driver’s seat. The fourth and fifth objectives in (17) ensure minimum suspension rattle space and can be achieved by minimizing the front and rear suspension deflections, respectively. Notice that second and third objectives in (17) do not fulfill any primary objective of the suspension system. These are only chosen for study 1 merely to be consistent in the comparison of the proposed MOHTS algorithm with other algorithms which published the results based on the objective function vector as given by (17). More importantly all five selected objectives take the absolute values of their respective quantities and are conflicting with each other; hence they suite the multiobjective optimization quite effectively.Firstly, four 2-objective optimization studies, as taken by Jamali et al. [21], are executed with the pairs of objectives as z¨c,z˙1, z¨c,z˙2, z¨c,zs-z1, and z¨c,zs-z2 using MOHTS and NSGA-II. For both MOHTS and NSGA-II, the population size is selected as 80, and the total number of generations is set to be 240. Furthermore for NSGA-II the crossover and mutation probabilities were taken as 0.9 and 0.1, respectively. The obtained Pareto front for the 2-objective optimization studies is presented in Figures 4–7. In the figures, points Ai and Bi are the extreme left and bottom trade-off points, respectively, of the Pareto front given by MOHTS and points Ci and Di are the extreme left and bottom trade-off points, respectively, of the Pareto front given by NSGA-II. In addition, points Ei and Fi refer to the intermediate trade-off points given by MOHTS and NSGA-II, respectively. These trade-off points are different from the extreme points but show significant design information which would be discussed individually later in this section.Figure 4
Pareto front of seat acceleration and front tire velocity in 2-objective optimization using MOHTS and NSGA-II.Figure 5
Pareto front of seat acceleration and rear tire velocity in 2-objective optimization using MOHTS and NSGA-II.Figure 6
Pareto front of seat acceleration and front suspension deflection in 2-objective optimization using MOHTS and NSGA-II.Figure 7
Pareto front of seat acceleration and rear suspension deflection in 2-objective optimization using MOHTS and NSGA-II.As can be seen from Figure4, overall MOHTS has produced a Pareto front with the better diversity among the optimal solutions and with the best nondominated set of solutions compared to that of NSGA-II. Extreme left trade-off point of MOHTS A1 has lower seat acceleration than its equivalent point C1 of NSGA-II but at the cost of higher front tire velocity. Similarly extreme bottom trade-off point B1 of MOHTS has lower front tire velocity than its match point D1 of NSGA-II, though in expense of higher seat acceleration. In addition to the extreme points on the Pareto front, the intermediate trade-off points are equally important. For instance, the optimal point of MOHTS E1 has lower seat acceleration than any optimal point of NSGA-II and lower front tire velocities than many optimal points of NSGA-II. Even more E1 registers 48.5% lower seat acceleration than B1 and 24.65% lower front tire velocity than A1. Likewise, F1 of NSGA-II has 44.34% lower seat acceleration than D1 and 4.12% lower front tire velocity than C1. The points E1 and F1 would not have been obtained without using Pareto based multiobjective optimization algorithms; thus such algorithms are able to help the decision maker by giving alternative design trade-off points. Similar to points E1 and F1, improvements in other intermediate trade-off points in comparison with their respective Pareto fronts’ extreme trade-off points are demonstrated in Table 3. It should be noted that the chosen intermediate trade-off points have the highest average percentage improvement (as compared to the rest of the intermediate trade-off points) in both the objectives with respect to the extreme worst points of the same objective. The dashes in Table 3 show that the chosen intermediate trade-off point has no relevance to the selected objective function and hence no improvement or deterioration can be measured.Table 3
Design information of intermediate trade-off points.
Trade-off point
Percentage improvement (positive) or deterioration (negative) of trade-off point in
In comparison with
z
¨
c
z
˙
1
z
˙
2
z
s
-
z
1
z
s
-
z
2
E
1
−77.08
24.65
—
—
—
A
1
48.5
−0.44
—
—
—
B
1
F
1
−0.17
4.12
—
—
—
C
1
44.34
−4.39
—
—
—
D
1
E
2
−17.24
—
1.53
—
—
A
2
65.82
—
−1.23
—
—
B
2
F
2
−2.66
—
7.58
—
—
C
2
6.59
—
−0.18
—
—
D
2
E
3
−12.37
—
—
9.28
—
A
3
5.73
—
—
−6.23
—
B
3
F
3
−15.48
—
—
10.49
—
C
3
3.25
—
—
−4.61
—
D
3
E
4
−20.27
—
—
—
39.81
A
4
0.26
—
—
—
−0.39
B
4
F
4
−15.53
—
—
—
31.99
C
4
0.77
—
—
—
−0.8
D
4Furthermore, as can be observed in Figures5, 6, and 7, MOHTS also yields better nondominated set of Pareto optimal solutions with better and wider spread of the optimal solutions than NSGA-II in the rest of the three 2-objective optimization studies.The comparison of the extreme trade-off points given by the proposed algorithm MOHTS, NSGA-II, MUGA [21], and a combined PSO-GA based MOEA [22] is presented in Tables 4–7. The points Ai′ and Bi′ are the extreme left and bottom points given by MUGA, and points Ai′′ and Bi′′ are the extreme left and bottom points given by combined PSO-GA based MOEA. As can be observed from the last two columns of the tables, MOHTS either outperforms or performs similar to NSGA-II but clearly outperforms the previous work by Bouazara [23], Jamali et al. [21], and Mahmoodabadi et al. [22]. Since Bouazara [23] decomposed the multiobjective optimization into single objective optimization by using the weighted coefficients, his study yielded a single optimum solution.Table 4
The comparison of the extreme trade-off points (as shown in Figure3) in terms of the values of objective functions and their corresponding design variables between present work (MOHTS), NSGA-II, and previous work by Bouazara [23] Jamali et al. [21], and Mahmoodabadi et al. [22].
Trade-off point
k
s
s (N/m)
k
s
1 (N/m)
k
s
2 (N/m)
c
s
s (Ns/m)
c
s
1 (Ns/m)
c
s
2 (Ns/m)
r (m)
z
¨
c (m/s2)
z
˙
1 (m/s)
A
1  (MOHTS)
77448.364
10000
19940.983
3991.464
500
1997.748
0.5
1.831
0.495
B
1  (MOHTS)
61311.317
20000
20000
4000
2000
1904.71
0.047
4.06
0.358
C
1  (NSGA-II)
146432.050
10000
10000
1369.147
1739.995
1798.133
0.415
2.129
0.388
D
1  (NSGA-II)
81907.529
19999.348
17542.899
3517.244
2000
1767.371
0.496
3.832
0.358
A
1
′  from [21]
146863
10039.2
10313.7
3941.18
1917.65
1958.82
0.5
2.807
0.417
B
1
′  from [21]
57451
10039.2
18196.1
2517.65
964.706
1476.47
0.3412
3.297
0.407
A
1
′
′  from [22]
146728.6
10006.86
10020.65
3964.849
1806.955
1947.486
0.498
2.794
0.415
B
1
′
′  from [22]
149700.2
10003.23
19907.76
1796.857
940.0627
548.968
0.499
3.331
0.405
P  from [24]
105630
14958
14985
2548
1956
1967
0.288
3.089
0.375Table 5
The comparison of the extreme trade-off points (as shown in Figure4) in terms of the values of objective functions and their corresponding design variables between present work (MOHTS), NSGA-II, and previous work by Bouazara [23], Jamali et al. [21], and Mahmoodabadi et al. [22].
Trade-off point
k
s
s (N/m)
k
s
1 (N/m)
k
s
2 (N/m)
c
s
s (Ns/m)
c
s
1 (Ns/m)
c
s
2 (Ns/m)
r (m)
z
¨
c (m/s2)
z
˙
2 (m/s)
A
2  (MOHTS)
76695.747
10000
19925.858
4000
500
1994.622
0.5
1.831
0.361
B
2  (MOHTS)
149997.760
18093.905
19998.477
1586.362
500
1998.436
0
6.282
0.351
C
2  (NSGA-II)
146158.890
10000.508
10000
2283.425
1776.585
1905.497
0.415
2.102
0.386
D
2  (NSGA-II)
122584.350
10159.024
19999.999
3893.357
1988.289
1999.973
0.270
2.311
0.357
A
2
′  from [21]
146863
10039.2
10313.7
3941.18
1917.65
1958.82
0.5
2.807
0.434
B
2
′  from [21]
52745.1
10000
10000
3635.249
670.588
1658.82
0.058
3.475
0.430
A
2
′
′  from [22]
131574.7
10027.98
10107.79
3292.868
1906.304
1935.256
0.478
2.818
0.433
B
2
′
′  from [22]
50504.37
11852.13
10000.4
3737.814
502.1251
1112.495
0.003
5.1074
0.428
P  from [24]
105630
14958
14985
2548
1956
1967
0.288
3.089
0.377Table 6
The comparison of the extreme trade-off points (as shown in Figure5) in terms of the values of objective functions and their corresponding design variables between present work (MOHTS), NSGA-II, and previous work by Bouazara [23], Jamali et al. [21], and Mahmoodabadi et al. [22].
Trade-off point
k
s
s (N/m)
k
s
1 (N/m)
k
s
2 (N/m)
c
s
s (Ns/m)
c
s
1 (Ns/m)
c
s
2 (Ns/m)
r (m)
z
¨
c (m/s2)
Z
S
-
Z
1 (m/s)
A
3  (MOHTS)
77446.455
10000
19982.171
3996.929
500.853
1998.217
0.5
1.83
0.079
B
3  (MOHTS)
92412.553
10000
20000
4000
2000
2000
0.5
2.182
0.068
C
3  (NSGA-II)
77397.693
10005.791
18820.197
3997.672
500
1914.903
0.5
1.867
0.08
D
3  (NSGA-II)
107565.51
10147.906
19940.144
3604.317
2000
1999.218
0.466
2.229
0.068
A
3
′  from [21]
146863
10039.2
10313.7
3941.18
1917.65
1958.82
0.5
2.807
0.09
B
3
′  from [21]
57058.8
10000
10000
3588.24
2000
1676.47
0.114
3.045
0.08
A
3
′
′  from [22]
147769.6
10012.38
10037.35
3926.99
1997.283
1879.508
0.498
2.796
0.0896
B
3
′
′  from [22]
50006.19
10000.09
10000.09
3999.446
1999.998
785.631
4.28
E
-
05
3.605
0.0866
P  from [24]
105630
14958
14985
2548
1956
1967
0.288
3.089
0.0844Table 7
The comparison of the extreme trade-off points (as shown in Figure6) in terms of the values of objective functions and their corresponding design variables between present work (MOHTS), NSGA-II, and previous work by Bouazara [23], Jamali et al. [21], and Mahmoodabadi et al. [22].
Trade-off point
k
s
s (N/m)
k
s
1 (N/m)
k
s
2 (N/m)
c
s
s (Ns/m)
c
s
1 (Ns/m)
c
s
2 (Ns/m)
r (m)
z
¨
c (m/s2)
Z
S
-
Z
2 (m/s)
A
4  (MOHTS)
78255.336
10000
19984.602
3996.583
502.137
1998.292
0.5
1.832
0.056
B
4  (MOHTS)
104827.01
10000
20000
4000
1999.516
2000
0.5
2.209
0.034
C
4  (NSGA-II)
92998.785
10002.354
19987.524
3995.824
601.103
1780.096
0.478
1.945
0.051
D
4  (NSGA-II)
114107.32
10354.485
19910.39
3637.269
2000
1983.220
0.485
2.264
0.035
A
4
′  from [21]
146863
10039.2
10313.7
3941.18
1917.65
1958.82
0.5
2.807
0.06
B
4
′  from [21]
92745.1
10000
18313.7
2941.18
576.471
1994.12
0.484
3.551
0.052
A
4
′
′  from [22]
147107.3
10001.45
10020.01
3941.79
1843.581
1999.742
0.499
2.791
0.059
B
4
′
′  from [22]
149232.8
10001.75
19998.5
3959.178
500.222
1999.908
0.499
3.655
0.05
P  from [24]
105630
14958
14985
2548
1956
1967
0.288
3.089
0.054The response of the seat acceleration, front suspension deflection, and rear suspension deflection of the Pareto optimal design points obtained from MOHTS and NSGA-II and that of the optimal point found by Bouazara [23] are as shown in Figures 8–10, respectively. Figure 8 shows the seat acceleration response of the half car ride model for the design points A3 of MOHTS, C3 of NSGA-II, and P of the work by Bouazara [23]. Figure 9 depicts the front suspension deflection response of the system for the design points which yield the minimum front suspension deflection, and these points are B3 of MOHTS, D3 of NSGA-II, and P of the work by Bouazara [23]. Similarly rear suspension deflection response of the system for the optimal points B4 of MOHTS, D4 of NSGA-II, and P of the work by Bouazara [23] is demonstrated in Figure 10.Figure 8
Seat acceleration response of the half car model for different optimal design points.Figure 9
Front suspension deflection response of the half car model for different optimal design points.Figure 10
Rear suspension deflection response of the half car model for different optimal design points.Lastly, 5-objective optimization study, which seeks the Pareto front for simultaneous optimization of all five objectives of (17), with the same algorithm parameters of 2-objective optimization study results, is investigated. Since 5-objective optimization study deals with the five-dimensional objective space and hence for better understanding of the obtained Pareto front, we use the idea proposed by [30] of using biobjective cross sections (slices) of the Pareto front. Such slices are known as decision maps and they give a clear picture of trade-offs between all the objectives. Such decision maps help the decision maker to better understand the obtained Pareto front. However this approach does not take into account the instability of the Pareto front. The same four pairs of 2-objective optimization as presented in study 1 are presented in study 2 as well with their original Pareto front being superimposed over the Pareto front of 5-objective optimization study. Figures 11–14 are the various planes of the Pareto front given by MOHTS and Figures 15–18 are the same planes of the Pareto front given by NSGA-II. In comparison with the 2-objective optimization study, 5-objective optimization study offers more choices for the decision maker. Figure 11 demonstrates the plane of z¨c,z˙1 from the 5-objective optimization study obtained by the proposed MOHTS algorithm and as can be seen Pareto front is widely spread as compared to 2-objective optimization Pareto front for the same two objectives. This 5-objective Pareto front of Figure 11 has also been depicted in Figures 12, 13, and 14 with a view from the planes z¨c,z˙2, z¨c,zs-z1, and z¨c,zs-z2, respectively. Similarly Figures 15–18 show the Pareto fronts of 5-objective optimization study with their corresponding 2-objective Pareto front superimposed in the planes of z¨c,z˙1, z¨c,z˙2, z¨c,zs-z1, and z¨c,zs-z2, respectively. In summary, the comparison of the proposed MOHTS algorithm with NSGA-II indicates that MOHTS has produced a widespread and better Pareto front than NSGA-II. We ran both MOHTS and NSGA-II for 10 runs and computed the average time taken to solve the optimization problem in study 1. MOHTS and NSGA-II took 1567 and 1364 seconds on average, respectively, showing NSGA-II edges over MOHTS in computational time required to solve the numerical problem. Apart from the Pareto front and computational time, parameter or performance sensitivity could be selected as a performance comparison criterion between MOHTS and NSGA-II. However tuning the algorithm parameters to obtain the best possible outcomes requires exhaustive studies and thus to be fair in making the comparison we choose the same mutation and crossover probabilities for both the algorithms as indicated earlier in this section.Figure 11
Pareto front of seat acceleration and front tire velocity in 5-objective optimization and 2-objective optimization using MOHTS for numerical study 1.Figure 12
Pareto front of seat acceleration and rear tire velocity in 5-objective optimization and 2-objective optimization using MOHTS for numerical study 1.Figure 13
Pareto front of seat acceleration and front suspension deflection in 5-objective optimization and 2-objective optimization using MOHTS for numerical study 1.Figure 14
Pareto front of seat acceleration and rear suspension deflection in 5-objective optimization and 2-objective optimization using MOHTS for numerical study 1.Figure 15
Pareto front of seat acceleration and front tire velocity in 5-objective optimization and 2-objective optimization using NSGA-II for numerical study 1.Figure 16
Pareto front of seat acceleration and front tire velocity in 5-objective optimization and 2-objective optimization using NSGA-II for numerical study 1.Figure 17
Pareto front of seat acceleration and front tire velocity in 5-objective optimization and 2-objective optimization using NSGA-II for numerical study 1.Figure 18
Pareto front of seat acceleration and front tire velocity in 5-objective optimization and 2-objective optimization using NSGA-II for numerical study 1.


## 4.2. Numerical Study 2

In contrast to the first numerical study, we choose more realistic approach in developing the objective functions for optimizing the design of a half car model which is considered with the same set of design variables, vehicle constants, and bounds on the design variables as chosen in study 1. As stated earlier, vehicle suspension system’s three opposing performance indices or objectives are used to evaluate the suspension system’s dynamic response. A better method to estimate the ride discomfort is to measure the root mean square (RMS) of the vertical acceleration felt by the passenger/driver. In this study RMS of the seat acceleration is chosen as the first objective. Equivalently road holding is determined through the RMS of tire deflection, and suspension rattle space is evaluated by measuring the RMS of suspension deflection. The RMS of front and rear suspension deflection and RMS of front and rear tire deflections are selected as other four objectives to form the 5-objective optimization problem in this study. The optimization problem is stated as follows.Find the design vector,X=ksscssks1cs1ks2cs2rT, which minimizes(18)FX=f1Xf2Xf3Xf4Xf5XTf1X=1T∫0Tz¨c2dtf2X=1T∫0Tzs-z12dt,f3X=1T∫0Tzs-z22dtf4X=1T∫0Tz1-zp12dt,f5X=1T∫0Tzs-zp22dtwith no constraints.Unlike the previous study, the Pareto fronts obtained from 2-objective studies are superimposed over the corresponding Pareto fronts yielded from the 5-objective studies to shorten the length of the paper. We again use the same slicing method for visualizing the Pareto front in two-dimensional space to better demonstrate the effectiveness of MOHTS and NSGA-II. Figures19–22 show the planes (slices) of f1,f2, f1,f3, f1,f4, and f1,f5, respectively, for the 5-objective optimization study using MOHTS and likewise Figures 23–26 depict the same planes in the same order for the 5-objective optimization using NSGA-II. The objective function and design variable values of the extreme points of the 2-objective studies’ Pareto fronts as shown in Figures 19–26 are given in Table 8. It is clear from the table that MOHTS is superior to NSGA-II in minimizing the objectives of 2-objective optimization studies and the minimum values of each objective are shown in bold letters.Table 8
The comparison of the extreme trade-off points (as shown in Figures18–25) in terms of the values of objective functions and their corresponding design variables between present work (MOHTS) and NSGA-II.
Trade-off point
k
s
s (N/m)
k
s
1 (N/m)
k
s
2 (N/m)
c
s
s (Ns/m)
c
s
1 (Ns/m)
c
s
2 (Ns/m)
r (m)
f1(m/s2)
f
2 (m)
f
3 (m)
f
4 (m)
f
5 (m)
A
5  (MOHTS)
62325.038
10000.000
20000.000
4000.000
500.000
1999.952
0.500
0.850
0.035
0.026
0.003
0.003
B
5  (MOHTS)
105868.150
10000.000
20000.000
4000.000
2000.000
2000.000
0.500
1.062
0.033
0.016
0.003
0.002
C
5  (NSGA-II)
111650.910
10004.447
13795.435
3676.384
1411.603
1970.654
0.434
1.017
0.037
0.022
0.003
0.002
D
5  (NSGA-II)
113125.130
10000.024
13996.260
3668.958
1477.650
1997.922
0.454
1.021
0.036
0.021
0.003
0.002
A
6  (MOHTS)
66552.042
10000.000
20000.000
4000.000
500.000
1967.310
0.500
0.850
0.035
0.025
0.003
0.003
B
6  (MOHTS)
104962.350
10000.000
20000.000
4000.000
2000.000
1979.535
0.500
1.063
0.033
0.016
0.003
0.002
C
6  (NSGA-II)
95623.254
10030.399
15187.577
1680.633
1478.532
1995.701
0.500
1.045
0.036
0.021
0.003
0.002
D
6  (NSGA-II)
132369.180
10409.407
19999.716
2728.990
2000.000
1730.412
0.450
1.138
0.034
0.016
0.003
0.003
A
7  (MOHTS)
63426.812
10000.002
19976.225
3999.546
500.357
1962.724
0.500
0.849
0.035
0.026
0.003
0.003
B
7  (MOHTS)
62028.590
10000.000
20000.000
4000.000
500.000
500.000
0.500
1.584
0.048
0.028
0.002
0.009
C
7  (NSGA-II)
57838.283
10434.762
19752.225
2306.963
736.508
1672.299
0.500
1.019
0.038
0.026
0.003
0.003
D
7  (NSGA-II)
58174.332
10430.336
20000.000
2334.326
838.576
500.000
0.500
1.463
0.045
0.025
0.002
0.007
A
8  (MOHTS)
62635.017
10000.000
20000.000
4000.000
500.000
2000.000
0.500
0.850
0.035
0.026
0.003
0.003
B
8  (MOHTS)
50000.000
20000.000
10000.000
4000.000
2000.000
2000.000
0.486
1.908
0.040
0.026
0.007
0.001
C
8  (NSGA-II)
65778.794
10000.000
18885.969
3710.928
500.052
1932.241
0.500
0.874
0.036
0.026
0.003
0.003
D
8  (NSGA-II)
50308.460
19948.963
10000.000
3521.482
1971.198
1972.075
0.486
1.900
0.040
0.026
0.007
0.001Figure 19
Pareto front of seat acceleration RMS and front suspension deflection RMS in 5-objective optimization and 2-objective optimization using MOHTS for numerical study 2.Figure 20
Pareto front of seat acceleration RMS and rear suspension deflection RMS in 5-objective optimization and 2-objective optimization using MOHTS for numerical study 2.Figure 21
Pareto front of seat acceleration RMS and front tire deflection RMS in 5-objective optimization and 2-objective optimization using MOHTS for numerical study 2.Figure 22
Pareto front of seat acceleration RMS and rear tire deflection RMS in 5-objective optimization and 2-objective optimization using MOHTS for numerical study 2.Figure 23
Pareto front of seat acceleration RMS and front suspension deflection RMS in 5-objective optimization and 2-objective optimization using NSGA-II for numerical study 2.Figure 24
Pareto front of seat acceleration RMS and rear suspension deflection RMS in 5-objective optimization and 2-objective optimization using NSGA-II for numerical study 2.Figure 25
Pareto front of seat acceleration RMS and front tire deflection RMS in 5-objective optimization and 2-objective optimization using NSGA-II for numerical study 2.Figure 26
Pareto front of seat acceleration RMS and rear suspension deflection RMS in 5-objective optimization and 2-objective optimization using NSGA-II for numerical study 2.As can be observed from Figures19 and 23, the Pareto front for the 5-objective optimization is quite similar for both MOHTS and NSGA-II and this is also true for the rest all the planes, thus, showing that the potential of the proposed MOHTS algorithm is as good as NSGA-II in finding the Pareto front solutions of the multiobjective optimization problem. Moreover, when 2-objective optimization Pareto fronts from MOHTS are compared for all the shown pairs of objectives with NSGA-II, it appears that the proposed MOHTS algorithm dominates NSGA-II in obtaining the better nondominated Pareto front with the widespread (diveresed) set of optimal solutions. The reason behind the success of MOHTS in yielding a better Pareto front can be attributed to the search mechanism of HTS, as it takes a good care in maintaining the diversified set of solutions in the search space with its better exploration capabilities, which in turn allows the nondominated sorting and crowding distance approach mechanisms to produce a widespread Pareto front. Furthermore, as exploitation capability has also been wisely implemented in HTS algorithm it helps the dominance sorting and diversity preserving mechanisms in producing a better nondominated set. In contrast to the first numerical study, the average computational time for MOHTS and NSGA-II in this study is 1545 and 1633 seconds, respectively. The reason behind this conflict can be attributed to the uncertainties in the numerical solution of the half car model’s dynamic equations of motion.


## 5. Conclusions

In this paper, a new multiobjective optimization algorithm called multiobjective heat transfer search (MOHTS) has been proposed which works on the search technic of recently developed metaheuristic HTS and on the elitist nondominated sorting for obtaining different nondomination levels and crowding distance approach to preserve the diversity among the optimal set of solutions of NSGA-II. The proposed algorithm has been applied to the MOOP of the vehicle dynamics to optimize the passive suspension parameters of the vehicle by employing half car lumped mass-spring-damper model. Since the chosen MOOP represents a complex optimization problem because of five conflicting objectives and ordinary differential governing equations of the system, it is assumed that it can test the real potential of the proposed MOHTS algorithm. Two optimization studies are presented with a different set of optimization objectives and the nature of the objective functions is higher order linear ordinary differential equations in both studies. Since the proposed multiobjective optimization algorithm (MOHTS) has been developed on the platform of NSGA-II for sorting the population based on the nondominancy and preserving diversity among the obtained set of optimal solutions, its performance at all the stages is compared with NSGA-II. The optimization studies reveal that in 2-objective optimization MOHTS is able to yield the better nondominated Pareto front with the widespread (diveresed) set of optimal solutions as compared to NSGA-II. Furthermore from the analysis of extreme trade-off points obtained from the 2-objective optimization numerical study 1 it was observed that MOHTS produces a better Pareto front of the optimal solutions than multiobjective uniform diversity genetic algorithm (MUGA) and combined PSO-GA based MOEA. In 5-objective optimization studies MOHTS has produced a better and equivalent Pareto front as compared to NSGA-II in numerical studies 1 and 2, respectively. Hence, the performance of the proposed multiobjective optimization algorithm MOHTS was tested on the optimization problem with the linear differential equations of higher order and was found to give better Pareto front as compared to NSGA-II, MUGA, and combined PSO-GA based MOEA.



---

*Source: 2034907-2017-05-03.xml*