Multi-Objective Optimization of Planetary Gearbox with Adaptive Hybrid Particle Swarm Differential Evolution Algorithm

Abstract

This paper considers the problem of constrained multi-objective non-linear optimization of planetary gearbox based on hybrid metaheuristic algorithm. Optimal design of planetary gear trains requires simultaneous minimization of multiple conflicting objectives, such as gearbox volume, center distance, contact ratio, power loss, etc. In this regard, the theoretical formulation and numerical procedure for the calculation of the planetary gearbox power efficiency has been developed. To successfully solve the stated constrained multi-objective optimization problem, in this paper a hybrid algorithm between particle swarm optimization and differential evolution algorithms has been proposed and applied to considered problem. Here, the mutation operators from the differential evolution algorithm have been incorporated into the velocity update equation of the particle swarm optimization algorithm, with the adaptive population spacing parameter employed to select the appropriate mutation operator for the current optimization condition. It has been shown that the proposed algorithm successfully obtains the solutions of the non-convex Pareto set, and reveals key insights in reducing the weight, improving efficiency and preventing premature failure of gears. Compared to other well-known algorithms, the numerical simulation results indicate that the proposed algorithm shows improved optimization performance in terms of the quality of the obtained Pareto solutions.

1. Introduction

Planetary gearboxes have a wide application in various mechanical systems, such as industrial drives, rotorcraft, automobiles, wind turbines, etc., where they can offer compact dimensions and higher power densities with less noise and higher torque-to-weight ratios, especially compared to standard parallel axis gear trains [1,2]. However, the design of such gearboxes involves multiple planet branches, which also reduces efficiency, where the multiple gear mesh regions in planetary gear sets dictate the overall efficiency. The main goal is to reduce the weight and power loss of the construction and produce gearbox with enhanced service life of the components. In this regard, in recent years a lot of research has been directed towards the optimal design of the gear trains [3,4]. In this way, researchers have considered the weight minimization of the single stage spur gears under the constraints, such as gear bending strength and shafts torsional strength [5]. Based on the obtained optimization results authors concluded that the ratio between weight of gears and space area has been improved by 48.8% and 33.6%, respectively, which confirmed the efficiency of the proposed method. Similarly, authors considered the weight optimization of the helical gear pair, where the design variables included: module of the gears, helix angle, number of teeth of the pinion gear and the width of gears [6]. The obtained results indicated that the proposed optimization method has improved performance compared to conventional gradient-based optimization algorithms. Furthermore, authors in [3] formulated the optimization problem on the single stage gearbox, where the objective function is the volume of gears, shaft and bearings. More recently, authors in [4] proposed a method to minimize the volume of the planetary gear trains with the aim to reduce space and mass. The proposed method successfully obtained the optimal value in benchmark test functions and the best parameters regarding the design of planetary gear trains compared to the well-known designs in the literature.

Majority of the above-mentioned papers deal with the problem of single-objective optimization (SOO). However, many practical engineering optimization problems have one or more conflicting objectives, which need to be simultaneously optimized. Regarding the gear train optimization, it is important to analyze the influence of the different contradictory objectives, such as weight or volume, power loss, center distance etc. Therefore, the multi-objective optimization (MOO) problems have received significant attention especially presently in the field of mechanical engineering applications, rotor-dynamic, electrical machine design and wireless communications, etc., [7,8,9,10,11,12]. It is a challenging task, which motivated researchers to have a growing interest in developing optimization methods for solving these problems. In this regard, researchers considered the optimization problem of the multi-speed gearbox using a multi-objective evolutionary algorithm involving four conflicting objectives [13]. The author applied non-dominated sorting genetic algorithm (NSGA-II) and showed the efficiency of proposed method in solving mixed integer optimization problem under several non-linear constraints in obtaining the optimal Pareto curve. Recently, authors considered the multi-objective optimization of two-stage helical gearbox, considering two conflicting objectives the volume of the gearbox and the total power loss [14]. Furthermore, besides the constraints such as bending stress and pitting stress, in this study authors employed tribological constraints such as scuffing and wear. Similarly, authors consider the two-stage spur gearbox which is formulated as the multi-objective optimization problem with three objectives, e.g., volume, power-output and center distance under several geometric and design constraints [9]. The obtained results concluded that the proposed method obtains better design solutions compared to those reported in the literature. Recently, authors considered the multi-objective optimization of the planetary gearbox with the same objectives [15]. In addition, several different constraints have been considered, such as regular strength requirement as well as the bearing selection and scuffing.

Generally, the MOO problems have been solved by converting them into the equivalent single-objective optimization problem using the conventional techniques, known as scalarization procedures, such as: weighted-sum and epsilon-constraint methods [16,17]. The widely applied weighted-sum method has good convergence properties; however, it can only be applied to obtain convex parts of the Pareto optimal front [18]. On the other hand, epsilon-constraint method can be applied not only to obtain convex parts of the Pareto set, but also for the non-convex parts [17]. However, there is no unified strategy in the literature to deal with epsilon-constraints values that always lead to feasible problem instances. To overcome these drawbacks, recently several parameter-free methods have been proposed in the literature that work directly with the multi-objective problem, without scalarization [19,20].

To successfully solve MOO problems, various optimization methods have been applied from conventional gradient-based methods, surrogate-based algorithms to nature-inspired evolutionary algorithms [20,21,22]. In recent years, researchers have applied different types of evolutionary algorithms (EAs) to solve multi-objective optimization problems [8,9,23]. However, since the conventional EAs are designed to be used for single-objective unconstrained optimization problems, they cannot be directly applied to solve constrained multi-objective optimization problems. To deal with the constraints, several different constraint handling techniques have been proposed in the literature, which can be broadly classified as: penalty functions, separate considerations of infeasible and feasible solutions and dealing with constraints using the MOO problem [24]. First method is the most widely employed technique to deal with constraints, where COP is transformed into the unconstrained optimization problem by introducing a penalty value term to the objective function, which is related to the constraints. The second method deals with the constraints in two stages. First, the algorithm optimizes constraints without taking into the account the objective functions. Then, in the second stage, when significant number of feasible solutions has been obtained the algorithms focuses on the optimization of objective functions. The newest technique, among the above-mentioned, treats the constraints of the MOO problem as the additional objectives, and therefore transforming the constrained optimization problem to unconstrained MOO problem.

In recent years, several metaheuristic algorithms inspired by nature, such as Genetic Algorithm (GA) [10], particle swarm optimization (PSO) [25], differential evolution (DE) [26], firefly algorithm (FA) [27], gray wolf optimization (GWO) [28] etc. have been applied to solve a wide array of engineering optimization problems. Among them the GA, as a promising optimization method from the group of EAs has been widely applied to solve different SOO problems regarding the gearbox optimization, such as weight minimization, power loss minimization etc., [3,4,5]. Besides the GA, several papers have been published that deal with SOO of gear trains using the algorithms, such as DE, PSO and GWO [28,29,30]. Regarding the MOO problems in gearbox design, most of the papers outlined previously are relying on the well-known NSGA-II and other improved variations of this algorithm to solve this complex optimization problem [10,13,14]. Recently, several papers appear in the literature that deal with the MOO problems using the hybrid variants of the well-known EAs [8,31]. Although a lot of work in the literature is devoted on developing more effective optimization algorithms and improving the performance of the existing ones, due to the no free lunch theorem of mathematical optimization the algorithm maximizing performance on one class of optimization problem will likely perform worse at others. The no free lunch theorem states that if no prior assumptions about the optimization problems can be made, averaging over all possible optimization problems, all optimization algorithms have equal performance [32], i.e., it is not possible to provide general-purpose universal optimization algorithm, rather only to provide a strategy that can outperform another on the specific group of optimization problems under consideration.

Among different EAs, DE emerged as the efficient algorithm which can successfully obtain global optimum of multimodal optimization problems. This algorithm has been initially proposed by Storn and Price [26] to solve real-valued unconstrained single-objective optimization problems. However, due to its simple structure and easy implementation it has since been expanded to solve mixed integer, constrained and multi-objective optimization problems [2,33]. The performance of DE algorithm depends in large measure on the appropriate choice of its main control parameters, such as scale factor F, crossover rate $CR$ and population size $N_P$ [34]. To address the issues with optimization performance, the researchers first considered the appropriate choice of the control parameter values for different type of objective function landscapes. However, it has been shown that single fixed control parameter value chosen for the certain problem is not adequate for the whole optimization process. Therefore, researchers developed several adaptive and self-adaptive strategies for changing the value of control parameters based on some information from the optimization process [34]. Furthermore, to successfully deal with MOO problems researchers proposed numerous variants of DE algorithm, which include external archive where the successfully non-dominated solutions found during the evolutionary process are stored and lately, during the evolution process, employed to direct the individuals towards global optimum [11,35].

Another widely applied algorithm from the group of EAs is the PSO algorithm. This method is inspired by the intelligence of the swarm of birds searching for food [25]. Initially, it has been proposed for solving unconstrained single-objective optimization problems in the continuous domain. In the last few decades, PSO algorithm has been applied to several practical engineering optimization problems, where it has shown good performance. However, the conventional PSO algorithm has some drawbacks in solving large scale complex optimization problems, such as premature convergence and trapping into the local optima [36]. As with other EAs, the optimization performance of PSO algorithm depends largely on the proper choice of the control parameters such as social $c_1$ and cognitive $c_2$ acceleration coefficients and population size $N_P$. In this way, to deal with above-mentioned drawbacks there has been several papers dealing with the appropriate choice of these parameters, among which the time-varying acceleration coefficients have been widely applied [36,37]. To improve the performance in solving the MOO problems researchers have introduced the archiving strategies into the PSO algorithm, where the best individuals are selected from a non-dominated external archive [11,38]. Furthermore, researchers studied the impact of different PSO topologies on the outcome of the optimization process [22].

To further enhance the optimization performance of PSO algorithm in solving MOO problems, in this paper an improved hybrid PSO and DE algorithm, called Multi-Objective Hybrid Particle Swarm Optimization Differential Evolution (MHPSODE) algorithm, has been proposed to deal with the complex MOO of planetary gearbox design. To address the issues that appear during the optimization of the conventional DE and PSO algorithms, in this paper a hybrid algorithm has been proposed which successfully combines the traits of each algorithm and overcomes the disadvantages. The hybridization has been performed by introducing the mutation operators DE/rand/1 and DE/current-to-best/1 from DE algorithm into the velocity update equation of the PSO algorithm, with the adaptive normalized population spacing parameter employed to select the appropriate mutation operator for the current optimization condition. The experimental results have demonstrated that the proposed MHPSODE algorithm is able to achieve better balance between convergence and diversity compared to well-known existing algorithms, such as NSGA-II and DEMO by adopting the proposed adaptive mutation mechanism.

In summary, the main contributions of this paper can be outlined as follows:

• The methodology for the MOO of the single stage planetary gearbox optimization, along with a comprehensive number of conflicting objective functions and several critical design constraints has been formulated in this paper.
• The theoretical formulation for the calculation of the planetary gearbox power efficiency has been outlined, as a significant objective in the formulated MOO problem. To evaluate the efficiency of the planetary gearbox and perform numerical analysis the appropriate numerical procedure has been presented.
• To improve the optimization performance, especially while solving complex multi-objective optimization problems, the hybrid algorithm between DE and PSO algorithms, called MHPSODE algorithm has been proposed. The hybridization between two well-known EAs has been achieved by introducing mutation operators from DE algorithm into the velocity update equation of the PSO algorithm. Based on an adaptive normalized population spacing parameter the proposed algorithm is able to choose appropriate modified mutation operator for the current optimization conditions.
• The performance of the proposed MHPSODE algorithm has been verified by comparing it with well-known multi-objective metaheuristic algorithms such as NSGA-II and DEMO on 10 benchmark MOO problems defined in CEC2009. The experimental results demonstrated that the proposed MHPSODE algorithm has significantly better optimization performance compared to other existing algorithms, in terms of both quality of the obtained Pareto solutions and convergence.

This paper is organized as follows. In Section 2 the formulation of the multi-objective planetary gear optimization problem has been given. Furthermore, in Section 2.1 the formulation of the constrains for the considered optimization problem have been outlined. Section 2.2 gives the mathematical formulation of the planetary gear train efficiency model. Section 3 presents the problem of constrained multimodal optimization. Furthermore, in Section 3.1 and Section 3.2 the theoretical formulation of the PSO and DE algorithms, has been outlined, respectively. In addition, Section 3.3 gives the formulation of the proposed hybrid HPSODE algorithm. Section 4 presents the results of the numerical simulation, where in Section 4.1 the focus is on the experiments dealing with planetary gearbox optimization, and in Section 4.2 the statistical comparison of the optimization performance of the proposed algorithm on the CEC2009 benchmark problems has been outlined. Finally, the conclusions are drawn in Section 5.

2. Problem Formulation

This paper considers the problem of MOO of the planetary gearbox with the aim to obtain the gear design parameters, which lead to construction with lower weight and volume while simultaneously having increased efficiency and service life. In this regard, Figure 1 illustrates the planetary gearbox considered for this MOO problem, which consists of a carrier, three planet gears (g) and a single sun gear with external gearing (a), input shaft, output shaft and a ring gear with internal gearing (b). The power comes thought the input shaft, on which the sun gear is mounted, then it is transferred by meshing engagement between gears to planet gears, which are mounted on the carrier, while the ring gear is kept fixed. Then the power is outputted thought the output shaft, which is directly connected to the carrier, as illustrated in Figure 1.

The materials, input power and speed, etc. and other parameters used in the design of the planetary gearbox in this study are outlined in

Table 1. Parameters of the planetary gearbox considered in the paper.

Parameters	Units	Symbol	Value
Input Power	[kW]	$P_a$	175
Input speed	[${\min }^{-1}$]	$n_a$	2750
Pressure angle	[${}^{\circ }$]	${\alpha }_n$	20
Gear material	18CrNi8
Gear surface Roughness	[$\mathsf{\mu }\mathrm{m}$]	$R_a$	0.8
Factor of safety against bending	[-]	${S_F}_{min}$	1.2
Factor of safety against pitting	[-]	${S_H}_{min}$	1.25
Number of planet gears	[-]	$n_w$	3

.

In this regard, in this paper the following objective functions have been formulated for the considered MOO optimization problem:

• Center distance

  $$f_1=a=\frac{m_n{z}_a}{cos\beta }\left(1+u_{a-g}^H\right)\frac{cos{\alpha }_t}{cos{\alpha }_{wt}},$$

  (1)

  where the involute of the pressure angle is defined as

  $$\mathrm{inv}{\alpha }_{wt}=2\frac{x_a+x_b}{z_a+z_b}tan{\alpha }_n+\mathrm{inv}{\alpha }_t.$$

  (2)
• Contact ratio

  $$f_2={\epsilon }_{\alpha \left(a-g\right)}=\frac{0.5\left(\sqrt{d_{a\left(a\right)}^2-d_{a\left(b\right)}^2}+\sqrt{d_{a\left(g\right)}^2-d_{a\left(g\right)}^2}\right)-asin{\alpha }_{wt}}{\pi m_{t}cos{\alpha }_t}.$$

  (3)
• Total volume of gears

  $$f_3=V\left(x\right)=\frac{\pi }{4}b\left[{d_{\left(a\right)}}^2+n_w\left(d_{a\left(b\right)}^2-D^2\right)+\left(d_{\left(g\right)}^2-d_s^2\right)\right].$$

  (4)
• Safety factor for bending stress

  $$f_4=S_{F\left(a\right)}\left(x\right)=\frac{{\left[{\sigma }_F\right]}_{M\left(a\right)}}{{\sigma }_{F\left(a\right)}}.$$

  (5)
• Safety factor for contact stress

  $$f_5=S_{H\left(a\right)}\left(x\right)=\frac{{\left[{\sigma }_H\right]}_{M\left(a\right)}}{{\sigma }_{H\left(a\right)}}.$$

  (6)
• Safety factor for bending stress for ring gear

  $$f_6=S_{F\left(b\right)}\left(x\right)=\frac{{\left[{\sigma }_F\right]}_{M\left(b\right)}}{{\sigma }_{F\left(b\right)}}.$$

  (7)

In the above equations, the tooth root stress for the sun gear is determined according to

$${\sigma }_{F,a}=\frac{F_t}{b{m}_n}{Y}_{F_a}{Y}_{S_a}{Y}_{\epsilon }{K}_A{K}_v{K}_{F\alpha },$$

(8)

while the critical root stress is obtained as

$${\left[{\sigma }_F\right]}_M={\sigma }_{Flim}{Y}_{ST}{Y}_{RrelT}.$$

(9)

In this regard, the safety factor against breakage can be obtained according to the equation

$$S_{F,a}=\frac{{\left[{\sigma }_F\right]}_{\left(a\right)}}{{\sigma }_{F\left(a\right)}}\ge S_{Fmin},$$

(10)

where the $S_{Fmin}$ is the minimum allowed safety factor value, which is given in Table 1.

Effective contact stress is determined according to the expression

$${\sigma }_H=Z_H{Z}_E{Z}_{\epsilon }\sqrt{\frac{F_t}{b{d}_{\left(a\right)}}\frac{u+1}{u}{K}_A{K}_v{K}_{H\alpha }},$$

(11)

where the contact stress for sun gear – planet gear pair is obtained as the minimum value, according to

$${\sigma }_{Ha}={\sigma }_{Hg}<{\left[{\sigma }_H\right]}_M=min\left\{{\left[{\sigma }_H\right]}_{\left(a\right)},{\left[{\sigma }_H\right]}_{\left(g\right)}\right\}.$$

(12)

The critical contact stress is determined as

$${\left[{\sigma }_H\right]}_M={\sigma }_{Hlim}{Z}_N{Z}_L{Z}_R,$$

(13)

while the safety factor form pitting is calculated as

$$S_{Ha,g}=\frac{{\left[{\sigma }_H\right]}_{\left(a,g\right)}}{{\sigma }_H}\ge S_{Hlim},$$

(14)

where the $S_{Hmin}$ is the minimum allowed safety factor value, which is given in Table 1. The different factors used in the above equations are given in the Appendix A.

2.1. Constraints Formulation

In this section, several design constraints regarding the strength requirements employed in the MOO problem considered in this paper have been outlined. Among others, the employed constraints include bending strength, pitting strength, assembly condition, etc. By implementing these constraints, the requirement for the gearbox with enhanced service life of the components is successfully satisfied.

2.1.1. Bending Constraints

The constraint regarding the gear tooth bending strength is defined as

$$g_{\left(1,2,3\right)}=\frac{{\left[{\sigma }_F\right]}_{M\left(a,g,b\right)}}{{\sigma }_{F\left(a,g,b\right)}}-S_F>0.$$

(15)

2.1.2. Pitting Constraints

The constraint against surface fatigue resistance has been defined as

$$\begin{array}{c}\hfill g_4=\frac{{\left[{\sigma }_H\right]}_{M\left(a,g\right)}}{{\sigma }_H}-S_H>0,\end{array}$$

(16)

$$\begin{array}{c}\hfill g_5=\frac{{\left[{\sigma }_H\right]}_{M\left(g,b\right)}}{{\sigma }_H}-S_H>0.\end{array}$$

(17)

2.1.3. Space Requirement

Mounting of the planet gears in the gearbox assembly demands the appropriate clearance between tip circles of gears in mesh to exist. Therefore, the following constraint needs to be satisfied

$$g_7=2asin\left(\frac{\pi }{n_w}\right)-f_z-d_{a-g}\ge 0.$$

(18)

where $f_z=0.5·m_n$ is the minimum clearance.

2.1.4. Assembly Condition

To avoid possible interference of teeth during the meshing process, the condition that must be satisfied is that simultaneous meshing of central sun gear with planet gears must be always satisfied. In this regard, the equality constraint is defined as

$$h_1=\frac{z_a{z}_b}{n_{w}D\left(z_g,z_b\right)}-i=0.$$

(19)

2.2. Planetary Gear Train Efficiency

The load dependent power losses arise during the operation of gears primarily as a consequence of relative motion between two gear tooth surfaces in contact under the normal load, where the sliding and rolling friction always exist. The power loss of the gear pair is expressed by means of efficiency, which is defined as a degree to which a system is successful in producing a desired result [10]. The efficiency of the gear pair along active length of the line of action (LOA) for a constant angular speed ${\omega }_1$ is obtained in the integration according to

$${\eta }_{gb}^H=\frac{1}{l_a}{\int }_A^D{\eta }_i\left(\xi \right)d\xi ,$$

(20)

in which the instantaneous efficiency is obtained as

$${\eta }_i\left(\xi \right)=\frac{T_2}{T_1}\frac{1}{u_{gb}^H},$$

(21)

where $T_1$ is the constant torque acting on the pinion $T_2$ is the torque acting on the driven gear, which needs to be determined, and $u_{gb}^H$ is the relative gear ratio, $l_a=\overline{A_2{E}_2}$ is the active length of the path of contact, $\xi$ is the coordinate along the LOA.

Based on the Coulomb’s Law of friction, the sliding friction forces acting on ith meshing tooth of the pinion and driven gear can be calculated as

$$\begin{array}{cc}\hfill & F_{\mu pi}={\mu }_i\left(\xi \right)F_{ni}\left(\xi \right)\mathrm{sign}\left(V_{si}\left(\xi \right)\right),i=1,2,\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & F_{\mu gi}=-{\mu }_i\left(\xi \right)F_{ni}\left(\xi \right)\mathrm{sign}\left(V_{si}\left(\xi \right)\right),i=1,2,\hfill \end{array}$$

(23)

where ${\mu }_i\left(\xi \right)$ is the time-varying friction coefficient of ith mesh tooth pair, $\mathrm{sgn}\left(·\right)$ is the signe function, and $F_{ni}$ denotes the normal force acting upon the two teeth in contact. According to the Figure 2, the direction of the sliding friction force should be opposite to that of the relative sliding velocity $V_{si}\left(t\right)$ between surfaces of two teeth in contact, and is always perpendicular to the LOA. In the process of the gear engagement, the tribological conditions change as a consequence of varying mesh properties and change in lubricant film thickness as the one surface of the tooth rolls over the other tooth. In this way, the friction coefficient varies with the angular position of each gear and the direction of friction force changes at the pitch point.

Therefore, Benedict and Kelley [39] used elasto-hydrodynamic lubrication (EHL) theory to explain the friction in gear teeth in contact, and based on experimental results proposed an equation for the prediction of the dynamic friction coefficient, which is given as follows

$${\mu }_i\left(\xi \right)=C_1{log}_{10}\left(\frac{C_2·F_{ni}\left(\xi \right)}{b}{\eta }_o{V}_{si}\left(\xi \right)V_{ei}^2\left(\xi \right)\right)i=1,2,$$

(24)

where b is the width of the gear, ${\eta }_0$ is the dynamic viscosity of the oil and the coefficients $C_1$ and $C_2$ are taken as $C_1=0.0127$ and $C_2=29.66$, respectively.

The intensity of the rolling friction is determined according to the expression

$$F_{Rpi}\left(\xi \right)=C·h_i\left(\xi \right)·b·\mathrm{sgn}\left(V_{ei}\right),i=1,2,$$

(25)

where the minimum film thickness is determined based on the expression given by Dowson and Higginson [40], as follows

$$h_i\left(\xi \right)=1.6·{{\alpha }_0}^{0.6}·{\left({\eta }_o·V_{ei}\right)}^{0.7}·E^{0.003}\frac{R^{0.43}}{F_n^{0.13}},$$

(26)

and ${\alpha }_0$ denotes viscosity-pressure coefficient of lubricant, R is the effective radius of curvature, and E is Young modulus of gear material.

Therefore, from the static equilibrium condition of forces acting upon gears in mesh, for $F_{n1}=F_{n2}$, the following expression for the normal force is obtained

$$F_n\left(\xi \right)=\frac{T_p+F_{Rp1}\left(\xi \right)\xi +\lambda \left(\xi \right)F_{Rp2}\left(\xi \right)l_{p1}\left(\xi \right)}{{\mu }_1\xi \mathrm{sign}\left(V_{s1}\left(\xi \right)\right)+2r_{bp}+\lambda \left(\xi \right){\mu }_2\left(\xi \right)l_{p1}\left(\xi \right)\mathrm{sign}\left(V_{s2}\left(\xi \right)\right)}.$$

(27)

The sliding and entraining velocities $V_{si}\left(\xi \right)$ and $V_{ei}\left(\xi \right)$, respectively, are determined according to the expression

$$\begin{array}{cc}\hfill & V_{si}\left(\xi \right)=v_{pi}\left(\xi \right)-v_{gi}\left(\xi \right)i=1,2,\hfill \\ \hfill & V_{ei}\left(\xi \right)=\frac{1}{2}\left[v_{pi}\left(\xi \right)+v_{gi}\left(\xi \right)\right]i=1,2,\hfill \end{array}$$

(28)

where the tangential velocities $v_{pi}\left(\xi \right)$ and $v_{gi}\left(\xi \right)$ of pinion and gear on any point of the profile of the gear are determined as follows

$$\begin{array}{cc}\hfill & v_{pi}\left(\xi \right)=l_{pi}\left(\xi \right)·{\omega }_{1}i=1,2,\hfill \\ \hfill & v_{gi}\left(\xi \right)=l_{gi}\left(\xi \right)·{\omega }_{2}i=1,2,\hfill \end{array}$$

(29)

Here ${\omega }_1$ and ${\omega }_2$ denote the rotational speed of the pinion and gear, respectively, and the $l_{pi}\left(\xi \right)$ and $l_{gi}\left(\xi \right)$ are the lengths given by

$$\begin{array}{cc}\hfill & l_{p1}\left(\xi \right)=\xi +p_b,\hfill \\ \hfill & l_{p2}\left(\xi \right)=asin{\alpha }_w+l_{p1}\left(\xi \right),\hfill \\ \hfill & l_{p3}\left(\xi \right)=asin{\alpha }_w+\xi ,\hfill \\ \hfill & l_{gi}\left(\xi \right)=l_a-l_{pi}\left(\xi \right),i=1,2,\hfill \end{array}$$

(30)

where the active length of LOA is obtained as

$$l_a=\left(r_{b1}+r_{b2}\right)tan{\alpha }_w.$$

(31)

Based on the obtained total mesh force, the torque $T_2$ is than calculated according to the equation

$$T_2\left(\xi \right)=2F_n{r}_{bp}+l_{g2}\left(F_{\mu p2}\left(\xi \right)+F_{Rp2}\right)+\lambda \left(\xi \right)l_{g1}\left(F_{\mu p1}\left(\xi \right)+F_{Rp1}\right).$$

(32)

In the above expressions, $\lambda \left(\xi \right)$ denotes dimensionless coefficient, which takes into the account changes caused by the single-tooth meshing and double-tooth meshing alternation during the engagements, according to the expression

$$\lambda \left(\xi \right)=\left\{\begin{array}{cc}0,& l_a-p_b\le mod\left(\xi ,p_p\right)\le p_b\\ 1,& \mathrm{otherwise}\end{array}\right..$$

(33)

However, the intensity of the normal force between two surfaces cannot be obtained directly from Equation (27), since the values of the coefficients ${\mu }_i\left(\xi \right)$ and $h_i\left(\xi \right)$ are not a priori known. Therefore, the iterative procedure is established, where in the first iteration it is assumed that the value of the normal force is

$$F_n^{\left(0\right)}=\frac{T_p}{2r_{bp}}.$$

(34)

This iterative procedure is repeated until the termination criterion

$$\left|\frac{F_n^{\left(i+1\right)}-F_n^{\left(i\right)}}{F_n^{\left(i\right)}}\right|\le \epsilon ,$$

(35)

is satisfied. Then the torque acting on the driven gear and instantaneous efficiency can be calculated according to the expressions (32) and (21), respectively. In this way, the pseudo-code of the proposed iterative procedure is presented in Algorithm 1.

Algorithm 1 Computational methodology of efficiency prediction model

Initialize parameters such as ${\eta }_o$, ${\alpha }_o$, $F_n^{\left(0\right)}$,$\epsilon$

while  $\xi \le l_a$  do

while  $\left|\frac{F_n^{\left(i+1\right)}-F_n^{\left(i\right)}}{F_n^{\left(i\right)}}\right| \le \epsilon$  do

Calculate coefficients ${\mu }_i\left(\xi \right)$ and $h_i\left(\xi \right)$ according to Equations (24) and (26), respectively

Determine normal force $F_n^{\left(i\right)}$ using Equation (27)

$i\leftarrow i+1$

end while

Calculate torque on the driven gear $T_2\left(\xi \right)$

Determine instantaneous efficiency using

${\eta }_i\left(\xi \right)=\frac{T_2}{T_1}\frac{1}{u_{gb}^H}$

end while

Calculate total efficiency

${\eta }_{gb}^H=\frac{1}{l_a}{\int }_A^D{\eta }_i\left(\xi \right)d\xi$

3. Multi-Objective Optimization

Generally, the MOO problem can be defined as an optimization problem with $m=1,\dots ,M$ objective functions that are simultaneously minimized. Most of the practical engineering optimization problems have some limitations regarding the range of the decision variables, and therefore, the optimization problems usually are further expanded with several type of constraints to ensure the physical feasibility of the solution. These problems, are referred to as constrained multi-objective optimization problems (CMOPs), which can be formulated in the mathematical form as follows

$$\begin{array}{cc}\hfill & min{f}_m\left({\mathbf{x}}_i\right)m=1,\dots ,M\wedge i=1,\dots N_P\hfill \\ \hfill & s.t.\hfill \\ \hfill & g_k\left({\mathbf{x}}_i\right)\ge 0,k=1,\dots ,K\hfill \\ \hfill & h_l\left({\mathbf{x}}_i\right)=0,l=1,\dots ,L\hfill \\ \hfill & x_{i,j}\in \left[x_j^{Lower},x_j^{Upper}\right],j=1,\dots ,n\hfill \end{array}$$

(36)

where M denotes the total number of objective functions to be optimized, $g_k\left({\mathbf{x}}_i\right)$ is k-th inequality constraint, K is the total number of inequality constraints, $h_l\left({\mathbf{x}}_i\right)$ is lth equality constraint, L is the total number of equality constraints, $x_j^{Lower}$ and $x_j^{Upper}$ denote the lower and upper bounds of the decision variables in the search space, respectively. The set $\mathbb{F}$, which satisfies all the constraints, is called feasible domain $\mathbb{F}=\left\{\left.\mathbf{x}\right|g_k\left({\mathbf{x}}_i\right)\ge 0\wedge h_l\left({\mathbf{x}}_i\right)=0\wedge x_{i,j}\in \left[x_j^{Lower},x_j^{Upper}\right]\right\}$.

Definition 1.

Let $\mathbf{x}$ be the point in the feasible solution space, $\mathbf{x}\in \mathbb{F}$. The point $\mathbf{x}$ is said to be Pareto optimal if and only if in the feasible domain there does not exist another point $\mathbf{y}\in \mathbb{F}$ for which $f\left(\mathbf{y}\right)\le f\left(\mathbf{x}\right)$, $f\left(\mathbf{y}\right)\ne f\left(\mathbf{x}\right)$.

Definition 2.

Dominated solution: The solution $\mathbf{x}$ dominates the solution $\mathbf{y}$ (denoted as $\mathbf{x}\prec \mathbf{y}$) if and only if $f_i\left(\mathbf{x}\right)\le f_i\left(\mathbf{y}\right)$$\forall i\in \left\{1,2,\dots ,m\right\}$ and $\exists j\in \left\{1,2,\dots ,m\right\}$ such that $f_j\left(\mathbf{x}\right)<f_j\left(\mathbf{y}\right)$.

Definition 3.

Non-dominated solution: The solution $\mathbf{x}$ and solution $\mathbf{y}$ are non-dominate with each other if and only if $\exists j\in \left\{1,2,\dots ,m\right\}$ for which $f_j\left(\mathbf{x}\right)<f_j\left(\mathbf{y}\right)$ and $\nexists i\in \left\{1,2,\dots ,m\right\}$ such that $f_i\left(\mathbf{x}\right)<f_i\left(\mathbf{y}\right)$.

Definition 4.

Pareto front: The set, which is composed of non-dominated Pareto optimal solutions, is called Pareto front (PF), denoted as $X=\left\{\left.\mathbf{x}\in \mathbb{F}\right|\nexists \mathbf{y}\in \mathbb{F}\right.$ such that $f\left(\mathbf{y}\right)>f\left(\mathbf{x}\right).$

Definition 5.

Ideal point: the point ${\mathbf{x}}_{ideal}$ is the point, which has the best objective function values in Pareto front ${\mathbf{x}}_{ideal}=min\mathbf{f}\left(\mathbf{x}\right)=min\left\{f_1\left(\mathbf{x}\right),\dots ,f_i\left(\mathbf{x}\right),\dots ,f_m\left(\mathbf{x}\right)\right\}$.

3.1. Multi-Objective Particle Swarm Optimization Algorithm

The PSO is the optimization method form the group of EAs, which is inspired by the behavior of birds in the swarm [25]. The PSO algorithm consists of $N_P$ particles in the swarm, which can be represented with the set $P^{\left(G\right)}=\left\{\left.{\mathbf{x}}_i\right|{\mathbf{x}}_i\in \mathbb{F}\right\},i\in \left\{1,2,\dots ,N_P\right\}$, where G denotes the generation index. Each particle in the swarm represents a candidate solution to the considered optimization problem. Two parameters, namely current position ${\mathbf{x}}_i^{\left(G\right)}$ and velocity ${\mathbf{s}}_i^{\left(G\right)}$ are associated with each particle. The position of the particle is influenced by the best position previously visited by the particle, denoted as ${\mathbf{x}}_{pbest,i}^{\left(G\right)}$ and the best position discovered by the entire swarm of particles ${\mathbf{x}}_{gbest}^{\left(G\right)}$. In this regard, at each iteration of the PSO algorithm the particle updates its velocity and then the position according to the expression

$$\begin{array}{cc}\hfill {\mathbf{s}}_i^{\left(G+1\right)}& ={\mathbf{s}}_i^{\left(G\right)}+c_{1}ran{d}_1\left({\mathbf{x}}_{pbest,i}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)+c_{2}ran{d}_2\left({\mathbf{x}}_{gbest}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right),\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill {\mathbf{x}}_i^{\left(G+1\right)}& ={\mathbf{x}}_i^{\left(G\right)}+{\mathbf{s}}_i^{\left(G+1\right)},\hfill \end{array}$$

(38)

where $c_1$ denotes the cognitive coefficient, $c_2$ is social coefficient, $ran{d}_1$ and $ran{d}_2$ denote distinct random numbers in the range $ran{d}_i\in \left(0,1\right),i=1,2$. For solving the MOO problems, the procedure of the PSO algorithm need to be modified, in terms of mechanism for updating the personal best and global best positions [41]. The personal best position represents the best solution found by individual particle in the previous generations. For the MOO problem the personal best position for each particle is changed if the current personal best solution is dominated by the current position, which is written with the following equation

$$\begin{array}{c}\hfill {\mathbf{x}}_{pbest,i}^{\left(G+1\right)}=\left\{\begin{array}{c}{\mathbf{x}}_{pbest,i}^{\left(G\right)},\mathrm{if}{\mathbf{x}}_{pbest,i}^{\left(G\right)}\prec {\mathbf{x}}_i^{\left(G\right)}\\ {\mathbf{x}}_i^{\left(G\right)},\mathrm{otherwise}\end{array}\right..\end{array}$$

(39)

On the other hand, the global best position ${\mathbf{x}}_{gbest}^{\left(G\right)}$ represents the best solution found in the entire swarm up to current iteration. Unlike with SOO problems, in MOO problems due to the contradictory nature of the multiple objectives it is impossible to choose single global best solution for the ${\mathbf{x}}_{gbest}^{\left(G\right)}$. Therefore, an external archive is introduced into the PSO algorithm, which captures a set of non-dominated solutions obtained in the optimization process. In each iteration, to evaluate the diversity of the particles in the external archive the crowding distance is calculated, and then the ${\mathbf{x}}_{gbest}^{\left(G\right)}$ is selected by applying the binary tournament [22,41].

3.2. Differential Evolution Algorithm

DE algorithm is an effective and robust optimization algorithm, initially proposed to solve unconstrained optimization problems in the continuous domain [26]. The optimization process of the DE algorithm consists of the following phases: initialization, mutation, crossover and selection.

The DE algorithms begins the optimization process with a population of $N_p$ individuals $\left\{\left.{\mathbf{x}}_i\right|x_{i,j}^{\left(G\right)}\in \left[x_j^{Upper},x_j^{Lower}\right]\right\}$$\forall i\in \left\{1,2,\dots ,N_P\right\}$, uniformly and randomly deployed over the feasible solution space, where each individual is represented by an n-dimensional vector ${\mathbf{x}}_i={\left[x_{i,1},\dots ,x_{i,j},\dots ,x_{i,n}\right]}^T$. Therefore, each individual ${\mathbf{x}}_i$ represents a potential solution to the optimization problem. Every component of the ith individual should be generated within the upper and lower bounds $x_{i,j}\in \left[x_j^{Upper},x_j^{Lower}\right]$. In this regard, the component of the ith vector is generated according to the following expression

$$\begin{array}{c}{x}_{i,j}^{\left(G\right)}=x_j^{Lower}+\left(x_j^{Upper}-x_j^{Lower}\right)\mathrm{rand}\left(0,1\right)\\ i=1,2,\dots ,N_P\wedge j=1,2,\dots ,n,\end{array}$$

(40)

where $rand\left(0,1\right)$ is a uniformly and randomly generated number in the range $\left(0,1\right)$.

The mutation operator in DE algorithm is considered to be a random disturbance term, which allows the algorithm to explore the search space and escape from local optima, and in this way keeping the diversity of the population. In the literature, exist several different mutation operators, which have different characteristic and can be suitable for certain type of the optimization problem or particular phase of the optimization process. In this way, the most commonly applied mutation operators include [34]:

• DE/rand/1

$${\mathbf{x}}_i^{\left(G\right)}={\mathbf{x}}_{r_1}^{\left(G\right)}+F\left({\mathbf{x}}_{r_2}^{\left(G\right)}-{\mathbf{x}}_{r_3}^{\left(G\right)}\right),$$

(41)

• DE/rand/2

$${\mathbf{x}}_i^{\left(G\right)}={\mathbf{x}}_{r_1}^{\left(G\right)}+F\left({\mathbf{x}}_{r_2}^{\left(G\right)}-{\mathbf{x}}_{r_3}^{\left(G\right)}\right)+F\left({\mathbf{x}}_{r_4}^{\left(G\right)}-{\mathbf{x}}_{r_5}^{\left(G\right)}\right),$$

(42)

• DE/best/1

$${\mathbf{x}}_i^{\left(G\right)}={\mathbf{x}}_{best}^{\left(G\right)}+F\left({\mathbf{x}}_{r_1}^{\left(G\right)}-{\mathbf{x}}_{r_2}^{\left(G\right)}\right),$$

(43)

• DE/current-to-best/1

$${\mathbf{x}}_i^{\left(G\right)}={\mathbf{x}}_i^{\left(G\right)}+F\left({\mathbf{x}}_{best}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)+F\left({\mathbf{x}}_{r_1}^{\left(G\right)}-{\mathbf{x}}_{r_2}^{\left(G\right)}\right),$$

(44)

where $r_1,r_2,r_3,r_4$ and $r_5$ represent integers from the set $\left\{1,2,\dots ,N_P\right\}\setminus i$ for which $r_1\ne r_2\ne r_3\ne r_4\ne r_5$, ${\mathbf{x}}_{best}$ denotes the individual with the minimum objective function value, and F represents the constant scaling factor $F\in \left[0,2\right]$.

After the mutation, in the DE algorithm the crossover operator is applied to increase diversity. According to the binomial crossover the trial vector ${\mathbf{u}}_i^{\left(G\right)}={\left[u_{i,1}^{\left(G\right)},\dots ,u_{i,n}^{\left(G\right)}\right]}^T$ is generating by mixing the jth component of the target $x_{i,j}^{\left(G\right)}$ and the corresponding mutant vector $u_{i,j}^{\left(G\right)}$, according to the expression

$$\begin{array}{cc}\hfill & u_{i,j}^{\left(G\right)}=\left\{\begin{array}{cc}{v}_{i,j}^{\left(G\right)},& \mathrm{if}ran{d}_j\le CR\vee j=j_{rand}\\ x_{i,j}^{\left(G\right)},& \mathrm{otherwise}\end{array}\right.,\hfill \\ \hfill & j=1,\dots ,n,\hfill \end{array}$$

(45)

where $ran{d}_j$ is the random number in the range $\left(0,1\right)$, index $j_{rand}$ is a randomly chosen integer from the set $\left\{1,2,\dots ,n\right\}$ and $CR$ is the crossover probability such that $CR\in \left[0,1\right]$.

Finally, after crossover the selection operator is applied to determine whether the individuals generated by the crossover and mutation operators are better than individuals in the previous generation. Therefore, based on the objective function value of the trial vector and objective function value of the target vector the selection is performed according to the expression

$${\mathbf{x}}_i^{\left(G+1\right)}=\left\{\begin{array}{cc}{\mathbf{u}}_i^{\left(G\right)},& \mathrm{if}f\left({\mathbf{u}}_i^{\left(G\right)}\right)\le f\left({\mathbf{x}}_i^{\left(G\right)}\right)\\ {\mathbf{x}}_i^{\left(G\right)},& \mathrm{otherwise}\end{array}\right.,$$

(46)

where ${\mathbf{x}}_j^{\left(G+1\right)}$ is the individual kept for the next generation.

3.3. Multi-Objective Hybrid Particle Swarm Optimization and Differential Evolution Algorithm

To successfully solve MOO problems, the EAs must effectively balance the convergence and diversity of the swarm during the optimization process. It has been previously shown that the PSO algorithm tends to be easily trapped into the local optima and to lose the diversity [36]. Therefore, to keep an effective balance between convergence and diversity the mutation operators from the DE algorithm have been introduced into the optimization procedure of the PSO algorithm. In the literature, it has been well established that the DE/rand1 and DE/rand/2 mutation operators have good exploration ability and hence, are able to successfully discover the region of global optimum [33]. On the other hand, mutation operators DE/best/1 and DE/current-to-best/1 focus their search in the region of the best solution, and therefore have better exploitation and faster convergence towards optimal solution [33].

In this regard, to improve the global exploration ability and prevent premature convergence to one of the local optima, the key elements from the DE/rand/1 mutation operator are included into the velocity update equation of the PSO algorithm. Therefore, the term $c_1\left({\mathbf{x}}_{pbest,i}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)$ in velocity equation has been substituted with the term $F\left({\mathbf{x}}_{r_2}^{\left(G\right)}-{\mathbf{x}}_{r_3}^{\left(G\right)}\right)$, and the velocity update equation has the following form

$${\mathbf{s}}_i^{\left(G+1\right)}={\mathbf{s}}_i^{\left(G\right)}+ran{d}_{1}F\left({\mathbf{x}}_{r_2}^{\left(G\right)}-{\mathbf{x}}_{r_3}^{\left(G\right)}\right)+c_{2}ran{d}_2\left({\mathbf{x}}_{gbest}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right).$$

(47)

On the other hand, to improve the exploitation ability the term $F\left({\mathbf{x}}_{best}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)$ from the DE/current-to-best/1 algorithm has been introduced into Equation (37) instead of the term $c_{2}ran{d}_2\left({\mathbf{x}}_{gbest}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)$. Furthermore, to maintain the required level of diversity, the term $F\left({\mathbf{x}}_{r_1}^{\left(G\right)}-{\mathbf{x}}_{r_2}^{\left(G\right)}\right)$ has been added to Equation (37), which becomes

$$\begin{array}{cc}\hfill & {\mathbf{s}}_i^{\left(G+1\right)}={\mathbf{s}}_i^{\left(G\right)}+c_{1}ran{d}_1\left({\mathbf{x}}_{pbest,i}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)+c_{2}F\left({\mathbf{x}}_{best}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)+F\left({\mathbf{x}}_{r_1}^{\left(G\right)}-{\mathbf{x}}_{r_2}^{\left(G\right)}\right).\hfill \end{array}$$

(48)

Based on the above, to obtain a suitable diversity and convergence, which can balance the global exploration and local exploitation abilities of the particles, the normalized population spacing $P{S}^{\left(G+1\right)}$ has been introduced into the algorithm. The value of the normalized population spacing of the ith particle in the population in the $\left(G+1\right)$ th iteration can be determined according to the expression

$$P{S}^{\left(G+1\right)}=\sqrt{\frac{1}{N_P-1}\sum _{i=1}^{N_P}{\left({\overline{d}}^{\left(G+1\right)}-d_i^{\left(G+1\right)}\right)}^2},$$

(49)

where $d_i^{\left(G+1\right)}$ denotes the minimum L1-norm distance between the ith particle and all other particles in the population, and ${\overline{d}}^{\left(G+1\right)}$ represents the average minimum L1-norm distance of all particles.

Therefore, employing the introduced normalized population spacing $P{S}^{\left(G+1\right)}$ to combine the developed Equations (47) and (48), the expression for velocity update of the PSO algorithm can be stated with the pseudo-code given in Algorithm 2.

Algorithm 2 Presudocode of the velocity update equation of the hybrid algorithm

if   $P{S}^{(G+1)}>0.75$  then

if  $rand>0.5$  then

${\mathbf{s}}_i^{\left(G+1\right)}={\mathbf{s}}_i^{\left(G\right)}+ran{d}_{1}F\left({\mathbf{x}}_{r_2}^{\left(G\right)}-{\mathbf{x}}_{r_3}^{\left(G\right)}\right)+c_{2}ran{d}_2\left({\mathbf{x}}_{gbest}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)$

else

${\mathbf{s}}_i^{\left(G+1\right)}={\mathbf{s}}_i^{\left(G\right)}+c_{1}ran{d}_1\left({\mathbf{x}}_{pbest,i}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)+c_{2}ran{d}_2\left({\mathbf{x}}_{gbest}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)$

end if

else if $P{S}^{(G+1)}\le 0.75$  then

if  $rand>0.5$  then

$\begin{array}{c}{\mathbf{s}}_i^{\left(G+1\right)}={\mathbf{s}}_i^{\left(G\right)}+c_{1}ran{d}_1\left({\mathbf{x}}_{pbest,i}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)+c_{2}F\left({\mathbf{x}}_{best}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)+F\left({\mathbf{x}}_{r_1}^{\left(G\right)}-{\mathbf{x}}_{r_2}^{\left(G\right)}\right)\end{array}$

else

${\mathbf{s}}_i^{\left(G+1\right)}={\mathbf{s}}_i^{\left(G\right)}+c_{1}ran{d}_1\left({\mathbf{x}}_{pbest,i}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)+c_{2}ran{d}_2\left({\mathbf{x}}_{gbest}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)$

end if

end if

According to the pseudo-code given in Algorithm 2, it can be observed that based on the value of the population spacing $P{S}^{\left(G+1\right)}$ the appropriate equation for the velocity update will be chosen. Therefore, for $P{S}^{\left(G+1\right)}>0.75$ the diversification of the population is uneven, and the exploitation ability should be enhanced. Therefore, the improved expression in Equation (47) will be applied with the probability of 0.5, otherwise Equation (37) is applied. On the other hand, smaller value of $P{S}^{\left(G+1\right)}$, e.g., $P{S}^{\left(G+1\right)}\le 0.75$, indicates that the local exploitation should be improved. Therefore, the expression (48) is applied with the probability 0.5, otherwise Equation (37) is applied.

To further increase diversity of the population the crossover operator from the DE algorithm has been incorporated into the PSO algorithm. Based on the above, the pseudo-code of the proposed MHPSODE iterative procedure is presented in Algorithm 3.

Algorithm 3 Pseudo-code of the proposed hybrid particle swarm differential evolution algorithm

Initialize the parameters $N_P$, n, $MaxIter$

Uniformly and randomly generate positions of the particles in the swarm

Initialize velocities of each particle in the swarm ${\mathbf{s}}_i^{\left(0\right)},\forall i\in 1,2,\dots ,N_P$

Set the initial ${\mathbf{x}}_{i,pbest}^{\left(0\right)}$ and ${\mathbf{x}}_{gbest}^{\left(0\right)}$ values

while $iter<MaxIter$ do

update the inertia weight $w^{\left(G\right)}$

for i = 1:$N_P$ do

Calculate the velocity

if $P{S}^{(G+1)}>0.5$ then

if $rand>0.5$ then

${\mathbf{s}}_i^{\left(G+1\right)}={\mathbf{s}}_i^{\left(G\right)}+ran{d}_{1}F\left({\mathbf{x}}_{r_2}^{\left(G\right)}-{\mathbf{x}}_{r_3}^{\left(G\right)}\right)+c_{2}ran{d}_2\left({\mathbf{x}}_{gbest}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)$

else

${\mathbf{s}}_i^{\left(G+1\right)}={\mathbf{s}}_i^{\left(G\right)}+c_{1}ran{d}_1\left({\mathbf{x}}_{pbest,i}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)+c_{2}ran{d}_2\left({\mathbf{x}}_{gbest}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)$

end if

else if $P{S}^{(G+1)}\le 0.5$ then

if $rand>0.5$ then

$\begin{array}{c}{\mathbf{s}}_i^{\left(G+1\right)}={\mathbf{s}}_i^{\left(G\right)}+c_{1}ran{d}_1\left({\mathbf{x}}_{pbest,i}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)+c_{2}F\left({\mathbf{x}}_{best}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)+F\left({\mathbf{x}}_{r_1}^{\left(G\right)}-{\mathbf{x}}_{r_2}^{\left(G\right)}\right)\end{array}$

else

${\mathbf{s}}_i^{\left(G+1\right)}={\mathbf{s}}_i^{\left(G\right)}+c_{1}ran{d}_1\left({\mathbf{x}}_{pbest,i}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)+c_{2}ran{d}_2\left({\mathbf{x}}_{gbest}^{\left(G\right)}-{\mathbf{x}}_i^{\left(G\right)}\right)$

end if

end if

Update the position of particles

${\mathbf{x}}_i^{\left(G+1\right)}={\mathbf{x}}_i^{\left(G\right)}+{\mathbf{s}}_i^{\left(G+1\right)}$

Evaluate crossover operator from the DE algorithm according to Equation (45)

Determine the personal best solution according to

${\mathbf{x}}_{pbest,i}^{\left(G+1\right)}=\left\{\begin{array}{c}{\mathbf{x}}_{pbest,i}^{\left(G\right)},\mathrm{if}{\mathbf{x}}_{pbest,i}^{\left(G\right)}\prec {\mathbf{x}}_i^{\left(G\right)}{\mathbf{x}}_i^{\left(G\right)},\mathrm{otherwise}\end{array}\right.$

Calculate the global best solution

end for

end while

4. Experimental Results

In this section, the experimental analysis using the numerical simulations has been caried out to determine the Pareto optimal solutions of the considered MOO problem of the planetary gearbox optimization and to analyze the improvements in terms of optimization performance of the proposed MHPSODE algorithm. In this regard, to analyze the optimization performance and to verify the performance of the modifications introduced in the proposed MHPSODE algorithm, in this paper, the statistical comparison has been performed between the proposed algorithm and several well-known algorithms, such as NSGA-II [42], MO_Ring_PSO_SCD [22] and DEMO [43] on CEC2009 benchmark problems.

4.1. The Planetary Gear Train Optimization

The results of the multi-objective planetary gearbox optimization problem, formulated in Section 2, have been presented in this subsection. In this regard, the minimization of the volume, center distance, safety factor for bending stress and contact stress, and minimization of power losses have been treated as the objectives. Results obtained from the formulated MOO problem have been compared with the reference gearbox outlined in the ISO/AGMA standards [44]. To perform comparison, the ISO VG 220 oil has been used in the numerical simulations, as referenced in the standard for the industrial gearbox.

Numerical results, such as the corresponding number of teeth on planet, central and ring gear, module as well as the profile shift coefficients obtained from the numerical simulation using the proposed MHPSODE and NSGA-II algorithms have been shown in

Table 2. The obtained values of design variables for the considered MOO problem of planetary gearbox and parameters of the reference industrial gearbox.

Algorithm	${\mathit{z}}_{\mathit{a}}$	${\mathit{z}}_{\mathit{g}}$	${\mathit{z}}_{\mathit{b}}$	${\mathit{m}}_{\mathit{n}}$[mm]	${\mathit{x}}_{\mathit{a}}$	${\mathit{x}}_{\mathit{g}}$	${\mathit{x}}_{\mathit{b}}$	${\mathit{\eta }}_{\mathbf{gb}}^{\mathit{H}}$	V[mm3]	a[mm]
Multi-objective optimization with gearbox center distance and efficiency as objectives
MHPSODE	22	44	45	6	0.25	0.251	0.25	0.9908	$2\times {10}^9$	150.5
NSGA-II	21	52	56	6	0.25	0.25	0.2545	0.9887	$2.2\times {10}^9$	221
Multi-objective optimization with gearbox center distance and bending stress
MHPSODE	22	32	102	6	0.25	0.25	0.4147	0.9735	$2.99\times {10}^9$	140
NSGA-II	22	37	105	6	0.25	0.25	0.274	0.9717	$3.98\times {10}^9$	180
Multi-objective optimization with gearbox center distance and contact stress
MHPSODE	18	22	45	10	0.5571	0.25	0.25	0.969	$3\times {10}^9$	138
NSGA-II	18	22	50	12	0.2848	0.25	0.494	0.965	$3.86\times {10}^9$	245
Multi-objective optimization with gearbox volume and efficiency as objectives
MHPSODE	17	22	45	10	0.25	0.25	0.41	0.9912	$2.45\times {10}^9$	199
NSGA-II	18	23	45	12	0.25	0.25	0.25	0.9689	$3.1\times {10}^9$	251
Industrial gearbox
-	20	37	94	9	-	-	-	0.97	$3.21\times {10}^9$	285

, where the values for the reference industrial gearbox are also presented for comparison. Due to the conflicting nature of multiple objectives, the result of the MOO problem is not a single unique solution, rather a set of solutions, e.g., a Pareto set. To extract a single solution, which can be compared to the reference values, the ideal solution must be determined. The ideal solution ${\mathbf{x}}_{ideal}$ represents the minimum best value for each individual objective, regardless of compliance with the rest of objective functions. Seeing that the ideal solution does not belong to the Pareto set, it is obtained as the approximation of the true ideal value from the obtained objective values on Pareto frontier. Therefore, the solutions, named compromise solution, of the considered MOO problem, presented in Table 2 are obtained as the results on the Pareto curve cosset to the ideal solution, where the measure for closeness is the Euclidean distance.

The Pareto optimal curve for the MOO problem obtained with the proposed MHPSODE and well-known NSGA-II algorithms, where the objectives are the power loss of the planetary gear set and center distance between sun and planet gears, is depicted in Figure 3.

From the Figure 3 it is observed that efficiency of the planetary gearbox and the center distance between sun and planet gears are conflicting objectives. The ideal solution for the center distance and power efficiency are 91 mm and 0.995%, respectively. Therefore, the compromise solution obtained by MHPSODE algorithm (marked with red diamond), as the solution which belongs to the Pareto set and is closest to the ideal solution is 150.5 mm for the center distance and 0.9908% for the gearbox efficiency.

Next, the center distance between sun and planet gears and the bending stress are considered to be the objectives of the MOO problem. Therefore, the corresponding Pareto curves obtained with the proposed and NSGA-II algorithms are depicted in in Figure 4.

As in the previous case, from the Figure 4 it is observed that considered objectives are conflicting. The ideal solution is 89.92 mm for center distance and 367 $\mathrm{N}/{\mathrm{mm}}^2$ for bending stress. Therefore, the compromise solution for this case is center distance 140 mm and bending stress 996.8 $\mathrm{N}/{\mathrm{mm}}^2$, obtained from the MHPSODE algorithm.

The obtained Pareto curves for the case when the considered objectives are the center distance and contact stress ${\sigma }_H$ are depicted in Figure 5.

Since, the considered objectives are conflicting, the ideal solution is determined where center distance is 40 mm and ${\sigma }_H$ is 615 $\mathrm{N}/{\mathrm{mm}}^2$. Therefore, we are adopting center distance of 138 mm and contact stress of 1070 $\mathrm{N}/{\mathrm{mm}}^2$ as the compromise solution obtained by MHPSODE algorithm.

Furthermore, the minimization of volume of the planetary gearbox, and therefore its mass, is of great importance in applications where the lightweight constructions are required, such as in aerospace industry. In this regard, the stated MOO problem has been carried out, where the considered objectives are the volume of the gears in the gearbox and the power loss of the planetary gearbox. Figure 6 shows the Pareto curves for this case obtained from proposed MHPSODE and NSGA-II algorithms.

From the results depicted in Figure 6, it can be observed that the considered objective functions are conflicting. Analyzing the obtained Pareto set, starting from the left to right and moving along the Pareto set, it can be observed that with the increase in the value of gearbox volume, leads to simultaneous decrease in the value of the planetary gearbox efficiency, and vice versa. Therefore, the ideal solution is volume of gears $1.28\times {10}^9$${\mathrm{mm}}^3$ and gearbox efficiency of 0.995%. The point on the Pareto front closest to the ideal solution, a compromise solution, is volume of gears $2.45\times {10}^9$${\mathrm{mm}}^3$ and efficiency of 0.9912% obtained form MHPSODE algorithm.

Analyzing the results shown in Table 2, and taking into the consideration the results presented in Figure 3, Figure 4, Figure 5 and Figure 6, it can be observed that the best solution is obtained using the proposed MHPSODE algorithm. The best solutions in terms of the gearbox efficiency, among outlined is obtained for the case when the gears volume and efficiency is considered. Compared to the industrial gearbox reference, it leads to the 23.6% improvement in gearbox volume, 30% decrease in center distance, which directly correlate with gearbox size, and 2.2% improvement in gearbox efficiency. In addition, it can be observed that the best reduction in center distance is obtained in the case when the contact stress is considered to be one of the objectives. Furthermore, from Figure 3, Figure 4, Figure 5 and Figure 6 it can be observed that the NSGA-II algorithm cannot obtain the complete Pareto front, which is especially true for non-convex Pareto fronts in Figure 3 and Figure 4. Comparing the best solutions in Table 2 obtained by MHPSODE and NSGA-II algorithms, it can be observed the improvement of 30.4% in center distance, 10% improvement of gearbox volume and 0.25% improvement of gearbox efficiency for the solutions obtained by proposed MHPSODE algorithm. The availability of the full Pareto front is of great importance for the designer, as it provides a range of solutions from which to choose the most appropriate, and gives an insight into how the change of one objective function affects the other one.

Furthermore, in Figure 7 the MOO problem was considered where center distance between sun and planet gear and volume have been taken as the objectives.

From the results in Figure 7, it can be concluded that between two objectives exists linear correlation with the correlation coefficient 0.96. Therefore, we conclude that when considering MOO of gears only one of the objectives should be considered. This can increase the efficiency and reduce the workload on the design process.

Furthermore, in this sub case the MOO has been performed taking into consideration 5 different ISO grade oils in the numerical simulations. The analysis has been carried out volume and efficiency as the objective functions to determine the best performing oil for the formulated problem. Therefore, Figure 8 shows the obtained Pareto fronts, where four oil types have been employed e.g., ISO VG 100, ISO VG 150, ISO VG 680 and ISO VG 460.

According to the results in Figure 8, it can be concluded that the ISO VG 460 oil performed better than other considered ISO grade oils and showed the minimum value for both considered objective functions. Comparing the compromise solutions obtained, the solution obtained using ISO VG 460 oil obtained 10% improvement in center distance and 1.2% improvement in power efficiency compared to the worst performing oil, e.g., ISO VG 100. The main reason for the poor performance of the other types of oil is lover viscosity, which lead to the lower power efficiency and increase of gearbox volume. This analysis emphasizes that the proper choice of lubricant does not only affect the power efficiency of gears, but also has impact on the overall dimensions of the gearbox.

4.2. The Benchmark Results

In this section, the optimization performance of the proposed MHPSODE algorithm has been evaluated and compared with well-known multi-objective evolutionary algorithms, such as NSGA-II [42], MO_Ring_PSO_SCD [22], and DEMO [43] on CEC2009 benchmark problems. The considered multi-objective optimization test instances released in CEC2009 benchmark consist of 10 unconstrained optimization problems and 3 problems that deal with constrained optimization [45]. In the following simulations the first 10 functions are taken into consideration. Furthermore, the DTLZ-I [46] test problem has also employed for analysis.

To analyze and compare the optimization performance of multiple algorithms, in this paper the Inverted Generational Distance (IGD) is employed as the performance metrics [45,47]. The IGD metric can be calculated based on the expression

$$IGD=\frac{{\sum }_{v_i\in P_A}d\left(v_i,P_A\right)}{\left|P^{*}\right|},$$

(50)

where $d\left(v_i,P_A\right)$ denotes the minimum Euclidean distance between $v_i$ and the points in the non-dominated solution set $P_A$ obtained by the tested algorithm, and $P^{*}$ is the set of uniformly distributed points in the solution space. To calculate IGD metrics, the 10,000 reference points have been sampled from the obtained Pareto fronts. It should be noted that smaller values of IGD metric points to better performing algorithm.

To obtain results for the statistical analysis all the experiments are independently run for 30 times and the standard deviation of the above-mentioned performance metrics has been calculated. To perform the statistical analysis, and determine the statistically significant differences between the obtained results, in this paper Wilcoxon Signed-Rank Test and Friedman test have been performed with the significance level $\alpha =0.05$.

First, to determine the significant differences between the performance of two algorithms, the Wilcoxon signed-rank test has been applied. Here, the sum of ranks for the problem, where first algorithm outperforms the second is denoted with $R^{+}$, while the sum of ranks, where second algorithm outperforms the first one is symbolized with $R^{-}$. In this regard, for the null hypothesis of the Wilcoxon signed-rank test it is assumed that: “there is no difference between the mean results of the two samples” [48]. Conversely, the alternative hypothesis is that: “there is a difference in the mean results of the two samples” [48]. The p value has been employed and compared to the significance level in the statistical analysis to test the hypothesis. Therefore, the null hypothesis can be rejected when $p\le \alpha$. According to the obtained results of the statistical analysis one of the following three signs $\left(+,-,\approx \right)$ have been given to the comparison of the algorithms. In this regard, the sign $\left(+\right)$ denotes that the first algorithm significantly outperforms the second algorithm, $\left(-\right)$ sign means that the first algorithm performed significantly worse than the first one, and $\left(\approx \right)$ denotes that the two algorithms taken for consideration have comparable performance.

Secondly, the Friedman test [48] has been performed to determine the significant differences in performance between all considered algorithms, which is achieved by comparing the ranks of all algorithms. In this way, the algorithm with a minimum rank value is denoted as the best performing algorithm, while the one with the highest rank is considered to be the worst performing algorithm. For the null hypothesis of the Friedman test it is assumed that: “there is no difference among the performance of all algorithms”, while the alternative hypothesis is “there is a difference among the performance of all algorithms” [48].

The results of numerical simulation computed over 30 independent runs over ten multi-objective test functions of the CEC2009 benchmark and DTLZ-I test problems are presented in

Table 3. The mean and standard deviation values of IDG metrics for different considered optimization algorithms.

Test Instances	DEMO	NSGA-II	MO_Ring_PSO_SCD	MHPSODE
	mean ± std
MMF1	$\underline{4.86\times {10}^{-2}\pm 6.48\times {10}^{-2}}$	$8.75\times {10}^{-2}\pm 1.42\times {10}^{-1}$	$6.43\times {10}^{-2}\pm 6.58\times {10}^{-2}$	$\mathbf{3.35}\times {\mathbf{10}}^{-\mathbf{2}}\pm \mathbf{3.86}\times {\mathbf{10}}^{-\mathbf{2}}$
MMF2	$1.89\times {10}^{-1}\pm 9.18\times {10}^{-2}$	$\underline{1.85\times {10}^{-1}\pm 7.75\times {10}^{-2}}$	$\mathbf{1.59}\times {\mathbf{10}}^{-\mathbf{2}}\pm \mathbf{5.32}\times {\mathbf{10}}^{-\mathbf{3}}$	$1.94\times {10}^{-1}\pm 7.17\times {10}^{-2}$
MMF3	$9.01\times {10}^{-1}\pm 2.98\times {10}^{-1}$	$2.46\times {10}^{-1}\pm 8.67\times {10}^{-2}$	$\underline{1.69\times {10}^{-1}\pm 1.40\times {10}^{-1}}$	$\mathbf{1.42}\times {\mathbf{10}}^{-\mathbf{1}}\pm \mathbf{1.68}\times {\mathbf{10}}^{-\mathbf{1}}$
MMF4	$2.48\times {10}^{-1}\pm 1.11\times {10}^{-1}$	$\mathbf{7.37}\times {\mathbf{10}}^{-\mathbf{2}}\pm \mathbf{3.07}\times {\mathbf{10}}^{-\mathbf{2}}$	$1.22\times {10}^{-1}\pm 6.79\times {10}^{-2}$	$\underline{9.25\times {10}^{-2}\pm 3.28\times {10}^{-2}}$
MMF5	$\mathbf{9.64}\times {\mathbf{10}}^{-\mathbf{3}}\pm \mathbf{1.56}\times {\mathbf{10}}^{-\mathbf{2}}$	$2.30\times {10}^{-2}\pm 3.01\times {10}^{-2}$	$6.13\times {10}^{-2}\pm 1.33\times {10}^{-1}$	$\underline{1.81\times {10}^{-2}\pm 5.46\times {10}^{-2}}$
MMF6	$4.62\times {10}^{-2}\pm 8.10\times {10}^{-2}$	$4.92\times {10}^{-2}\pm 6.94\times {10}^{-2}$	$\underline{4.43\times {10}^{-2}\pm 8.14\times {10}^{-2}}$	$\mathbf{4.36}\times {\mathbf{10}}^{-\mathbf{2}}\pm \mathbf{3.83}\times {\mathbf{10}}^{-\mathbf{2}}$
MMF7	$\underline{1.23\times {10}^{-2}\pm 4.06\times {10}^{-3}}$	$3.96\times {10}^{-2}\pm 5.86\times {10}^{-2}$	$6.51\times {10}^{-2}\pm 7.11\times {10}^{-2}$	$\mathbf{1.12}\times {\mathbf{10}}^{-\mathbf{2}}\pm \mathbf{3.60}\times {\mathbf{10}}^{-\mathbf{3}}$
MMF8	$2.40\times {10}^{-1}\pm 1.66\times {10}^{-1}$	$1.45\times {10}^{-1}\pm 1.00\times {10}^{-1}$	$\mathbf{8.58}\times {\mathbf{10}}^{-\mathbf{2}}\pm \mathbf{9.26}\times {\mathbf{10}}^{-\mathbf{2}}$	$\underline{9.22\times {10}^{-2}\pm 8.35\times {10}^{-2}}$
MMF9	$2.06\times {10}^{-1}\pm 2.11\times {10}^{-1}$	$1.25\times {10}^{-1}\pm 9.39\times {10}^{-2}$	$\underline{7.44\times {10}^{-2}\pm 4.93\times {10}^{-2}}$	$\mathbf{5.99}\times {\mathbf{10}}^{-\mathbf{2}}\pm \mathbf{4.03}\times {\mathbf{10}}^{-\mathbf{2}}$
MMF10	$4.57\times {10}^{-1}\pm 3.37\times {10}^{-1}$	$2.50\times {10}^{-1}\pm 1.21\times {10}^{-1}$	$\underline{1.21\times {10}^{-1}\pm 1.34\times {10}^{-1}}$	$\mathbf{8.41}\times {\mathbf{10}}^{-\mathbf{2}}\pm \mathbf{9.24}\times {\mathbf{10}}^{-\mathbf{2}}$
DTLZ-I	$1.69\times {10}^0\pm 1.97\times {10}^0$	$2.56\times {10}^1\pm 1.86\times {10}^1$	$3.31\times {10}^1\pm 7.71\times {10}^1$	$\mathbf{1.46}\times {\mathbf{10}}^{\mathbf{0}}\pm \mathbf{1.80}\times {\mathbf{10}}^{\mathbf{0}}$

in the form of mean and standard deviation values of IGD metrics. Here, the best results obtained are marked in bold and second best are underlined.

Based on the results in Table 2, it can be concluded that the proposed MHPSODE algorithm performed the best on functions $f_1$, $f_3$, $f_6$, $f_7$, $f_9$ and $f_{10}$, while it performed the second best on $f_4$, $f_5$ and $f_8$. It can be observed that only on $f_2$ the proposed algorithm did not perform satisfactorily, where MO_Ring_PSO_SCD is the best performing algorithm while NSGA-II is the second best. Furthermore, it can be observed that MHPSODE achieved the best performance on DTLZ-I test problem.

To perform the statistical pair-wise comparison of the optimization performance of the proposed MHPSODE and other considered algorithms, the obtained results have been analyzed using Wilcoxon signed-rank test. In this regard, in

Table 4. The statistical results of the Wilcoxon’s test between HPSODE and other considered algorithms on CEC2009 benchmark.

Algorithms	${\mathit{R}}^{+}$	${\mathit{R}}^{-}$	p-Value	Dec.
MHPSODE vs. DEMO	30,818	14,333	0.00	+
MHPSODE vs. MO_Ring_PSO_SCD	24,799	20,352	0.1345	≈
MHPSODE vs. NSGA-II	27,930	17,221	0.000	+

the statistical results of applying Wilcoxon’s signed-rank test on CEC2009 benchmark problems have been presented.

From the results in Table 4, it can be concluded that the proposed MHPSODE algorithm has significantly better performance than DEMO and NSGA-II algorithms. However, in comparison with the MO_Ring_PSO_SCD algorithm it is observed that proposed algorithm achieves similar performance, while in this case MHPSODE algorithm still provides higher $R^{+}$ values than $R^{-}$ in all considered cases.

Furthermore, the numerical experiments have been carried out to study the effectiveness of hybridization proposed in this paper by comparing the overall performance of all considered algorithms. In this regard, the average ranks according to Friedman for the considered algorithms for different objective functions of CEC2009 benchmark problems have been presented in

Table 5. Average ranks computed using the Friedman statistical test for all algorithms across all problems using CEC2009 benchmark, with the significance level $\alpha =0.05$.

Algorithm	MMF1	MMF2	MMF3	MMF4	MMF5	MMF6	MMF7	MMF8	MMF9	MMF10	Mean Ranking	Rank
MHPSODE	2.93	1.43	1.66	2.13	2.5	1.8	3.4	1.96	2.00	1.73	2.16	1
MO_Ring_PSO_SCD	2.76	1.56	2.00	2.47	3.13	2.63	3.26	2.00	2.3	1.87	2.4	2
DEMO	2.5	3.5	2.6	1.53	2.73	2.6	2.13	2.63	3.03	2.97	2.62	3
NSGA-II	1.8	3.5	3.73	3.86	1.63	2.97	1.2	3.4	2.66	3.43	2.82	4
Friedman p value	0	0	0	0	0	0	0	0	0	0

.

2 In this table, the rank of the best performing algorithm is denoted in bold, and the second best is underlined. From the results in Table 5, it can be observed that the p values computed using the Friedman test for different objective functions are less than significance $\alpha =0.05$. This means that there is significant difference between the performance of different considered algorithms, which agrees with the conclusions drawn from statistical comparison using Wilcoxon’s test. Therefore, the obtained results of the statistical comparison demonstrate the effectiveness of the hybrid algorithm proposed in this paper.

Finally, the average computation times of the considered algorithms in searching for the global optimum of the MOO problem have been determined on the computer with 1.6 GHz CPU and 8 GB of RAM. Therefore, the comparison between the average computation times of the MHPSODE, MO_Ring_PSO_SCD, DEMO and NSGA-II algorithms in solving considered CEC2009 benchmark problems are shown in

Table 6. Average computation times of the considered algorithms on CEC2009 benchmark functions.

Algorithm	MMF1	MMF2	MMF3	MMF4	MMF5	MMF6	MMF7	MMF8	MMF9	MMF10	Total avg. Time
MHPSODE	$1.79\times {10}^2$	$2.08\times {10}^2$	$7.72\times {10}^1$	$1.50\times {10}^1$	$1.04\times {10}^1$	$1.77\times {10}^1$	$1.65\times {10}^2$	$2.39\times {10}^1$	$8.83\times {10}^1$	$1.10\times {10}^1$	$7.96\times {10}^1$
NSGA-II	$3.47\times {10}^1$	$2.16\times {10}^1$	$3.01\times {10}^1$	$1.87\times {10}^1$	$5.59\times {10}^1$	$4.33\times {10}^1$	$3.86\times {10}^1$	$1.90\times {10}^1$	$1.86\times {10}^1$	$1.90\times {10}^1$	$2.99\times {10}^1$
DEMO	$2.94\times {10}^0$	$2.89\times {10}^0$	$3.21\times {10}^0$	$2.79\times {10}^0$	$4.13\times {10}^0$	$3.34\times {10}^0$	$3.26\times {10}^0$	$5.79\times {10}^0$	$3.21\times {10}^0$	$3.75\times {10}^0$	$3.53\times {10}^0$
MO_Ring_PSO_SCD	$2.25\times {10}^2$	$1.86\times {10}^2$	$1.61\times {10}^2$	$1.39\times {10}^2$	$2.50\times {10}^2$	$2.60\times {10}^2$	$2.38\times {10}^2$	$1.88\times {10}^2$	$1.93\times {10}^2$	$1.76\times {10}^2$	$2.02\times {10}^2$

.

From the results presented in Table 6, it can be observed that the DEMO algorithm has the fastest implementation among the considered algorithms, while the MO_Ring_PSO_SCD algorithm has the worst implementation. The NSGA-II algorithm has consistent average execution times among all considered MOO functions, while the execution time of proposed MHPSODE algorithm depends on the complexity of the problem which is solved. Since the considered planetary gearbox optimization problem is not a real-time task, and based on the conclusions drawn from the analysis of the performance of the algorithms on MOO problems, the higher execution time of the proposed MHPSODE provides a compromise between the required design parameters of planetary gearbox and average computation time.

5. Conclusions

In this paper, the multi-objective optimization of the planetary gear train has been considered with multiple objective functions, such as center distance, volume of the gears, gearbox efficiency, bending stress etc. The theoretical formulation and numerical procedure for the calculation of the planetary gearbox power efficiency has been formulated as an important objective in the considered problem. To solve this complex MOO problem in this paper the hybrid algorithm, named MHPSODE, as the hybridization between DE and PSO algorithms has been proposed. The hybridization has been performed by introducing mutation operators from DE algorithm into the velocity update equation of the PSO algorithm, where the appropriate mutation operator is applied based on an adaptive normalized population spacing parameter.

Using the proposed algorithm, the non-convex Pareto frontier has been obtained, which represents the relationship between conflicting objectives. The Pareto front obtained in analysis is convenient for the designer to select the appropriate parameters of the construction. The results were validated and compared with the parameters of the reference industrial gearbox, which showed the maximum 30% reduction in center distance, 23.6% improvement in gearbox volume and 2.2% improvement in efficiency. In addition, to determine the best performing oil for the formulated problem the analysis has been performed with multiple ISO grade oils, which concluded that ISO VG 460 achieves the best trade-off between efficiency and gearbox dimensions. Furthermore, to evaluate the optimization performance of the proposed MHPSODE algorithm the performance has been compared on the CEC 2009 benchmark problems. Compared to the well-known algorithms, the proposed method obtained the better convergence and more accurate Pareto front on the considered optimization problems. Therefore, the proposed hybrid algorithm applied on the considered MOO problem of planetary gearbox optimization can provide the design parameters which lead to improved efficiency, reduction in the total weight and enhanced service life of the planetary gearbox.