Multi-objective design optimization of an engine accessory drive system with a robustness analysis

Journal Pre-proof

Multi-objective design optimization of an engine accessory drive system with a robustness analysis H. Zhu , Y. Hu , W.D. Zhu , W. Fan , B.W. Zhou PII: DOI: Reference:

S0307-904X(19)30546-3 https://doi.org/10.1016/j.apm.2019.09.016 APM 13020

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

8 April 2019 13 August 2019 3 September 2019

Please cite this article as: H. Zhu , Y. Hu , W.D. Zhu , W. Fan , B.W. Zhou , Multi-objective design optimization of an engine accessory drive system with a robustness analysis, Applied Mathematical Modelling (2019), doi: https://doi.org/10.1016/j.apm.2019.09.016

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Inc.

Highlights 

A hybrid mutation particle swarm optimization (HMPSO) algorithm is presented.



Non-parametric statistical tests on ten benchmark functions are conducted.



Key structure parameters of an engine accessory drive system are optimized via the HMPSO for vibration control.



Optimization on an example engine accessory drive system show that the system vibration is greatly improved after optimization.



Suggestions on design of the engine accessory drive system are given according to a robustness analysis.

Corresponding author. Tel.:+86 17394999829. E-mail addresses: [email protected] (H. Zhu)

Multi-objective design optimization of an engine accessory drive system with a robustness analysis

H. Zhu1 , Y. Hu2 , W.D. Zhu3 , W. Fan1 , B.W. Zhou2 1. Department of Mechanics and Engineering Science, Sichuan University, Chengdu 610065, China 2. State Key Laboratory of Mechanical Transmissions , Chongqing University, Chongqing 400044, China 3. Department of Mechanical Engineering, University of Maryland, Baltimore County, Baltim ore, MD 21250, USA

Abstract Designing a good engine accessory drive system becomes a hard work with its increasingly complicated configuration and high demands on its dynamic characteristics. In this work, a hybrid 2

mutation particle swarm optimization (HMPSO) algorithm is presented to optimize the key structure parameters of an engine accessory drive system for its vibration control. The superiority of the HMPSO algorithm against several other concerned metaheuristic algorithms in terms of solution quality and stability are verified by non-parametric statistical tests on ten benchmark functions. The design problem of the engine accessory drive system is a multi-objective optimization problem; the weighted sum method and main target method are applied to convert it to a single-objective one. Optimization on an example engine accessory drive system using the HMPSO algorithm demonstrates obvious improvement in system vibration after optimization. A robustness analysis is conducted to identify the robustness of dynamic responses of the engine accessory drive system with respect to small variations of the design variables relative to the optimal design in the design space, and suggestions on design of an engine accessory drive system are given according to it.

Keywords: optimal design; engine accessory drive system; vibration control; hybrid mutation particle swarm optimization; robustness analysis

1. Introduction An engine accessory drive system is extensively used in vehicles to power engine accessories to provide enough pressures, torques, and electricity for accessory components such as the water pump, 3

alternator, air conditioner, and power steering pump, and so on [1]. Torsional vibrations of accessory pulleys and transverse vibrations of belt spans in an engine accessory drive system have greatly been concerned by engineers [2-6]. Large vibrations may appear when the engine accessory drive system is poorly designed. Traditional design of an engine accessory drive system largely relies on experience and experiments. Several original design plans need to be firstly given by designers from experience and then a good one will be picked out from them by experimental validation. This design method is expensive and time consuming, and a globally optimal design cannot be strictly guaranteed by this method. Hence, it is of great importance to put forward an effective approach to guide engineers to design an engine accessory drive system and avoid undesirable vibration problems. Design optimization has been used in solving different engineering vibration problems like vibration isolation problem of active suspension systems [7], eigenfrequency separation problems of two- material structures [8] and vibration control of the floating wind turbine [9], and so on. However, only a few researchers use design optimization methods to reduce vibrations of engine accessory drive systems [6,10,11,12]. A sequential quadratic programming and Kuhn-Tucker method are used by Nouri and Zu [6] to optimize an engine accessory drive system for its torsional vibration reduction. Zhu et al. [10] used a sequence quadratic programming method with oval pulley parameters selected as design variables to minimize the overall torsional vibration amplitude of a timing belt drive system with an oval cogged pulley. Hou et al. [11] presented an approach by incorporating sensitivity information into a classic coordinate alternating procedure to find an optimum design. These methods mentioned above are deterministic and gradient-based methods. Hence, a global optimum is hard to be directly obtained by these methods for such a system with 4

many local optima, and original values of design parameters usually need to be changed several times until a most acceptable optimal solution is obtained. Recently, Zhu et al. [12] presented a method to simultaneously modify natural frequencies of an engine accessory drive system via configuration modification to match natural frequencies of the system with its working conditions so that the system resonance conditions are successfully avoided. This is a rough design method and the globally optimal design cannot be guaranteed as well. Numerous global optimization algorithms, such as the genetic algorithm (GA) [13,14], particle swarm optimization (PSO) algorithm [15], firefly algorithm (FA) [16], simulated annealing (SA) algorithm [17], teaching- learning-based optimization (TLBO) [18], and mine blast algorithm (MBA) [19], have been proposed to overcome these drawbacks. As a global optimization algorithm, the PSO algorithm is efficient in computation, easy for implementation, and reliable to search for global optima of various engineering design problems [20-23]. Unlike traditional deterministic optimization methods such as the steepest descend method, quasi-Newton method, and interior-reflective Newton method, differentiable conditions of objective functions are not required to be met for the PSO algorithm. It is capable of solving an optimization problem with strong nonlinearity, high dimensions, and various constraints [24]. However, the standard PSO algorithm is easy to be premature and stick to local optima. To increase the diversity of the particle swarm and/or increase convergence rates, many hybrid optimization algorithms with some operators or other algorithms incorporated into the PSO have been presented [25-30]. Novitasari et al. [25] and He and Wang [26] presented a hybrid algorithm that combines the SA with PSO algorithm to solve constrained optimization problems and optimize a Support Vector Machine, respectively. Wang and Yin [27] introduced a ranking selection strategy 5

into the basic PSO to self-adaptively control optimal solution search behaviors of the particle swarm that generates a new algorithm named the ranking selection-based PSO (RSPSO). The crossover operators and/or mutation operators used in genetic algorithms (GAs) have been used in many studies by putting them into the PSO or combining them with the PSO to yield new algorithms like the modified PSO (MPSO) [28], quantum-behaved PSO using mutation operators with a Gaussian distribution (G-QPSO) [29], and straightforward PSO (SPSO) with a logistic chaotic mutation operator [30], etc. These operators improve swarm diversity and prevent premature convergence and stagnation of the PSO algorithms. Garg [31] recently used the parameter free penalty function to solve various constrained optimization problems. The hybrid optimization algorithms discussed above have been used to solve various specific engineering optimization problems. Motivated by the advantages of the hybrid optimization algorithms mentioned above and to overcome the shortcomings of the present design methods, a hybrid mutation PSO (HMPSO) algorithm is presented here to solve the design problem of an engine accessory drive system for its vibration control. It is a multi-objective optimization problem according to the requirements of the design problem. Hence, the weighted sum method (WSM) and main target method (MTM) are first used to convert the multi-objective optimization problem to a single-objective one and two different objective functions are obtained as a result. A modified inertia factor is adopted to balance the convergence rate and global optima search ability, and a hybrid mutation strategy with six mutation operators are introduced to diversify the swarm in the HMPSO algorithm. Ten benchmark functions are used to test its performance in terms of the optimum solution quality and stability, and the non-parametric statistical tests are conducted to test the statistical significance between the HMPSO algorithm and several other existing metaheuristic algorithms. The HMPSO algorithm and several 6

other existing metaheuristic optimization algorithms are then used to solve the design problem of an example engine accessory drive system. The optimal results from all the algorithms are discussed and compared. A robustness analysis is finally conducted to investigate robustness of dynamic responses of the optimal configuration with respect to variations of design variables in the design space, and suggestions on design of an engine accessory drive system are given according to it. The method presented in this paper can provide guidelines for engineers in designing a good engine accessory drive system.

2. System model Fig. 1 shows a schematic of a prototypical engine accessory drive system. As seen in Fig.1, it contains a tensioner, a serpentine belt, and n pulleys connected to their corresponding accessories, with n=8 here. The tensioner is made up of a tensioner arm with a preloaded torsional spring assembled inside and a tensioner pulley pinned to it. The tensioner arm is pivoted to the frame and the tensioner pulley is in contact with its two adjacent belt spans. Zhu et al. [12, 32] built a hybrid dynamic model for this system. A Cartesian coordinate system (X, Y) and local coordinate system

 xˆ , yˆ  j

j

are established as references to describe geometric positions and motions of system

components. Pulleys whose radii are denoted by Rˆ j  j  1, 2,, n  are numbered by turn in a counterclockwise direction as shown in Fig. 1 and the tensioner pulley is numbered as the i-th pulley where 1
 j  j  1, 2, , n  , tensioner arm rotation t and belt transverse vibrations wˆ j  xˆ , tˆ  . Each pulley rotates with a moment of inertia Jˆ j and viscous damping coefficient Dˆ j within its pulley bearings, 7

and subjects to an accessory torque Mˆ j . The tensioner arm rotates with a moment of inertia Jˆarm , viscous damping coefficients Dˆ t within its pivot bearings, torsional stiffness of the torsional spring

Kˆ r , and an arm length lˆt . Each belt span is modeled as a uniformly moving Bernoulli- Euler beam with an elastic modulus E, a cross-section area A, a bending stiffness EI and a constant moving speed cˆ . Each belt span is pretended with a tension Pˆj at equilibrium and a dynamic tension Pˆd , j during operation. Equations of transverse vibrations of each belt span are derived from the Hamilton’s principle [32,33]:

ˆ

 2 wˆ j  2 wˆ j  2 wˆ j  4 wˆ j 2 ˆ ˆ ˆ ˆ ˆ ˆ  2 c  Pj   c  Pd , j  EI  0, j  1, 2, , n tˆ 2 xˆtˆ xˆ 2 xˆ 4





(1)

Based on a pseudo-static stretch assumption, the dynamic tension Pˆd , j of each belt span is [4] ˆ EA  ˆ ˆ   1 l j wˆ 2 dx  , j  1, 2,..., n; j  i  1, i Pˆd , j  R   R j 1 j 1 j,x j   j j 2 0 lˆj   ˆ EA  ˆ ˆ   l  sin   1 li1 wˆ 2 dx  Pˆd ,i 1  R   R i i t t 1 i 1, x i 1   i 1 i 1 2 0 lˆi 1   ˆ EA  ˆ ˆ   l  sin   1 li wˆ 2 dx  Pˆd ,i  R   R i 1 i 1 t t 2 i,x i   i i 2 0 lˆi  

(2)

(3)

(4)

The dynamic tension in each belt span due to its transverse vibration is usually small and can be neglected so that the last terms in Eqs. (2)-(4) are omitted. Boundary conditions of all belt spans are wˆ j |lˆ  0, j  1, 2, , n, j  i  1, i j

wˆ i 1 |0  wˆ i |lˆ  0, wˆ i 1 |lˆ  lˆtt cos 1 , wˆ i |0  lˆt t cos 2 i 1

i

EI

 2 wˆ j  0, tˆ  xˆ 2j

 0, EI

   0, j  1, 2, , n

 2 wˆ j lˆj , tˆ xˆ 2j

(5) (6) (7)

Equations of pulleys and tensioner arms are [32]





Jˆ j j  D j j  Rˆ j Pˆd , j  Pˆd , j 1  Mˆ j , j  2,3,, n

8

(8)

 3 wˆ i 1 2 2 ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ   J tt  Dtt   kt lt  ki 1lt sin 1  ki lt sin  2  ltt  EI | ˆ lˆt cos 1 xˆ 3 xˆ li1  3 wˆ  EI 3i |xˆ 0 lˆt cos  2  kˆi 1lˆt Rˆi 1ˆi 1  Rˆiˆi sin 1  kˆilˆt Rˆiˆi  Rˆi 1ˆi 1 sin  2 xˆ wˆ i 1 wˆ     Pˆi 1  Pˆd ,i 1  ˆ c 2 |xˆ lˆ  ˆ c i 1 |xˆ lˆ  lˆt cos 1 i  1 i  1 xˆ xˆ   wˆ i wˆ     ˆ c 2  Pˆd ,i  Pˆi |xˆ 0  ˆ c i |xˆ 0  lˆt cos  2  0 ˆ xˆ t  















(9)



where Jˆt  Jˆarm  mlˆt is the moment of inertia of the tensioner, in which m is the mass of the tensioner pulley, kˆ j  EA / lˆj and kˆt  Kˆ r / lˆt 2 . By taking the Galerkin procedure using an accurate spatial discretization method [32], the system model is discretized and the partial differential equations are transformed to ordinary differential equations (ODEs). Torsional vibrations of the accessory pulleys and tensioner arm, and transverse vibrations of the belt spans can then be calculated by solving the resulting ODEs using an ODE solver. The accuracy of the dynamic model of the engine accessory drive system has already been verified in our previous work in Ref. [32], and thus it will not be discussed here again.

3. Design problem 3.1. Introduction of the design problem Identifying a most suitable structure for an engine accessory drive system is a challenging design problem with its increasingly complicated configuration and high demands on its dynamic characteristics. A bad design in practical working conditions can cause serious vibration problems, which can result in some undesirable effects such as large dynamic tension fluctuations due to pulley torsional vibrations and belt vibrations, shortened belt life due to periodical stress variations, belt slip, and noise problems. The design task of a best engine accessory drive syste m is to simultaneously minimize 9

transverse vibration amplitudes of all belt spans and torsional vibration amplitudes of all accessory pulleys and the tensioner arm in various operational conditions. In addition, the belt tension should be kept within a range to ensure the conveyer ability of the belt, and avoid belt slip on pulleys and early belt fatigue. The minimum belt tension is assured by a pretension of the belt provided by a pre-torqued tensioner.

3.2. Objective function Dynamic characteristics of an engine accessory drive system are greatly affected by geometric structures of accessory pulleys and their installation conditions, and the specification of the tensioner. Generally speaking, geometric structures and assembling positions of accessories and their corresponding pulleys are given and it is not convenient to change them. Tensioner structure can be flexibly adjusted and its structure optimization is the most effective way to improve the overall dynamic characteristics of an engine accessory drive system [12]. Therefore, design variables are selected from structure parameters of the tensioner including the initial installation angle of the tensioner arm ˆt 0 , torsional stiffness of the tensioner torsional spring Kˆ r , tensioner arm length lˆt , and initial torque of the tensioner Tˆt 0 . An objective function that should meet all design targets needs to be established to assess the behavior of each selection of design variables when performing an optimization method [34]. The optimization objective demands minimizing transverse vibration amplitude of each belt span and torsional vibration amplitudes of each pulley and the tensioner arm at the same time. Obviously, it is a multi-objective optimization problem here. The design sub-objectives are conflicting as minimum transverse vibration amplitudes of all belt spans and minimum angular vibration amplitudes of all accessory pulleys actually cannot be achieved at the same time for a real engine accessory drive 10

system in different operational conditions. The WSM and MTM are adopted here to deal with the multi-objective optimization problem by converting it to a single-objective one. According to the design problem mentioned above, the objective function for optimization based on the WSM is given as below:





1  min 1 ˆt 0 , Kˆ r , lˆt , Tˆt 0   n n  2 w 2 p 2 1 ˆt 0 , Kˆ r , lˆt , Tˆt 0   ˆj  Aj   ˆm  Am   ˆt  At   j 1 m2 



(10)



w where A j , Amp and At are vibration amplitudes of the j-th belt span, m-th pulley and tensioner arm,

respectively; ˆj , ˆm and ˆt are corresponding weighting factors of each term. From Eq. (10), the objective function 1 is aimed at finding a set of design variables that can minimize the weighted quadratic sum of vibration amplitudes of all the components in the engine accessory drive system. Similarly, the objective function for optimization based on the MTM is given as below:





 2  min  2 ˆt 0 , Kˆ r , lˆt , Tˆt 0    2 ˆt 0 , Kˆ r , lˆt , Tˆt 0  max A1w , A2w , ., Anw , A2p , A3p ,, Anp , At 







(11)



where the objective function  2 is aimed to find a set of design variables to minimize the maximum vibration amplitudes of all the components in the engine accessory drive system. Considering the limit of installation space, conveyer ability and reliability of the belt, and to avoid motion interference as well as belt slippage, the design variables are subjected to the following constraints: ˆtl0  ˆt 0  ˆtu0  l  Kˆ r  Kˆ r  Kˆ ru s.t.  l u lˆt  lˆt  lˆt  ˆl ˆ ˆu Tt 0  Tt 0  Tt 0

(12)

where the superscripts u and l stand for upper and lower limits of each design variable, respectively. 11

4. Optimization algorithm 4.1. Basic concepts of the PSO algorithm There are various optimization algorithms that can be roughly divided into two categories: a stochastic algorithm and a deterministic one. The stochastic algorithm such as the PSO algorithm, which is demonstrated to be capable of finding global optima with high efficiency and reliability, is widely adopted to solve various engineering design problems [21,24,35-37]. The PSO algorithm consists of a swarm of particles randomly moving in a design space. The position of each particle in the design space represents a possible solution of a specific design optimization problem. Every particle has a speed and spreads in the design space. The position and speed of each particle are continuously changed by v i    1   v i     c1r1  Pi     xi      c2 r2  Pg     xi    

xi    1  xi     v i    1

(13) (14)

where v i    and xi    represent the speed and position of the i-th particle at the  -th time step, respectively; Pi    represents the historical best position of the i-th particle up to the  -th time step and Pg    represents the global best position of the population by now;  represents an inertia factor; random numbers r1 and r2 are limited in a range

 0, 1 ; two accelerating factors

c1

and c2 are constants that scale influences of Pi    and Pg    , respectively.

4.2. HMPSO algorithm The HMPSO algorithm is presented here by introducing a modified inertia factor in the basic PSO algorithm to balance the convergence rate and space search ability during an iteration process, and six hybrid mutation operators are also introduced to diversify the swarm in the algorithm. 12

4.2.1 Adaptive inertia factor The inertia factor  that has a great influence on the search process is a constant in the basic PSO algorithm. A large inertia factor can improve the global search ability of the algorithm in the design space while a small one can enhance its local search ability. To balance the convergence rate and space searching ability of the PSO algorithm during a searching process, a linearly decreasing dynamic inertia factor or accelerating factor is usually adopted [24,38,39]. Clerc [40] presented a constriction factor K to control the convergence speed. In this work, an modified inertia factor     , which is nonlinearly decreasing during the iteration process, is used to control the search speed:      2tot 

     min  max  min  cos 

(15)

where max and min are the maximum and minimum search velocities, respectively, and tot is the total number of iterations. 4.2.2 Hybrid mutation operators The global optimum search ability of the basic PSO algorithm mainly relies on mutual interaction (social learning) and influence of individual particles (cognitive learning). Particles move towards the currently global best position of the swarm in each iteration. A particle can escape from a local optimum with the help of neighboring particles. But if most of its neighboring particles are limited to a local extreme point, it is attracted to the trap of the local optimum, and as a result, premature convergence of the algorithm and the stagnation phenomenon [41] occur. Mutation operators that are often used in the GA [42,43] are widely adopted to refine the basic PSO algorithm [41,44,45]. Six hybrid mutation operators are introduced in the basic HMPSO algorithm to overcome 13

drawbacks of the basic PSO algorithm. These mutation operators permit the swarm to escape from a local optimum and search in different regions within the solution space (design space). Mutations of particular particles diversify the population and thus premature convergence and the stagnation phenomenon are prevented. Therefore, the presented HMPSO algorithm has the capability to automatically balance its global and local searching capabilities with a good convergence performance. Once a particle satisfies the mutation criterion below, it will be picked out and its position will be adjusted via mutation operators in each iteration: fi     f g    fg    



(16)

where f i    is the current fitness of the particle i in the λ-th iteration, f g    is the best fitness of the swarm corresponding to its global best position,  is an infinitesimal to avoid a zero     is the mutation probability which is linearly  tot 

denominator and  =max  max  min   decreasing with iteration.

To diversify the population of particles, the following six hybrid mutation operators are used to adjust the position of each selected particle: (a) Simple random mutation operator k

where

k

x il and

k

  xi    1  k xil  sin  rand     k xiu  k xil  2 

(17)

x iu are vectors of lower and upper position bounds of the particle i in the k-th

dimension, respectively, and rand is a uniformly distributed random number in [0, 1]. (b) Random mutation operator based on the current particle position k

xi    1  k xi      k xiu  k xil  sin  rand -0.5   14

(18)

(c) Random mutation operator based on the best particle position k

xi    1  k Pi      k xiu  k xil  sin  rand -0.5  

(19)

(d) Random mutation operator based on the global best position k

xi    1  k Pg      k xiu  k xil  sin  rand -0.5  

(20)

(e) Logistic chaotic mutation operator [46] k

xi    1  1  g     k xil  g    

(21)

where   k xil  z    1   k xiu  k xil  is the mutation size, in which z    1 = z    1  z     is a logistic chaotic sequence [46], with  being a controlling factor and z  0  being the initial value of the logistic chaotic sequence; and g    is defined by g     1     1 /   , in which  is 

a coefficient for controlling the shrinkage rate. (f) Mirror reflection mutation operator [47] k

xi    1  k xiu  k xil  k xi   

(22)

Once the condition in Eq. (16) is met, the position of the particle i is updated by mutation. Which one of the above six mutation operators is selected to adjust the current particle position is determined according to the value of a random probability distribution function p      0,1 , as shown in Table 1, where selection probabilities of the six mutation operators satisfy

0  p1  p2  p3  p4  p5  1 . The mutation process is redone if the position of a particle after mutation is located outside the design space.

4.2.3 Implementation procedure of the HMPSO algorithm The main implementation procedure of the HMPSO algorithm presented is shown in Fig. 2 and it is briefly described as below. 15

Step 1: Defining the fitness function and design space according to the design problem. Step 2: Setting initial values of the parameters in the HMPSO algorithm including the population size M, inertia factor  , accelerating factors c1 and c2 , mutation probability  , selection probability of each mutation operator p1, p2 , p3 , p4 , and p5 , and controlling factors of the logistic chaotic mutation operator  and  . Step 3: Swarm initialization: randomly generating a swarm population with a size of M and the initial position of each particle is given by xi  0   xil  rand   xiu  xil 

(23)

Evaluating the fitness function value of each initially generated particle. Ranking the positions of the particles according to the evaluation and identifying the initial best particle position Pi  0  and global best and worst positions Pg  0  and Pw  0  , respectively. Step 4: Updating the current position xi    and speed v i    of the particle i according to Eqs. (13), (14), and (16). Step 5: Evaluating the current fitness of each particle, and updating the best particle position

Pi    and global best and worst positions Pg    and Pw    . Step 6: Updating positions of selected particles according to the mutation criterion to diversify the swarm. If the new particle’s position is located beyond the lower and upper position bounds, the mutation process is repeated until it is located within the range  xli , xui  . Step 7: Evaluating the current fitness of the newly found particle after mutation, and comparing it with the fitness values of the global best particle and worst particle. If the new particle is better than the best particle, it will replace the best particle. Otherwise, the worst particle is replaced by it. Step 8: Repeating steps 4-7 until the ending criterion is met. 16

4.3. Benchmark function test 4.3.1 Benchmark test functions In this section, ten benchmark test functions, which have been extensively used to evaluate performances of new optimization algorithms [47-49], are adopted to assess performance of the HMPSO algorithm. The ten benchmark test functions are listed in Table 2 and the formulations of the last five benchmark test functions are given in the Appendix. The Griewank function is a multimodal one with lots of local minima and the remainder are simple unimodal ones [36,39].

4.3.2 Parameters setting The dynamic inertia factor of the HMPSO algorithm is identified by Eq. (15), where its lower limit min and upper limit max are set to 0.4 and 0.9, respectively. A large inertia factor represents a high search speed. At the beginning of the optimization process, the largest inertia factor is used so that more space is searched to improve global search efficiency. The inertia factor is then gradually reduced to the minimum to facilitate convergence of the algorithm. Both of the two accelerating factors

c1 and c2 are set to 1.65. The mutation probability is

linearly decreasing from 0.7 to 0.15 (i.e., max  0.7 and min  0.15 ) and the selection probabilities of the six mutation operators are p j  j / 6, where j  1, 2,

,5 . Values of the controlling factors

of the logistic chaotic mutation operator are     4 , and the initial value z  0  of the logistic chaotic sequence is randomly generated in the interval

 0, 1

and z  0   0.25, 0.5, 0.75, 1 .

The swarm size M and maximum number of iterations S for the HMPSO algorithm and dimension of all benchmark functions D are set to 20, 10000 and 30, respectively. 4.3.3 Results of HMPSO on benchmark functions Fig. 3 presents convergence histories of the proposed HMPSO algorithm on six benchmark 17

functions. Fig. 4(a) and 4(b) show the behaviors of the particles in the basic PSO and HMPSO, respectively. We can observe that in the HMPSO, the particles are more active and moving permanently to explore the search space compared with the particles in the basic PSO. The stagnation is alleviated by the introduced hybrid mutation operators that introduce diversity into the population. This process allows the swarm to escape from the local optima and to search in more different zones within the design space. As a result, the HMPSO has the automatic balance ability between global and local search abilities and thus it has better convergence than the basic PSO. The proposed HMPSO algorithm is compared with seven other famous metaheuristic algorithms including PSO with time-varying acceleration coefficients (PSO-TVAC) [49], global PSO (GPSO) [50], KPSO [40,41], dynamic multi-swarm PSO (DMS-PSO) [51], evolution strategy with covariance matrix adaptation (CMA-ES) [52], the newly proposed straightforward PSO (SPSO) [36] and MBA [19], in solution quality and stability. Parameters of these compared algorithms are set according to their references. The average and standard deviation (SD) values of the optimal solutions gained over 30 independent runs are detailed in Table 3. Performances of these algorithms on one of the ten benchmark test functions are ranked according to their mean values or SD values. Ranks of each algorithm are shown in brackets and the best results are shown in bold in Table 3. According to the ranks, the proposed HMPSO has better performances against other metaheuristic algorithms on the Sphere, Griewark, Ackley, HappyCat, MS and ESF6 functions in terms of both mean optimal solutions and SD values. On the Quadric and Weierstrass functions, the HMPSO is only second to the SPSO and DMS-PSO, and performs better than the other algorithms. On the Quadric Noise function, HMPSO, MBA, GPSO and PSO-TVAC have similarly good performances in terms of both mean and SD values that are better than those of the other algorithms. On the Rosenbrock function, HMPSO, CMS-EA and SPSO have similarly good performances that are better than those of the other algorithms. From the above analyses, the HMPSO has the best comprehensive 18

performance among all these algorithms on the ten benchmark functions according to their average ranks in terms of both the mean and SD values. Non-parametric statistical tests are popularly used to compare performances of two algorithms on multiple common problems by conducting ranking-based transformations [53,54], because no pre-assumptions are required on the probability distributions of the sampled data being assessed. In order to test the statistical significance, the Wilcoxon signed-rank test is used here to non-parametrically detect the significant differences between the statistical results of the HMPSO and those of the remainder algorithms using a 95% confidence interval. Null hypothesis is that there is no significant difference between the mean or SD values of these two algorithms. Table 4 summarizes the Wilcoxon signed-rank test results. R+ and R− in Table 4 represent sums of positive and negative ranks, respectively, and the p- value indicates the probability of rejecting the null hypothesis. It is seen that only the p- value of the pair comparison HMPSO vs CMA-ES in terms of the mean value is larger than 0.05, which means that the null hypothesis is true (i.e., there is no significant difference between the HMPSO and CMA-ES), while the HMPSO has significant differences with the other eight algorithms in terms of both mean and SD values. Therefore, it can be concluded that the HMPSO has better comprehensive performances than MBA, KPSO, GPSO, PSO-TVAC, DMS-PSO and SPSO and has similar performances with the CMA-ES at a significance level  =0.05 .

5. Design optimization of an example engine accessory drive system An example here concerns the design problem of a real- life engine accessory drive system with eight accessory pulleys and a tensioner, as shown in Fig. 1. Original configuration parameters of the example system are the same as those in Tables 4 and 5 in Ref. [32]. An engine with four cylinders is 19

used to provide powers for all accessories. For a four-cylinder engine, its firing frequency can be calculated by 𝜔𝑓 =

2𝜋𝑛𝑐𝑠 60

4

𝜋𝑛𝑐𝑠

2

15

× =

, where 𝑛𝑐𝑠 is the rotational speed of the CS pulley. Because of

angular speed fluctuations of the CS pulley, the engine accessory drive system suffers from a series of harmonic excitations with frequencies of integer multiples of 𝜔𝑓 /2. Only the dominant excitation with a single frequency 𝜔𝑓 and amplitude 1.1342 deg in the harmonic excitations applied on the CS pulley are considered for simplicity of performing the HMPSO algorithm. Modal analysis of the example engine accessory drive system shows that the first natural frequency with a torsional vibration mode of pulleys is 51.03 Hz [32]. When the angular speed of the CS pulley is 1531 rpm, the engine excitation frequency is 𝜔𝑓 = 51 Hz, which almost equals to 51.03 Hz. Hence resonance occurs in the engine accessory drive system in this occasion. The HMPSO algorithm is thus used to minimize vibrations of the engine accessory drive system via configuration parameters optimization.

5.1. Parameters setting Constraints on the design variables are given by

80o  ˆt 0  180o  10 N  m / rad  Kˆ r  80 N  m / rad s. t.  0.02 m  lˆt  0.1 m  ˆ 10 N  m  Tt 0  60 N  m



The position vector of each particle is defined as xi = ˆt 0 , Kˆ r , lˆt , Tˆt 0

(27)



T

according to the design



u u u u u variables, and subsequently the design space are defined by its upper limit xi = ˆt 0 , Kˆ r , lˆt , Tˆt 0



l l l l l and lower limit xi = ˆt 0 , Kˆ r , lˆt , Tˆt 0



T

of



T

xi according to the constraints on the design variables.

According to the objective functions in Eqs. (10) and (11), the fitness functions of the HMPSO

20





algorithm are defined as f a ˆt 0 , Kˆ r , lˆt , Tˆt 0 





H1 /  2n  and fb ˆt 0 , Kˆ r , lˆt , Tˆt 0  H 2 , which

correspond to the objective functions 1 and  2 , respectively. A small fitness function value means a good particle position. Values of the weighting coefficients

ˆm , m  2, 3,

ˆj , j  1, 2,

,n ,

, n , and ˆt in H1 are all equal to 1. Parameters in the HMPSO algorithm are set

according to experience and numerous numerical experiments. A swarm with a population size 60 and maximum number of iterations 80 in the HMPSO algorithm is used here. Other parameters in the HMPSO algorithm are the same as those set in section 4.3.2.

5.2. Optimization results The HMPSO algorithm with the parameters set in section 5.1 is used here to solve the design problem of the example engine accessory drive system. The optimal solutions can be obtained by performing the implementation procedures of the HMPSO algorithm presented in section 4.2.3. Five other existing metaheuristic algorithms, including the LPSO [39], KPSO [40], CMA-ES [52], SPSO [36] and GA with use of the tournament selection strategy [53], are also used as comparisons to confirm good performances of the HMPSO. The constriction factor of the KPSO K is 0.729, the crossover probability of GA Pc is linearly decreasing with iterations from 0.9 to 0.4, and other initial parameters of the LPSO, CMA-ES, SPSO, GA and KPSO are the same as those of the HMPSO. All algorithms are programmed using MATLAB and the best solution is gained over 30 independent runs. Table 5 presents comparison of the best optimal solutions given by the six optimization algorithms using the objective functions 1 and  2 . As seen from Table 5, the best solution is found by the SPSO when the objective function 1 is used, while the best solution is found by the HMPSO and SPSO when the objective function  2 is used. The performance of the KPSO is the 21

worst when the objective function 1 is used and that becomes GA when the objective function

 2 is used. Statistical results provided by the six algorithms for this design problem are compared in terms of the worst, mean and best solutions as well as the SD values and numbers of fitness function evaluations (NFFEs) when the best optimal solution is found as shown in Table 6. It can be seen from Table 6 that the HMPSO with either the objective function 1 or  2 used provides the best solution with the minimum SD value compared with the other algorithms. The SPSO has similar performances in terms of the best and mean solutions with the HMPSO but converges a little slower than the HMPSO. The convergence rate of the HMPSO is lower than that of the LPSO and KPSO due to the hybrid mutation operators and higher than GA and CMA-ES. Note that premature convergence easily occurs in the LPSO and KPSO in this numerical example, which demonstrates the importance and effectiveness of the hybrid mutation operators in improving swarm diversity in the HMPSO in return. The HMPSO, CMA-ES and SPSO have close SD values with the same order of magnitude that are smaller than those of the other algorithms when the objective function 1 is used. The HMPSO and SPSO have close SD values with the same order of magnitude that are smaller than those of the remainder algorithms when the objective function  2 is used. Fig. 5 presents the fitness function values f a and fb versus the number of iterations obtained from the HMPSO using the objective functions 1 and  2 , respectively. It is observed that the HMPSO converges fast to the optimal solution with less than 10 iterations for the design optimization problem of the engine accessory drive system with either the objective function 1 or  2 being used. Vibration amplitudes of all the accessory pulleys, belt spans and the tensioner arm after optimization are presented in Table 7. All these five algorithms using either the objective function 22

1 or  2 demonstrate good optimization results: vibration amplitudes of all the accessory pulleys and belt spans as well as the tensioner arm are greatly reduced after optimization. As an objective function,  2 overall performs a little better than 1 in terms of torsional pulley vibrations, while

1 performs much better than  2 in terms of transverse belt vibrations. Dynamic responses of the engine accessory drive system before optimization are compared with those after optimization by the HMPSO algorithm with the objective function 1 being used. Angular displacements of all the accessory pulleys before and after optimization are shown in Figs. 6(a) and 6(b), respectively. From Fig. 6(a), it is seen that resonance takes places in the torsional responses of the accessory pulleys and tensioner arm with high torsional vibration amplitudes before optimization. The maximum peak-to-peak value of the vibration curve of each component in Fig. 6(a) is defined as its vibration amplitude. As shown in Figs. 6(a) and 6(b), the maximum vibration amplitude of all accessory pulleys and the tensioner arm is reduced from 29 mm to 3.4092 mm after optimization. Transverse vibrations of belt spans 7 and 8 before and after optimization are demonstrated in Figs. 7 and 8, respectively. It is seen that vibration amplitudes of the belt spans 7 and 8 decrease from 3.6476 mm to 0.2419 mm and from 8.178 mm to 0.2419 mm after optimization, respectively. The obvious improvement of the system vibration is in fact due to frequency detuning: the first natural frequency of the engine accessory drive system with a torsional mode of pulleys after optimization is changed to 36.6 Hz, which is tuned away from the engine excitation frequency 51 Hz. Dynamic tension fluctuation amplitudes in all belt spans are also enormously decreased due to system vibration improvement as shown in Figs. 9(a) and 9(b).

5.3. Robustness analysis The actual configuration of an engine accessory drive system may have a small difference with 23

the optimal one due to machining and installation errors, belt elongation, belt and bearing abrasion, and deformations of shafts, and so on, which are inevitable in practical usage. Robustness analyses are conducted here to investigate robustness of dynamic responses of the engine accessory drive system with respect to variations of the design variables relative to the optimal design in the design space. The robustness of the dynamic performances of the engine accessory drive system with respect to changes of 𝑙̂𝑡 is firstly analyzed as shown in Fig. 10(a). In Fig. 10(a), 𝑙̂𝑡 has deviations relative to the optimal value from −0.01 m to +0.04 m (i.e., 𝑙̂𝑡 varies from 0.01 m to 0.05 m) while the other structure parameters are kept at their optimal values. It is seen from Fig. 10(a) that both angular vibration amplitudes of the accessory pulleys and tensioner arm and transverse vibration amplitudes of the belt spans are sensitive to changes of 𝑙̂𝑡. Small deviations of 𝑙̂𝑡 (either increase or decrease) have a notably detrimental impact on rotational pulley vibrations but transverse belt vibrations monotonically increase fast with 𝑙̂𝑡. Similarly, the robustness of the dynamic performances of the engine accessory drive system with respective to changes of other structur e parameters including 𝜃̂𝑡0 , ̂𝑟 and 𝑇̂𝑡0 are analyzed as presented in Figs. 10(b)–(d), respectively. As shown in Fig. 10(b), 𝐾 deviations of 𝜃̂𝑡0 from −15° to +22° (i.e., 𝜃̂𝑡0 ranges from 93° to 130° ) have a significant effect on transverse belt vibrations but have a small effect on rotational pulley vibrations. It is seen from Figs. 10(c) and 10(d) that the dynamic performance of the engine accessory drive system has ̂𝑟 and 𝑇̂𝑡0 , respectively. Both rotational vibration strong robustness with respect to changes of 𝐾 amplitudes of all accessory pulleys and the tensioner arm and transverse vibration amplitudes of the ̂𝑟 as shown in Fig. 10(c). Although both belt and pulley vibration belt spans slowly increase with 𝐾 amplitudes almost remain constant when 𝑇̂𝑡0 changes from 20 N·m to 60 N·m as shown in Fig. 24

10(d), increasing 𝑇̂𝑡0 will increase belt tension. This means that 𝑇̂𝑡0 is flexibly adjustable in the design space, but it must be emphasized that too small or too large 𝑇̂𝑡0 should be avoided to prevent belt slip on pulleys and prolong belt usage life. From the above analyses, it can be concluded that the dynamic performance of the engine ̂𝑟 , and not sensitive to accessory drive system is mostly sensitive to 𝑙̂𝑡 and 𝜃̂𝑡0 , a little sensitive to 𝐾 ̂𝑟 𝑇̂𝑡0 . Therefore, the optimal configuration parameters 𝑙̂𝑡 and 𝜃̂𝑡0 should be strictly assured; both 𝐾 and 𝑇̂𝑡0 have a wide range to adjust in the design space and thus can be flexibly adjusted according to its practical working conditions.

6. Conclusion A good design of the engine accessory drive system requires minimizing rotational vibrations of all accessory pulleys and the tensioner arm and transverse vibrations of all belt spans at the same time. It is a multi-objective optimization problem according to the requirements of the design problem. Unlike the traditional design method that greatly relies on experience and experiments, a HMPSO algorithm is proposed here for the first time to solve the design problem of the engine accessory drive system. Based on our previous works, key structure parameters of the tensioner ̂𝑟, 𝑇̂𝑡0 , 𝜃̂𝑡0 and 𝑙̂𝑡 are defined as design variables. The WSM and MTM are applied in including 𝐾 the HMPSO algorithm to convert the multi-objective optimization problem of the engine accessory drive system to a single-objective one with corresponding objective functions 1 and  2 , respectively. In addition, a modified inertia factor is used in the HMPSO algorithm to balance the convergence rate and space search ability during the iteration processes and a hybrid mutation strategy with six hybrid mutation operators is introduced to diversify the swarm and avoid premature 25

convergence. Ten benchmark functions are adopted to assess the behavior of the proposed HMPSO algorithm and the non-parametric statistical test results validate the superiority of the HMPSO algorithm against several other concerned existing metaheuristic algorithms in terms of solution quality and stability. The HMPSO algorithm is then used to solve the design optimization problem of an example engine accessory drive system and the results are compared with those obtained from other five existing metaheuristic optimization algorithms. Comparison demonstrates that the HMPSO algorithm provides better optimal solutions with less SD values than other algorithms except the SPSO when either the objective function 1 or  2 is used. Hence the overall superiority of the HMPSO algorithm in global optima search ability and stability is verified again. The objective function  2 overall performs a little better than the objective function 1 in terms of rotational pulley vibrations, while 1 performs much better than  2 in terms of transverse belt vibrations. Considering of the fact that the actual configuration of an engine accessory drive system may have a small difference with the optimal one, a robust analysis is conducted for the first time to calculate robustness of dynamic responses of the engine accessory drive system with respect to small variations of the design variables relative to the optimal design in the design space. It is demonstrated that 𝑙̂𝑡 and 𝜃̂𝑡0 have a great impact on dynamic performances of the engine accessory ̂𝑟 and 𝑇̂𝑡0 . Hence, it is drive system, while the system vibrations are not sensitive to variations of 𝐾 suggested that the optimal configuration parameters 𝑙̂𝑡 and 𝜃̂𝑡0 should be strictly assured, while ̂𝑟 and 𝑇̂𝑡0 can be flexibly adjusted in the design space according to the practical working 𝐾 conditions of the engine accessory drive system. The methods and conclusions presented in this work can provide well guidelines for designers and engineers in the design of an engine accessory drive 26

system.

Acknowledgments The authors gratefully acknowledge the support of the National Natural Science Foundation of China (Nos. 51805339 and 11802188) and Fundamental Research Funds for Central Universities in China (No. YJ201825).

Appendix. Formulations of the last five Benchmark functions f6

f10

Ackley’s Function:

 1 D 2 1 D  f 6  x   20exp  0.2 xi   exp   cos  2 xi    20  e    D i 1   D i 1  

(A.1)

HappyCat Function:

f7  x  

D

x D i 1

1/4

2 i



D D 1  2 0.5 x  xi   0.5   i  D i 1 i 1 

(A.2)

Weierstrass Function: km  km k  k   f8  x       a cos  2 b  xi  0.5     D  a k cos  2 b k  0.5 , a  0.5, b  0.3, km  20 (A.3) i 1  k  0 k 0  D

Modified Schwefel’s (MS) Function:





 1/2 , if zi  500  zi sin zi  2 zi  500    g  zi     500  mod  zi ,500   sin 500  mod  zi ,500   , if zi  500 , 10000 D  2  zi  500      mod  zi ,500   500 sin  mod  z i ,500   500   , if zi  500   10000 D 









D

f9  x   4 1 8 . 9 8 2 D 9  g i z i 1

Expanded Scaffer’s F6 (ESF6) Function:

27

,z  i x 

i

4 2 0 . 9 6 8 7 4 6 2 2 7 5 0 3(A.4) 6

g  x, y   0.5 

sin 2





x 2  y 2  0.5

1  0.001 x

2

 y2 



2

,

f10  x   g  x1 , x2   g  x2 , x3   ...  g  xD 1 , xD   g  xD , x1 

(A.5)

References [1] R.L. Cassidy, S.K. Fan, R.S. MacDonald, W.F. Samson, Serpentine extended life accessory drive. SAE paper (1979) No. 790699. [2] W.-B. Shangguan, X.-K. Zeng, Modeling and Validation of Rotational Vibration Responses for Accessory Drive Systems—Part II: Simulations and Analyses, ASME Journal of Vibration & Acoustics 135 (2013) 031003/1-13. [3] C.-J. Shieh, W.-H. Chen, Effect of angular speed on behavior of a V-belt drive system, International Journal of Mechanical Sciences 44 (2002) 1879-1892. [4] L. Kong, R.G. Parker, Equilibrium and Belt-pulley Vibration Coupling in Serpentine Belt Drives, ASME Journal of Applied Mechanics 70 (2003) 739-749. [5] R.G. Parker, Efficient Eigensolution, Dynamic Response, and Eigensentivity of Serpentine Belt Drives, Journal of Sound and Vibration 270 (2004) 15-38. [6] M. Nouri, J.W. Zu, Dynamic analysis and optimization of tensioner in automotive serpentine belt drive systems, ASME 2002 Design engineering technical conferences and computer and information in engineering conference, Montreal, Canada. [7] Z. Li, L. Zheng, Y. Ren, et al., Multi-objective optimization of active suspension system in electric vehicle with In-Wheel-Motor against the negative electromechanical coupling effects, Mechanical Systems and Signal Processing 116 (2019) 545-565. [8] J.S. Jensen, N.L. Pedersen, On maximal eigenfrequency separation-material in two- material structures: the 1D and 2D scalar cases, Journal of Sound and Vibration 289 (2006) 967-986. [9] J. He, X. Jin, S.Y. Xie, et al., Multi-body dynamics modeling and TMD optimization based on 28

the improved AFSA for floating wind turbines, Renewable Energy 141 (2019) 305-321. [10] H. Zhu, W.D. Zhu, Y. Hu, X. Wang, Periodic Response of a Timing Belt Drive System With an Oval Cogged Pulley and Optimal Design of the Pitch Profile for Vibration Reduction, ASME Journal of Computational and Nonlinear Dynamics 13(1) (2017) 011014/1-13. [11] Z.C. Hou, Y.X. Lao, Q.H Lu, Sensitivity analysis and parameter optimization for vibration reduction of undamped multi- ribbed belt drive systems, Journal of Sound and Vibration 317 (2008) 591–607. [12] H. Zhu, Y. Hu, W.D. Zhu, Modification of natural frequencies of an automotive belt drive system based on eigen-sensitivity analysis of its configuration parameters, European Journal of Mechanics A/Solids, 67 (2018) 137-156. [13] J.H. Holland, Outline for a logical theory of adaptive systems, Journal of the Association for Computing Machinery 9 (3) (1962) 297-314. [14] J.D. Bagley, The behavior of adaptive systems which employ genetic and correlation algorithms, Dissertation Abstracts International, University of Michigan, 1967. [15] J. Kennedy, R. Eberhart, Particle Swarm Optimization, IEEE International Conference on Neural Networks (ICNN), Perth, WA, Australia, 4 (1995) 1942-1948. [16] X.S. Yang, Firefly algorithms for multimodal optimization, Stochastic algorithms: foundations and applications, Springer Berlin Heidelberg, 5792 (2009) 169-178. [17] S. Kirkpatrick, C.D. Gelatt, M.P. Vecchi, Optimization by simulated annealing, Science 220 (1983) 671–680. [18] R.V. Rao, V.J. Savsani, D.P. Vakharia, Teaching- learning-based optimization: A novel method for constrained mechanical design optimization problems, Computer-Aided Design 43 (2011) 303-315. [19] A. Sadollah, A. Bahreininejad, H. Eskandar, M. Hamdi, Mine blast algorithm: a new population 29

based algorithm for solving constrained engineering optimization problems, Applied Soft Computing 13 (2013) 2592–2612. [20] M. Gouttefarde, C.M. Gosselin, Analysis of the wrentch-closure workspace of planar parallel cable-driven mechanisms, IEEE Transactions on Robotics 22 (2006) 434-445. [21] A.A. Nagra, F. Han, Q.H. Ling, An improved hybrid self- inertia weight adaptive particle swarm optimization algorithm with local search, Engineering Optimization 517(7) (2019) 1115-1132. [22] X.Y. Kou, G.T. Parks, S.T. Tan, Optimal design of functionally graded materials using a procedural model and particle swarm optimization, Computer-Aided Design 44(4) (2012) 300-310. [23] A. Mortazavi, V. Toğan, Simultaneous size, shape, and topology optimization of truss structures using integrated particle swarm optimizer, Structural and Multidisciplinary Optimization 54(4) (2016) 715-736. [24] T.B. Joshua, J. Xin, K.A. Sunil, Optimal design of cable-driven manipulators using particle swarm optimization, ASME Journal of Vibration & Acoustics 8 (2016) 041003/1-8. [25] D. Novitasari, I. Cholissodin, W.F. Mahmudy, Hybridizing PSO With SA for Optimizing SVR Applied to Software Effort Estimation, TELKOMNIKA Indonesian Journal of Electrical Engineering 14(1) (2016) 245-253. [26] Q. He, L. Wang, A hybrid particle swarm optimization with a feasibility-based rule for constrained optimization, Applied Mathematics and Computation 186 (2007) 1407-1422. [27] J. Wang, Z. Yin, A ranking selection-based particle swarm optimizer for engineering design optimization problems, Structural and Multidisciplinary Optimization 37 (2008) 131-147. [28] J.-J. Lei, J. Li, A modified particle swarm optimization for practical engineering optimization, 2009 Fifth International Conference on Natural Computation, 3 (2009) 177-180. [29] L.D.S. Coelho, Gaussian quantum-behaved particle swarm optimization approaches for 30

constrained engineering design problems, Expert Systems and Applications 37(2) (2010) 1676-1683. [30] M.J. Mahmoodabadi, M. Bisheban, An online optimal linear state feedback controller based on MLS approximations and a novel straightforward PSO algorithm, Transactions of the Institute of Measurement and Control 36(8) (2014) 1132-1142. [31] H. Garg, A hybrid PSO–GA algorithm for constrained optimization problems, Applied Mathematics and Computation 274 (2016), 292 – 305. [32] H. Zhu, Y. Hu, W.D. Zhu, Dynamic responses of an engine front-end accessory belt drive system with pulley eccentricities via two spatial discretization methods, Proceedings of the Institution of Mechanical Engineers, Part D: Journal of Automobile Engineering 232(4) (2018) 482-498. [33] R.G. Parker, Y. Lin, Parametric instability of axially moving media subjected to multifrequency tension and speed fluctuations, ASME Journal of Applied Mechanics 68 (2001) 49-57. [34] C. Fonseca, P. Fleming, An overview of evolutionary algorithms in multiobjective optimization, Evolutionary Computation 3 (1995) 1-18. [35] T. Ray, P. Saini, Engineering design optimization using a swarm with an intelligent information sharing among individuals, Engineering Optimization 33 (2001) 735-748. [36] M.J. Mahmoodabadi and M. Bisheban, An online optimal linear state feedback controller based on MLS approximations and a novel straightforward PSO algorithm, Transactions of the Institute of Measurement and Control 36 (2014) 1132–1142. [37] S. He, E. Prempain, Q. H. Wu, An improved particle swarm optimizer for mechanical design optimization problems, Engineering Optimization 36 (2004) 585-605. [38] N.B. Guedria, Improved accelerated PSO algorithm for mechanical engineering optimization problems, Applied Soft Computing 40 (2016) 455-467.

31

[39] Y. Shi, R. Eberhart, Empirical study of particle swarm optimization. IEEE International Conference on Evolutionary Computation, Washington, USA, 1999, 1945-1950. [40] M. Clerc, The Swarm and The Queen: Towards a deterministic and Adaptive Particle Swarm Optimization, IEEE Congress on Evolutionary Computation (CEC), 1999. [41] A.A.A. Esmin, S. Matwin, HPSOM: A Hybrid Particle Swarm Optimization Algorithm with Genetic Mutation, International Journal of Innovative Computing, Information and Control 9(5) (2013) 1919-1934. [42] A.T. Dosdoğru, M. Göçken, F. Geyik, Integration of genetic algorithm and Monte Carlo to analyze the effect of routing flexibility, International Journal of Advanced Manufacturing Technology 81(5-8) (2015) 1379-1389. [43] G. Zhou, Z.D. Ma, A.G. Cheng, G.Y. Li, J. Huang, Design optimization of a runflat structure based on multi–objective genetic algorithm, Structural and Multidisciplinary Optimization 51(6) (2015) 1363-1371. [44] S.H. Ling, H.H.C. Iu, K.Y. Chan, H.K. Lam, B.C.W. Yeung and F.H. Leung, Hybrid particle swarm optimization with wavelet mutation and its industrial applications, Systems, Man, and Cybernetics, Part B: Cybernetics, IEEE Transactions on 38(3) (2008) 743-763. [45] N. Kumar, D.P. Vidyarthi, A novel hybrid PSO-GA meta-heuristic for scheduling of DAG with communication on multiprocessor systems, Engineering with Computers 32(1) (2016) 35-47. [46] R. Caponetto, L. Fortuna, S. Fazzino, Chaotic sequencesto improve the performance of evolutionary algorithms, IEEE International Conference on Evolutionary Computation 7(3) (2003) 289-304. [47] J.-B. Cheng, Q.-W. Ye, Y. Zhou, X.-H. Cao, New Multi-Mutation Particle Swarm Optimization algorithm, Computer Engineering and Applications (in Chinese) 43(7) (2007) 59-61. [48] J. Kennedy, Stereotyping: Improving particle swarm performance with cluster analysis, IEEE International Conference on Evolutionary Computation 2 (2000) 1507-1512. 32

[49] A. Ratnaweera, K.H. Saman algamuge, C.W. Harry, Self-organizing hierarchical particle swarm optimizer with time- varying acceleration coefficients, IEEE Transactions on Evolutionary Computation 8(3) (2004) 240-255. [50] Y. Shi, R.C. Eberhart, A modified particle swarm optimizer, IEEE International Conference on Evolutionary Computation (1998) 69-73. [51] J.J. Liang, P.N. Suganthan, Dynamic multi- swarm particle swarm optimizer with local search. IEEE International Conference on Evolutionary Computation 1 (2005) 522-528. [52] N. Hansen, S.D. Muller, P. Koumoutsakos, Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES), Evolutionary Computation 11(1) (2003) 1–18. [53] J. Derrac, S.García, D. Molina, F. Herrera. A practical tutorial on the use of nonparametric statistical tests as a methodology for comparing evolutionary and swarm intelligence algorithms, Swarm and Evolutionary Computation 1(1) (2011) 3-18. [54] L.M. Li, K.D. Lu, G.Q. Zeng, L. Wu, M.R. Chen, A novel real-coded population-based extremal optimization algorithm with polynomial mutation: A non-parametric statistical study on continuous optimization problems, Neurocomputing 174 (2016) 577-587.

33

θˆ5

Y PSP

yˆ 5

X θˆ 6

xˆ 6 yˆ 6

Jˆ 7

WP

wˆ i 1

θˆ 7

θˆ3

Jˆ 6 Idler 2

wˆ i

β1

xˆ 4

Idler 1 Jˆ

xˆ3

3

θˆ8

Jˆ8 Ten

yˆ 3 yˆ 4

ALT Jˆ 4

β2 θt 0

xˆ8 yˆ θt 8

xˆ 7

yˆ 7

Jˆ5

xˆ5

xˆ 2

θˆ 4

yˆ 2

CS θˆ1

Jˆ1

yˆ1

AC xˆ1

Jˆ 2

θˆ 2

cˆ

Fig. 1 Schematic of a typical accessory belt drive system: CS—crank shaft pulley, AC—air conditioner pulley, PSP—power steering pump pulley, WP—water pump pulley, and Ten—tensioner pulley. [25]

34

Start Parameter setting

Initialize a swarm

=1 Update xi    and v i    via Eq. (15), (16), and (18)

 1 No

Evaluate the fitness

Output optimal Yes results

Meet termination criterion?

New swarm

Update Pi    and Pg   

Update Pi    and Pg   

Mutation

Fig. 2 Flowchart of the HMPSO algorithm.

(a) x 10

6

600

Function values

Sphere Quadric Rosenbrock

3

2

1

0 0

Quadric Noise Griewark HappyCat

500

4

Function values

5

(b)

400 300 200 100

200

400

600

800

1000

0 0

1200

Number of iterations

2000

4000

6000

8000

Number of iterations

Fig. 3 Convergence histories of the HMPSO algorithm on benchmark functions.

35

10000

(a)

(b)

100

100

particle 1 particle 3 particle 10

90

3

80

70

Function values 10

Function values 10

3

80

60 50 40 30

70 60 50 40 30

20

20

10

10

0 0

particle 1 particle 3 particle 10

90

10

20

30

40

50

60

70

80

90

100

0 0

10

20

Number of iterations

30

40

50

60

70

80

90

100

Number of iterations

Fig. 4(a) Behaviors of the particles in the basic PSO and (b) the HMPSO.

1.4

2.1

fa

b

2

fb 1.2

1.9

1.1

1.8

1

1.7

0.9 0

10

20

30

40

50

60

70

Fitness function value, f

Fitness function value, f

a

1.3

1.6 80

Number of iterations Fig. 5 Values of the fitness functions f a and fb for the best particle position obtained from the HMPSO algorithm with the objective functions 1 and  2 versus the number of iterations, respectively.

36

(a)

(b) 2.5

Angular displacements of accessories (mm)

Angular displacements of accessories (mm)

15

10

5

0

-5

-10

-15 0

AC ALT PSP WP Ten Arm

0.05

0.1

0.15

0.2

0.25

0.3

0.35

2 1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5 0

0.4

AC ALT PSP WP Ten Arm

0.02

0.04

0.06

0.08

0.1

0.12

Time (s)

Time (s)

Fig. 6 Angular displacements of accessory pulleys and the tensioner arm: (a) before optimization and (b) after optimization.

37

3 2 1 0 -1 0

-2 50 100 150

400

) m (m

7

-3 380

420

Time (m s)

200

440

460

x

Transverse displacement (mm)

(a)

250

Transverse displacement (mm)

(b)

0.1 0.05 0 -0.05 -0.1 -0.15 0

0 10

50 20

100 30

40 Time (ms)

150 50

60

200 70

)

x7 (mm

250

Fig. 7 Transverse displacement of the belt span 7: (a) before and (b) after optimization.

38

Transverse displacement (mm)

(a)

5 4 3 2 1 0 -1 -2 -3 -4 -5 380

100

390

400

410

420

25 430

440

Time (ms)

450

460

x8

50

(m m )

75

0

Transverse displacement (mm)

(b) 0.1 0.05 0 -0.05 -0.1 -0.15 100 75

Tim 50 e( ms )

25 0

0

10

20

30

40

50

60

x8 (mm)

Fig. 8 Transverse displacement of the belt span 8: (a) before and (b) after optimization.

39

(a)

(b)

8000

2600

6000

2400 2200

Dynamic tension (N)

4000

Dynamic tension (N)

belt span 1 belt span 2&3 belt span 4 belt span 5&6 belt span 7&8

2000 0 -2000 -4000 -6000 -8000 0

belt span 2 belt span 2 & 3 belt span 4 belt span 5 & 6 belt span 7 & 8 0.05 0.1

2000 1800 1600 1400 1200 1000 800

0.15

0.2

0.25

0.3

0.35

0

Time (s)

0.05

0.1

0.15

Time (s)

Fig. 9 Dynamic tensions in belt spans: (a) before optimization, and (b) after optimization.

40

0.2

(b)

Optimal value 0.2

1

0.15

0.1

0.5

0.05

0 0.01

0.015

0.02

0.025

0.03

0.035

0.04

0.045

0 0.05

0.9

0.06

Optimal value

0.05

0.03

0.3

0.02

0.15

0.01

95

100

105

110

AC Idler 1 ALT Idler 2

Span 7 Span 8

0.05

0.55

0.5

0.04

0.45

0.03

0.4

0.02

0.35 Optimal value 35

40

45

50

Torsional stiffness of tensioner, Kˆ r

55

Angular vibration amplitudes of pulleys (rad)

PS WP Ten Arm

Vibration amplitudes of belt spans (mm)

Angular vibration amplitudes of pulleys (rad)

0.6

30

115

120

125

Installation angle of tensioner arm, ˆt 0

0 130

(d)

0.07

25

0.6

0.45

(c)

0.01 20

0.75

0.04

Length of tensioner arm, lˆt

0.06

WP Ten Arm Span 7 Span 8

AC Idler 1 ALT PS Idler 2

0.065

0.055

PS WP Ten Arm

AC Idler 1 ALT Idler 2

Span 7 Span 8

0.045

0.43

0.41

0.035

0.39 Optimal value

0.025

0.015 20

0.3 60

0.45

0.37

25

30

35

40

45

50

Initial torque of tensioner, Tˆ0

55

Vibration amplitudes of belt spans (mm)

0.25

0.07

Vibration amplitudes of belt spans (mm)

Span 7 Span 8

AC Idler 1 ALT PS Idler 2 WP Ten Arm

Angular vibration amplitudes of pulleys (rad)

1.5

Vibration amplitudes of belt spans (mm)

Rotational vibration amplitudes of pulleys (rad)

(a)

0.35 60

Fig. 10 Robustness of rotational vibrations of accessory pulleys and the tensioner arm, and transverse vibrations of belt spans to changes of design variables: (a) length of the tensioner arm 𝑙̂𝑡, (b) ̂𝑟, and (d) initial installation angle of the tensioner arm 𝜃̂𝑡0 , (c) torsional stiffness of the tensioner 𝐾 torque of the tensioner 𝑇̂𝑡0 .

41

Table 1 Selection strategy of a mutation operator

p   Mutation operator

0  p1

p1  p2

p2  p3

p3  p4

p4  p5

p5  1

(a)

(b)

(c)

(d)

(e)

(f )

Table 2 Benchmark functions Functions

Formula D

Sphere

f1  x    xi2

Griewank

f2  x  

f  xo 

xo

Search range D

0, 0,

, 0

0

D

0, 0,

, 0

0

D

0, 0,

, 0

0

 1.28, 1.28

0, 0,

, 0

0

2 2 f5  x    100  xi 1  xi2    xi  1   i 1 

10, 10

1, 1,

Ackley

f6  x 

10, 10

0, 0,

, 0

0

HappyCat

f7  x 

 100, 100

D

0, 0,

, 0

0

Weiertrass

f8  x 

 100, 100

D

0, 0,

, 0

0

MS

f9  x 

 100, 100

D

0, 0,

, 0

0

ESF6

f10  x 

 100, 100

D

0, 0,

, 0

0

Quadric

 100, 100

i 1





1 D 2 D cos xi / i  1  xi   4000 i 1 i 1

D  i  f3  x      x 2j  i 1  j 1 

 100, 100

D

Quadric Noise

f 4  x    i xi4  rand

D

i 1

D 1

Rosenbrock

 600, 600

D

D

, 1

0

xo is the theoretical optimal solution and f  x o  is the corresponding function value.

42

Table 3 Comparison of the statistical results given by different optimization algorithms for the test functions Function

CMA-ES

MBA

KPSO

GPSO

PSO-T VAC

DMS-PSO

SPSO

HMPSO

Mean

4.50E-22(7)

4.03E-06(8)

4.47E-29(6)

1.98E-53(5)

8.36E-108(3)

3.85E-54(4)

5.20E-143(2)

1.62E-157(1)

SD

2.01E-21(7)

1.93E-06(8)

1.13E-28(6)

7.08E-53(5)

2.97E-107(3)

1.75E-53(4)

1.61E-142(2)

4.42E-157(1)

Mean

0(1)

1.49E-02(5)

1.10E-02(3)

2.37E-02(8)

1.91E-02(6)

1.31E-02(4)

2.35E-02(7)

0(1)

SD

0(1)

2.33E-02(7)

1.60E-02(4)

2.57E-02(8)

1.50E-02(3)

1.73E-02(5)

2.27E-02(6)

0(1)

Mean

4.49E-08(4)

2.34E-03(5)

1.86E+01(7)

6.45E-02(6)

3.68E-09(3)

4.75E+01(8)

1.47E-22(1)

3.26E-09(2)

SD

2.01E-07(4)

8.69E-03(5)

3.07E+01(7)

9.46E-02(6)

5.01E-09(3)

5.64E+01(8)

4.18E-22(1)

2.72E-09(2)

Quadric

Mean

5.27E-03(1)

9.84E-03(5)

1.49E-02(7)

7.77E-03(4)

6.89E-03(3)

1.10E-02(6)

3.40E-02(8)

6.10E-03(2)

Noise

SD

1.21E-02(7)

3.73E-03(4)

5.66E-03(6)

2.42E-03(1)

3.71E-03(3)

3.94E-03(5)

2.08E-02(8)

3.34E-03(2)

Rosenbrock

Mean

1.05E+0(1)

3.54E+01(8)

2.18E+01(5)

2.81E+01(6)

1.56E+01(4)

3.23E+01(7)

1.60E+0(2)

1.73E+0(3)

SD

4.70E+0(3)

3.62E+01(8)

1.15E+01(4)

2.46E+01(7)

2.41E+01(5)

2.41E+01(5)

2.06E+0(1)

2.39E+0(2)

Mean

1.22E-10(2)

2.00E+01(4)

2.00E+01(4)

2.00E+01(4)

2.00E+01(4)

2.54E+0(3)

2.00E+01(4)

4.91E-15(1)

SD

5.72E-11(3)

0(1)

1.81E-06(6)

3.78E-08(5)

2.38E-08(4)

4.70E+0(8)

2.74E-06(7)

1.54E-15(2)

Mean

4.68E-01(3)

6.20E-01(6)

3.33E+02(8)

5.84E-01(5)

2.76E-01(2)

1.21E+0(7)

5.73E-01(4)

2.40E-01(1)

SD

9.95E-02(2)

1.41E-01(4)

1.84E+02(8)

1.47E-01(5)

1.26E-01(3)

2.88E+0(7)

1.75E-01(6)

5.54E-02(1)

Mean

3.28E+0(8)

2.14E-13(5)

2.27E-13(6)

1.76E-13(4)

1.69E-13(3)

0(1)

1.89E-03(7)

3.55E-16(2)

SD

9.34E+0(8)

1.95E-13(6)

4.22E-14(4)

3.18E-14(3)

4.44E-14(5)

0(1)

6.55E-03(7)

1.59E-15(2)

Mean

3.82E-04(1)

3.82E-04(1)

2.11E+03(8)

2.53E+02(7)

3.82E-04(1)

8.31E+0(6)

3.83E-04(5)

3.82E-04(1)

SD

8.54E-13(2)

2.55E-07(5)

1.22E+03(8)

3.38E+02(7)

2.09E-11(3)

1.89E+01(6)

2.44E-07(4)

3.64E-13(1)

Mean

1.30E-03(2)

1.14E+01(5)

1.17E+01(6)

1.43E+01(8)

9.11E+0(3)

9.53E+0(4)

1.30E+01(7)

0(1)

SD

4.93E-03(2)

3.11E-01(4)

7.64E-01(6)

2.39E-01(3)

9.52E-01(8)

8.81E-01(7)

5.74E-01(5)

0(1)

(3)

(5.2)

(6)

(5.7)

(3.2)

(5)

(4.7)

(1.5)

(3.9)

(5.2)

(5.9)

(5)

(4)

(5.6)

(4.7)

(1.5)

Sphere

Griewark

Quadric

Ackley

HappyCat

Weierstrass

MS

ES-F6

Average rank Mean SD

Table 4 Wilcoxon signed-rank test results Mean

SD

Comparison

R

R

p-value

Comparison

R+

R−

p-value

HMPSO vs CMA-ES

7

29

0.119

HMPSO vs CMA-ES

0

45

0.007

HMPSO vs MBA

0

45

0.007

HMPSO vs MBA

1

54

0.007

HMPSO vs KPSO

0

55

0.005

HMPSO vs KPSO

1.5

53.5

0.005

HMPSO vs GPSO

0

55

0.005

HMPSO vs GPSO

0

66

0.008

HMPSO vs PSO-TVAC

0

45

0.007

HMPSO vs PSO-TVAC

0

55

0.004

HMPSO vs DMS-PSO

1

54

0.007

HMPSO vs DMS-PSO

1

54

0.006

HMPSO vs SPSO

4

51

0.016

HMPSO vs SPSO

4

51

0.016

+

−

43

Table 5 Best solutions given by the five optimization algorithms

1

2

a

OD

a

DV HMPSO

GA

SPSO

CMA-ES

LPSO

KPSO

HMPSO

SPSO

GA

𝜃̂𝑡0

108.8416

112.5958

108.5521

108.7035

108.842

107.8019

97.6544

97.6543

̂𝑟 𝐾

22.812

23.5236

22.9836

20

25.906

42.9783

20

𝑙̂𝑡

0.022

0.0232

0.0222

0.02316

0.022

0.0222

𝑇̂𝑡0

44.9311

22.5895

25

25.1612

43.1883

𝑓𝑎/𝑓𝑏

0.9417

0.9482

0.9414

0.9438

0.9435

a

CMA-ES

LPSO

103.9605

90

97

87.1210

130

20

22.4737

20

20

21.6155

48

0.025

0.025

0.0241

0.026

0.026

0.0269

0.05

47.3673

25.02

25

21.9124

25

60

31.4755

30

0.9548

1.6169

1.6169

1.6296

1.6171

1.6190

1.6212

-

DV and aOD represent the design variables and original design, respectively.

Table 6 Statistical results given by the five optimization algorithms Algorithms

Worst

Mean

Best

SD

NFFEs

GA

0.9847

0.9636

0.9482

2.18E-02

3200

CMA-ES

0.9466

0.9442

0.9438

8.02E-04

3600

LPSO

0.9774

0.9532

0.9435

9.60E-03

2460

KPSO

0.9781

0.9695

0.9548

1.46E-02

2220

SPSO

0.9426

0.9419

0.9414

6.03E-04

3280

HMPSO

0.9421

0.9420

0.9417

3.24E-04

2761

GA

1.8640

1.6431

1.62960

8.02E-02

3300

CMA-ES

1.6223

1.61923

1.61714

6.18E-04

2640

LPSO

1.9022

1.6218

1.61900

5.50E-03

2160

KPSO

1.7325

1.6520

1.61912

5.41E-02

2100

SPSO

1.61733

1.61707

1.61693

1.06E-05

2590

HMPSO

1.61741

1.61706

1.61692

3.18E-05

2408

1

2

44

KPSO

Table 7 Statistical vibration amplitudes of all accessory pulleys and belt spans, and the tensioner arm that result from optimization of the five optimization algorithms Algorithms

1

2

AC

ALT

Idle 1

PSP

Idle 2

WP

Ten

Arm

Span 7

Span 8

GA

1.3275

1.5805

1.6573

1.4671

1.3235

1.2795

0.8932

0.9290

0.1507

0.3900

CM A-ES

1.3266

1.5781

1.6546

1.4627

1.3180

1.2699

0.8892

0.9232

0.2619

0.2535

LPSO

1.3548

1.6252

1.7061

1.4735

1.2879

1.2007

0.8799

0.8306

0.2433

0.2432

KPSO

1.3528

1.6250

1.7065

1.4893

1.3203

1.2414

0.9022

0.8689

0.2852

0.2229

SPSO

1.3493

1.6156

1.6955

1.4669

1.2893

1.2063

0.8788

0.8431

0.2556

0.2387

HM PSO

1.3543

1.6240

1.7046

1.4713

1.2833

1.1954

0.8777

0.8267

0.2419

0.2419

GA

1.3105

1.5545

1.6296

1.4734

1.3675

1.3703

0.9155

1.0066

0.4608

0.1240

CM A-ES

1.3001

1.5396

1.6171

1.4845

1.4180

1.4312

0.9407

1.1138

0.9627

0.3795

LPSO

1.3005

1.5406

1.6190

1.4864

1.4217

1.4351

0.9426

1.1158

1.0212

0.4026

KPSO

1.3081

1.5510

1.6259

1.4759

1.3791

1.3844

0.9279

1.0928

1.0928

0.500

SPSO

1.3013

1.5410

1.6169

1.4824

1.4188

1.4218

0.9320

1.0667

0.7120

0.1019

HM PSO

1.3013

1.5410

1.6169

1.4824

1.4100

1.4218

0.9320

1.0677

0.6793

0.0972

5.6763

10.4528

11.6896

12.810

13.8330

14.4596

7.1336

4.1725

3.6476

8.178

BOa a

BO represents the vibration amplitudes of all pulleys and belt spans 7 and 8 before optimization

45