# Load Sharing Multiobjective Optimization Design of a Split Torque Helicopter Transmission

Authors: Chenxi Fu; Ning Zhao; Yongzhi Zhao
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Split torque designs can offer significant advantages over the traditional planetary designs for helicopter transmissions. However, it has two unique properties, gap and phase differences, which result in the risk of unequal load sharing. Various methods have been proposed to eliminate the effect of gap and promote load sharing to a certain extent. In this paper, system design parameters will be optimized to change the phase difference, thereby further improving load sharing. A nonlinear dynamic model is established to measure the load sharing with dynamic mesh forces quantitatively. Afterwards, a multiobjective optimization of a reference split torque design is conducted with the promoting of load sharing property, lightweight, and safety considered as the objectives. The load sharing property, which is measured by load sharing coefficient, is evaluated under multiple operating conditions with dynamic analysis method. To solve the multiobjective model with NSGA-II, an improvement is done to overcome the problem of time consuming. Finally, a satisfied optimal solution is picked up as the final design from the Pareto optimal front, which achieves improvements in all the three objectives compared with the reference design.


---

## Body

## 1. Introduction

The drive system of a rotorcraft must meet especially demanding requirements of high reduction ratio, high safety and reliability, lightweight, and little vibration. White [1] proposed a promising alternative design to the common planetary transmissions for helicopters, known as the split torque or split path arrangement, which offers two parallel paths for transmitting torque from the engine to the rotor. Kish [2] and Krantz [3, 4] pointed out that although a split torque design can offer significant advantages over the commonly used planetary design, the split torque design is considered even risk: there might have been unequal torques in the two parallel paths, which cause excessive wear in one of the paths and renders the split torque system ineffective. Therefore, the main problem in the design of split torque transmission lies in how to ensure the torque splitting equally between different paths.Actually the problem of load sharing exists in almost all types of multipath transmission system, such as split torque and planetary designs. There is a locked loop of gearing in a multipath design, and each of the gear meshes in the loop will be engaged at the same time, only if the gears are assembled with the required relative relationship, which is considered a particular tooth timing relationship. However, due to the manufacturing and assembly errors, the particular tooth timing relationship can not precisely be met in a real multipath transmission system. Therefore, there might be a gap (or several gaps) at one (or several) of the gear mesh locations under a nominal light load. Usually this kind of gap is considered the main cause of unequaled load in multipath transmission systems [3].To analyze the load sharing property of a multipath transmission system, a lot of research works have been conducted. The main research contents and methods include defining the load sharing coefficient based on mesh forces, which are evaluated through dynamic or static methods, to measure the load sharing property of a multipath transmission system, and then effects of different factors such as errors and operating conditions can be studied. Hayashi et al. [5] put forward a method to measure the load sharing of planetary transmissions and concluded that the dynamic load sharing differs greatly from static condition. Therefore, the dynamic analysis method is widely used because it reflects the real load sharing property under operating condition and it is naturally adopted in this paper. Krantz [3] firstly studied the dynamics of a split torque transmission system and concluded that the loads and motions of the two power paths differ although the system has symmetric geometry. Kahraman [6] investigated the load sharing property of a planetary transmission system with a nonlinear dynamic model established, which takes manufacturing and assembly errors into consideration, and finds that the operating conditions also having great effects on load sharing. Guo et al. [7] studied the influences of related factors on load sharing property of the wind turbine planetary gears based on its characteristics. Kahraman [6] derived the relationship between dynamic load sharing coefficient and static load sharing coefficient based on dynamic method, which makes static analyzing meaningful [8]. Bodas and Kahraman [9] and Singh [10] studied the influences of manufacturing and assembly errors on load sharing property of planetary transmissions, with 2D and 3D static contact models adopted, respectively. Afterwards experimental studies [11, 12] have also been conducted. Ligata et al. [13] established a discrete model of planetary transmissions to study its load sharing property. Then Singh [14, 15] investigated the model deeply and obtains multiple load sharing coefficients under varies of manufacturing errors and torques by defining the load sharing map.Benefiting from the centrosymmetry of design, the planetary multipath transmission systems are always in capacity of automatical load sharing, which offers an advantage over the split torque design. To promote the load sharing property of split torque transmission system, some methods have been proposed to compensate for or minimize the effect of the gap, including floating gears [16], quill shafts [2], and “clocking angle” [4, 17]. The floating gears arrangement permits the input pinion to float until gear loads are balanced between the two paths. There are two kinds of quill shafts load sharing devices: the conventional quill shafts and the one based on elastomeric elements. The conventional quill shafts assemble intermediate shafts with some torsional flexibility so as to minimize the difference in torque split between paths, whereas the latter one add some materials with a lower elastic modulus in the compound shafts to achieve the same purpose. The “clocking angle” method considers the “clocking angle” as a design parameter to adjust and optimize the load sharing, and the “clocking angle” is adjusted by varying the thicknesses of shim packs that axially positioned the compound shafts.Further research of Krantz [3] indicated that even though the manufacturing and assembly errors of a split torque transmission are precisely controlled and the gap is eliminated completely, there still unequaled loads in the two paths under real operating conditions. Actually the unique phase difference property of the split torque transmission system results in desynchrony between the two paths, which causes significant effects on load sharing. However, the phase differences are not independent design parameters, which are connected to the system geometry parameters. Therefore, the load sharing property of a split torque transmission can be promoted by adjusting the phase differences through optimizing the system geometry parameters. Since it reflects the real load sharing of the split torque transmission under dynamic condition, the load sharing coefficient used in this paper is evaluated by solving the nonlinear dynamic model.Despite of the load sharing property, the drive system of a rotorcraft must also meet the demanding requirements of lightweight and high safety [2, 3]. Savsani et al. [18] and Thompson et al. [19] reduced the mass of a gear pair and a multistage gear system, respectively, through optimizing the gear parameters. Kumar et al. [20] optimized a single pair of gear transmission with promoting of load capacity of gears considered as the objective. Based on the above considerations, a multiobjective optimization design of a split torque transmission system is conducted, with the promoting of load sharing property, lightweight, and safety considered as the objectives. The load sharing property, which is measured by load sharing coefficient, is evaluated under multiple operating conditions with dynamic analysis method.Deb et al. [21] proposed the improved nondominated sorting genetic algorithm (NSGA-II) based on multiobjective evolutionary algorithm for multiobjective optimization problems, which has been proved a simple and effective method [22]. In this paper, NSGA-II is adopted to solve the multiobjective optimization model. However, there are large numbers of nonlinear dynamic equations to be solved under multiple operating conditions when evaluating the fitness, and the solving time is not acceptable in engineering. Therefore, an improvement has been done to NSGA-II to solve the problem of time consuming: prediction strategy is used in the fitness evaluation step so as to avoid the evaluation of load sharing property which is computationally very expensive.


## 2. Nonlinear Dynamics Model of a Split Torque Transmission System

### 2.1. Modeling the Gear Mesh

Regarding the spur gears as special helical gears with helix angle of 0 degree, the general 3D pinion-gear meshing model can be established, as shown in Figure1. Assuming that the geometry is not affected by deflections (small displacements hypothesis) and provided that mesh elasticity can be transferred onto the base plane, a rigid-body approach can be employed. The pinion and the gear can therefore be assimilated to two rigid cylinders with 4 degrees of freedom each, which are connected by a stiffness element (or a distribution of stiffness elements) and a damp element. From a physical point of view, the 8 degrees of freedom of a pair represent the generalized displacements of 3 translational degrees (along axes x, y, and z) and 1 rotational degree (around axis z) each. Here, the angle ψ between axis x and the center line of the gear pair OpOg→ is defined as direction angle.Figure 1
Model of gear contact.In this model, flank deviations are taken into consideration. Conventionally flank deviations relative to the perfect geometry are positively defined in the direction of the outer normal and they are supposed to be small enough so that tooth contacts remain on theoretical base planes. Then the deviation of a gear pair can be defined as sum of flank deviations of the pinion and gear. Each theoretical contact line on the base plane is discretized in elementary cells centered in one point, whose deviation is defined as the normal distance on the base plane between a point of the pinion and a point of the gear that would be in contact for perfect geometries [23]. In fact, due to the influence of deviations, contact and deflection do not occur at all the point in the theoretical contact line. The real contact status of the gear pair can be obtained only when the contact status of all the point in all the contact lines is solved.As Figure1 illustrates, when a gear pair is nominally engaged, contact occurs at one certain point M∗, which has the maximum deviation, namely, eM∗. Select any point M in the theoretical contact line different from M∗, and name the flank deviations of M in the pinion and gear as efpM and efgM, respectively. The flank deviation is positive for an excess of material and negative when some material is removed from the ideal geometry. Therefore, the deviation in point M can be expressed as(1)eM=efpM+efgM.The gear model is shown in Figure1, and pinion and gear are assimilated to rigid cylinders with four-degrees-of-freedom connected by series of stiffness and damp. The generalized displacements can be expressed as(2)q=xp,yp,zp,θp,xg,yg,zg,θgT,where xp, yp, zp are translational displacements of the pinion along axes x, y, z, respectively. θp is rotational displacement of the pinion around axis z; xg, yg, zg are translational displacements of the gear along axes x, y, z, respectively. θg is rotational displacement of the gear around axis z.Since generalized displacements are small quantities, it can therefore be represented by infinitesimal translations and rotations. The perturbationq on pinion and gear gives a normal approach δM on the base plane relative to rigid body positions:(3)δM=κTq=cos⁡βbsin⁡α-ψcos⁡βbcos⁡α-ψsin⁡βbRbpcos⁡βb-cos⁡βbsin⁡α-ψ-cos⁡βbcos⁡α-ψ-sin⁡βbRbgcos⁡βbT·xpypzpθpxgygzgθg,where κ is a projective vector depending on gear geometry which projects displacements q onto the base plane; βb is the base helix angle (the value is 0 to spur gears); α is the pressure angle; Rbp, Rbg are the base radius of pinion and gear.Contact and deflection inM occur only if the normal approach δM is larger than the initial deviation δeM, and in such case the deflection ΔM can be evaluated as(4)ΔM=δM-δeM,where δeM=eM∗-eM, eM∗, eM is the deviation in M∗, M. Else if the normal approach δM is less than the initial deviation δeM, then the deflection ΔM=0.The elemental mesh force transmitted from the pinion onto the gear at one point of contactM is(5)dFM=kMdl·ΔM,where kM is the mesh stiffness at point M per unit of contact length which is calculated according to ISO6336 and dl is the elemental contact length.The total mesh forceFl can be deducted by integrating over the time-varying and deflection dependent contact length:(6)Fl=∫lkMΔMdl=∫lkMdl·κTq-∫lkMδeMdl.


### 2.2. Gear Train Arrangement

The full arrangement of a split torque or split-path transmissions [3] is depicted in Figure 2, which has two stages. (1) One is first reduction stage. The first stage is where torque is split between the input pinion and the two output gears. Usually helical gears are used. (2) The second stage is high torque reduction stage. The output shaft is driven by a gear which is driven simultaneously by two spur pinions, each coaxial to the gear in the first reduction stage.Figure 2
Full arrangement of a split torque transmission system.As shown in Figure2, the input pinion meshes with two gears, offering two paths to transfer power to the output gear; therefore, the whole system is divided into two paths. The two power paths are identified as L and R, with R to the right of L. The first-stage gear and second-stage pinion combination are collectively called the compound gear. The compound gear and gear shaft combination are called the compound shaft.As shown in Figure3, a right-hand Cartesian coordinate system is established such that the z-axis is coincident with the output gear shaft, the positive y-axis extends from the input gear center to the output pinion center, and the input gear drives clockwise. The first-stage pinion, gear (L), and gear (R) are marked with gears 1~3; the second-stage pinion (L), pinion (R), and gear are marked with gears 4~6. The input shaft, compound shaft (L), compound shaft (R), and output shaft are marked with axes 1~4, and the center of which is marked with O1, O2, O3, O4. Define the angle between O1O2→ and O1O3→ as shaft angle Φ1 and the angle between O4O2→ and O4O3→ as shaft angle Φ2.Figure 3
Coordinates of a split torque transmission system.A unique property of a split torque transmission is the phase relationships of the meshes. The input pinion drives two gears simultaneously, the length of the arc along the pitch circle joining the two pitch pointsa, b (the length of arc ab in Figure 3) is probably not an integer multiple of the circular pitch, and then the two meshes will not pass through the pitch point at the same instant of time, which results in the desynchrony of the two paths. It is the same with the second stage. In order to describe the unique property of a split torque transmission, the concept of phase difference of a stage is defined as quotient of the length of arc, joining the two pitch points and the circular pitch. The phase differences of the two stages are expressed as(7)Δφ1=ab⌢pt1=d1Φ1/2πmt1=Z1Φ12π,Δφ2=cd⌢pt2=d6Φ2/2πmt2=Z6Φ22π,where Δφ1,2 is the phase difference of the first and second stages; ab⌢, cd⌢ are the length of arcs ab, cd; pt1, pt2 are the circular pitch of the first and second stages; mt1, mt2 are the transverse module of the first and second stages; di is the pitch diameter of gear i; Zi is the tooth number of gear i.There are 6 gears with 4 meshing pairs in a split torque transmission system. Each gear can be assimilated to a rigid cylinder with 4 degrees of freedom. Therefore, there are 24 degrees of freedom in total. The generalized displacements of the whole system can be expressed as(8)x=x1,y1,z1,θ1,x2,y2,z2,θ2,…,x6,y6,z6,θ6T,where the subscripts 1~6 correspond to gears 1~6.


### 2.3. Equations of Dynamics

The general form of dynamic equations of a gear transmission system is(9)Mx¨+Cx˙+Kx=F,where M is the generalized mass matrix, C is the damp matrix, K is the stiffness matrix, and F is the load vector.With regard to the split torque transmission system, the general model can be specified to(10)Mx¨+Cx˙+Km+Kb+Kcx=F0+Fe,where Km is the mesh stiffness term, Kb is the supporting stiffness term, Kc is the coupling stiffness term, F0 is the constant torques term, and Fe is the additional force term caused by gear deviations.Here we give the principle to set up all the matrices and vector in the model.


#### 2.3.1. Mass MatrixM

Mass matrixM is a 24×24 diagonal matrix, which is expressed as(11)M=diag⁡m1,m1,m1,I1,m2,m2,m2,I2,…,m6,m6,m6,I6,where mi, i=1~6 is the mass of gears 1~6, Ii, i=1~6 is the rotational inertia around z-axis.


#### 2.3.2. Stiffness MatrixK

The stiffness matrixK consists of mesh stiffness Km, supporting stiffness Kb, and coupling stiffness Kc, and it is also a 24×24 symmetric matrix:(12)K=Km+Kb+Kc.(i) Meshing Stiffness Matrix
K
m
. The mesh stiffness matrix of the whole system Km comes from mesh stiffness of all the gear pairs; it can therefore be obtained by assembling the entire mesh stiffness submatrix together:(13)Km=∑n=14RijTKm,ijRiji,j=1,2;1,3;4,6;5,6,where Km,ij=∫lkijMdl·κijκijT is the nonlinear and time-dependent mesh stiffness submatrix of gear pair i, j, with κij: the projective vector of gear pair i and j:(14)Rij=columnicolumnj04×4⋯hhhhR1i⋯hh04×4⋯04×404×4⋯hhhh04×4⋯hhR2j⋯04×42×6is assemble matrix, with R1i=R2j=diag1,1,1,1.(ii) Supporting Stiffness Matrix
K
b
. The supporting stiffness matrix Kb comes from all the gears’ supporting stiffness in all the translational motion freedoms:(15)Kb=diag⁡kx1,ky1,kz1,0,kx2,ky2,kz2,0,…,kx6,ky6,kz6,0,where kxi, kyi, kzi, i=1~6 is the supporting stiffness of gear i in x, y, z directions.(iii) Coupling Stiffness Matrix
K
c
. The coupling stiffness matrix Kc of the whole system comes from the two compound shafts, and it can be obtained by assembling the two coupling stiffness submatrix together:(16)Kc=∑n=12RijTKc,ijRij,i,j=2,4;3,5,where(17)Kc,ij=kb,ij000-kb,ij000kb,ij000-kb,ij00ka,ij000-ka,ij0ks,ij000-ks,ijkb,ij000kb,ij00symka,ij0ks,iji,j=2,4;3,5is the coupling stiffness submatrix of compound gear i, j, with ks,ij, kb,ij, ka,ij: the torsional stiffness, bending stiffness, axial stiffness between the compound gears i and j.


#### 2.3.3. Damp MatrixC

Here Rayleigh’s damping is adopted:(18)C=aM+bK,with a, b: two constants to be adjusted from experimental results and experience.


#### 2.3.4. Load VectorF

The load vectorF consists of torques vector F0 produced by input and output torques and additional force vector Fe caused by gear deviations:(19)F=F0+Fe.The torques vectorF0 is expressed as(20)F0=0,0,0,Tin,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,ToutT.The additional force vectorFe caused by gear deviations is given by(21)Fe=∑n=14RijT∫lkijMδeijMdL·κiji,j=1,2;1,3;4,6;5,6.


## 2.1. Modeling the Gear Mesh

Regarding the spur gears as special helical gears with helix angle of 0 degree, the general 3D pinion-gear meshing model can be established, as shown in Figure1. Assuming that the geometry is not affected by deflections (small displacements hypothesis) and provided that mesh elasticity can be transferred onto the base plane, a rigid-body approach can be employed. The pinion and the gear can therefore be assimilated to two rigid cylinders with 4 degrees of freedom each, which are connected by a stiffness element (or a distribution of stiffness elements) and a damp element. From a physical point of view, the 8 degrees of freedom of a pair represent the generalized displacements of 3 translational degrees (along axes x, y, and z) and 1 rotational degree (around axis z) each. Here, the angle ψ between axis x and the center line of the gear pair OpOg→ is defined as direction angle.Figure 1
Model of gear contact.In this model, flank deviations are taken into consideration. Conventionally flank deviations relative to the perfect geometry are positively defined in the direction of the outer normal and they are supposed to be small enough so that tooth contacts remain on theoretical base planes. Then the deviation of a gear pair can be defined as sum of flank deviations of the pinion and gear. Each theoretical contact line on the base plane is discretized in elementary cells centered in one point, whose deviation is defined as the normal distance on the base plane between a point of the pinion and a point of the gear that would be in contact for perfect geometries [23]. In fact, due to the influence of deviations, contact and deflection do not occur at all the point in the theoretical contact line. The real contact status of the gear pair can be obtained only when the contact status of all the point in all the contact lines is solved.As Figure1 illustrates, when a gear pair is nominally engaged, contact occurs at one certain point M∗, which has the maximum deviation, namely, eM∗. Select any point M in the theoretical contact line different from M∗, and name the flank deviations of M in the pinion and gear as efpM and efgM, respectively. The flank deviation is positive for an excess of material and negative when some material is removed from the ideal geometry. Therefore, the deviation in point M can be expressed as(1)eM=efpM+efgM.The gear model is shown in Figure1, and pinion and gear are assimilated to rigid cylinders with four-degrees-of-freedom connected by series of stiffness and damp. The generalized displacements can be expressed as(2)q=xp,yp,zp,θp,xg,yg,zg,θgT,where xp, yp, zp are translational displacements of the pinion along axes x, y, z, respectively. θp is rotational displacement of the pinion around axis z; xg, yg, zg are translational displacements of the gear along axes x, y, z, respectively. θg is rotational displacement of the gear around axis z.Since generalized displacements are small quantities, it can therefore be represented by infinitesimal translations and rotations. The perturbationq on pinion and gear gives a normal approach δM on the base plane relative to rigid body positions:(3)δM=κTq=cos⁡βbsin⁡α-ψcos⁡βbcos⁡α-ψsin⁡βbRbpcos⁡βb-cos⁡βbsin⁡α-ψ-cos⁡βbcos⁡α-ψ-sin⁡βbRbgcos⁡βbT·xpypzpθpxgygzgθg,where κ is a projective vector depending on gear geometry which projects displacements q onto the base plane; βb is the base helix angle (the value is 0 to spur gears); α is the pressure angle; Rbp, Rbg are the base radius of pinion and gear.Contact and deflection inM occur only if the normal approach δM is larger than the initial deviation δeM, and in such case the deflection ΔM can be evaluated as(4)ΔM=δM-δeM,where δeM=eM∗-eM, eM∗, eM is the deviation in M∗, M. Else if the normal approach δM is less than the initial deviation δeM, then the deflection ΔM=0.The elemental mesh force transmitted from the pinion onto the gear at one point of contactM is(5)dFM=kMdl·ΔM,where kM is the mesh stiffness at point M per unit of contact length which is calculated according to ISO6336 and dl is the elemental contact length.The total mesh forceFl can be deducted by integrating over the time-varying and deflection dependent contact length:(6)Fl=∫lkMΔMdl=∫lkMdl·κTq-∫lkMδeMdl.


## 2.2. Gear Train Arrangement

The full arrangement of a split torque or split-path transmissions [3] is depicted in Figure 2, which has two stages. (1) One is first reduction stage. The first stage is where torque is split between the input pinion and the two output gears. Usually helical gears are used. (2) The second stage is high torque reduction stage. The output shaft is driven by a gear which is driven simultaneously by two spur pinions, each coaxial to the gear in the first reduction stage.Figure 2
Full arrangement of a split torque transmission system.As shown in Figure2, the input pinion meshes with two gears, offering two paths to transfer power to the output gear; therefore, the whole system is divided into two paths. The two power paths are identified as L and R, with R to the right of L. The first-stage gear and second-stage pinion combination are collectively called the compound gear. The compound gear and gear shaft combination are called the compound shaft.As shown in Figure3, a right-hand Cartesian coordinate system is established such that the z-axis is coincident with the output gear shaft, the positive y-axis extends from the input gear center to the output pinion center, and the input gear drives clockwise. The first-stage pinion, gear (L), and gear (R) are marked with gears 1~3; the second-stage pinion (L), pinion (R), and gear are marked with gears 4~6. The input shaft, compound shaft (L), compound shaft (R), and output shaft are marked with axes 1~4, and the center of which is marked with O1, O2, O3, O4. Define the angle between O1O2→ and O1O3→ as shaft angle Φ1 and the angle between O4O2→ and O4O3→ as shaft angle Φ2.Figure 3
Coordinates of a split torque transmission system.A unique property of a split torque transmission is the phase relationships of the meshes. The input pinion drives two gears simultaneously, the length of the arc along the pitch circle joining the two pitch pointsa, b (the length of arc ab in Figure 3) is probably not an integer multiple of the circular pitch, and then the two meshes will not pass through the pitch point at the same instant of time, which results in the desynchrony of the two paths. It is the same with the second stage. In order to describe the unique property of a split torque transmission, the concept of phase difference of a stage is defined as quotient of the length of arc, joining the two pitch points and the circular pitch. The phase differences of the two stages are expressed as(7)Δφ1=ab⌢pt1=d1Φ1/2πmt1=Z1Φ12π,Δφ2=cd⌢pt2=d6Φ2/2πmt2=Z6Φ22π,where Δφ1,2 is the phase difference of the first and second stages; ab⌢, cd⌢ are the length of arcs ab, cd; pt1, pt2 are the circular pitch of the first and second stages; mt1, mt2 are the transverse module of the first and second stages; di is the pitch diameter of gear i; Zi is the tooth number of gear i.There are 6 gears with 4 meshing pairs in a split torque transmission system. Each gear can be assimilated to a rigid cylinder with 4 degrees of freedom. Therefore, there are 24 degrees of freedom in total. The generalized displacements of the whole system can be expressed as(8)x=x1,y1,z1,θ1,x2,y2,z2,θ2,…,x6,y6,z6,θ6T,where the subscripts 1~6 correspond to gears 1~6.


## 2.3. Equations of Dynamics

The general form of dynamic equations of a gear transmission system is(9)Mx¨+Cx˙+Kx=F,where M is the generalized mass matrix, C is the damp matrix, K is the stiffness matrix, and F is the load vector.With regard to the split torque transmission system, the general model can be specified to(10)Mx¨+Cx˙+Km+Kb+Kcx=F0+Fe,where Km is the mesh stiffness term, Kb is the supporting stiffness term, Kc is the coupling stiffness term, F0 is the constant torques term, and Fe is the additional force term caused by gear deviations.Here we give the principle to set up all the matrices and vector in the model.


### 2.3.1. Mass MatrixM

Mass matrixM is a 24×24 diagonal matrix, which is expressed as(11)M=diag⁡m1,m1,m1,I1,m2,m2,m2,I2,…,m6,m6,m6,I6,where mi, i=1~6 is the mass of gears 1~6, Ii, i=1~6 is the rotational inertia around z-axis.


### 2.3.2. Stiffness MatrixK

The stiffness matrixK consists of mesh stiffness Km, supporting stiffness Kb, and coupling stiffness Kc, and it is also a 24×24 symmetric matrix:(12)K=Km+Kb+Kc.(i) Meshing Stiffness Matrix
K
m
. The mesh stiffness matrix of the whole system Km comes from mesh stiffness of all the gear pairs; it can therefore be obtained by assembling the entire mesh stiffness submatrix together:(13)Km=∑n=14RijTKm,ijRiji,j=1,2;1,3;4,6;5,6,where Km,ij=∫lkijMdl·κijκijT is the nonlinear and time-dependent mesh stiffness submatrix of gear pair i, j, with κij: the projective vector of gear pair i and j:(14)Rij=columnicolumnj04×4⋯hhhhR1i⋯hh04×4⋯04×404×4⋯hhhh04×4⋯hhR2j⋯04×42×6is assemble matrix, with R1i=R2j=diag1,1,1,1.(ii) Supporting Stiffness Matrix
K
b
. The supporting stiffness matrix Kb comes from all the gears’ supporting stiffness in all the translational motion freedoms:(15)Kb=diag⁡kx1,ky1,kz1,0,kx2,ky2,kz2,0,…,kx6,ky6,kz6,0,where kxi, kyi, kzi, i=1~6 is the supporting stiffness of gear i in x, y, z directions.(iii) Coupling Stiffness Matrix
K
c
. The coupling stiffness matrix Kc of the whole system comes from the two compound shafts, and it can be obtained by assembling the two coupling stiffness submatrix together:(16)Kc=∑n=12RijTKc,ijRij,i,j=2,4;3,5,where(17)Kc,ij=kb,ij000-kb,ij000kb,ij000-kb,ij00ka,ij000-ka,ij0ks,ij000-ks,ijkb,ij000kb,ij00symka,ij0ks,iji,j=2,4;3,5is the coupling stiffness submatrix of compound gear i, j, with ks,ij, kb,ij, ka,ij: the torsional stiffness, bending stiffness, axial stiffness between the compound gears i and j.


### 2.3.3. Damp MatrixC

Here Rayleigh’s damping is adopted:(18)C=aM+bK,with a, b: two constants to be adjusted from experimental results and experience.


### 2.3.4. Load VectorF

The load vectorF consists of torques vector F0 produced by input and output torques and additional force vector Fe caused by gear deviations:(19)F=F0+Fe.The torques vectorF0 is expressed as(20)F0=0,0,0,Tin,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,ToutT.The additional force vectorFe caused by gear deviations is given by(21)Fe=∑n=14RijT∫lkijMδeijMdL·κiji,j=1,2;1,3;4,6;5,6.


## 2.3.1. Mass MatrixM

Mass matrixM is a 24×24 diagonal matrix, which is expressed as(11)M=diag⁡m1,m1,m1,I1,m2,m2,m2,I2,…,m6,m6,m6,I6,where mi, i=1~6 is the mass of gears 1~6, Ii, i=1~6 is the rotational inertia around z-axis.


## 2.3.2. Stiffness MatrixK

The stiffness matrixK consists of mesh stiffness Km, supporting stiffness Kb, and coupling stiffness Kc, and it is also a 24×24 symmetric matrix:(12)K=Km+Kb+Kc.(i) Meshing Stiffness Matrix
K
m
. The mesh stiffness matrix of the whole system Km comes from mesh stiffness of all the gear pairs; it can therefore be obtained by assembling the entire mesh stiffness submatrix together:(13)Km=∑n=14RijTKm,ijRiji,j=1,2;1,3;4,6;5,6,where Km,ij=∫lkijMdl·κijκijT is the nonlinear and time-dependent mesh stiffness submatrix of gear pair i, j, with κij: the projective vector of gear pair i and j:(14)Rij=columnicolumnj04×4⋯hhhhR1i⋯hh04×4⋯04×404×4⋯hhhh04×4⋯hhR2j⋯04×42×6is assemble matrix, with R1i=R2j=diag1,1,1,1.(ii) Supporting Stiffness Matrix
K
b
. The supporting stiffness matrix Kb comes from all the gears’ supporting stiffness in all the translational motion freedoms:(15)Kb=diag⁡kx1,ky1,kz1,0,kx2,ky2,kz2,0,…,kx6,ky6,kz6,0,where kxi, kyi, kzi, i=1~6 is the supporting stiffness of gear i in x, y, z directions.(iii) Coupling Stiffness Matrix
K
c
. The coupling stiffness matrix Kc of the whole system comes from the two compound shafts, and it can be obtained by assembling the two coupling stiffness submatrix together:(16)Kc=∑n=12RijTKc,ijRij,i,j=2,4;3,5,where(17)Kc,ij=kb,ij000-kb,ij000kb,ij000-kb,ij00ka,ij000-ka,ij0ks,ij000-ks,ijkb,ij000kb,ij00symka,ij0ks,iji,j=2,4;3,5is the coupling stiffness submatrix of compound gear i, j, with ks,ij, kb,ij, ka,ij: the torsional stiffness, bending stiffness, axial stiffness between the compound gears i and j.


## 2.3.3. Damp MatrixC

Here Rayleigh’s damping is adopted:(18)C=aM+bK,with a, b: two constants to be adjusted from experimental results and experience.


## 2.3.4. Load VectorF

The load vectorF consists of torques vector F0 produced by input and output torques and additional force vector Fe caused by gear deviations:(19)F=F0+Fe.The torques vectorF0 is expressed as(20)F0=0,0,0,Tin,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,ToutT.The additional force vectorFe caused by gear deviations is given by(21)Fe=∑n=14RijT∫lkijMδeijMdL·κiji,j=1,2;1,3;4,6;5,6.


## 3. Case Study and Load Sharing Property

The load sharing properties of a split torque transmission are how equivalent the load allocated between the two paths, which can be measured by load sharing coefficient. According to [9], the load sharing coefficient kls is expressed as(22)kls=max⁡F-L,F-RF0,where F-L, F-R are the average mesh forces of the first stage in left and right paths; F0 is the static mesh force of the first stage. The larger load sharing coefficient kls, the worth load sharing properties.It is mentioned in Section1 that there are two reasons which cause the unequal torques in the split torque transmission system: (1) the gap at one of the four gear mesh locations caused by manufacturing and assembly errors, which directly results in the difference of deformations between the two paths; (2) the unique phase difference properties of the split torque transmission system, which result in the desynchrony between the two paths. Since various split torque load sharing methods have been proposed to compensate for or minimize the gap, the effect of this gap caused by manufacturing and assembly errors will be ignored in this paper. The main factor which causes the unequal torques studied here is the phase difference only.According to [3], a split torque transmission design for helicopter is introduced as the reference design. The load sharing properties of the reference design will be studied first, and then an optimization will be conducted to promote its load sharing. The geometry parameters of the reference design have been listed in Table 1, other parameters are as follows: input power Pin is 372.85 kW, input speed nin is 8780 r/min, center distance H between input and output shaft is 293.45 mm, and material of gears is SAE 9310. An example of flank deviation surface for a teeth is shown in Figure 4, and the flank deviation is considered the sum of cumulative pitch error and profile error, where the cumulative pitch error is simulated by normally distributed constant and profile error is simulated by sine surface. The flank deviation of different teeth is simulated differently and the periodicity of the gear pair is taken into consideration.Table 1
The geometry parameters of the reference design.
Parameters
First stage
Second stage
Pinion/gear
Pinion/gear
Tooth number
32/124
27/176
Module/mm
1.590
2.54
Pressure angle/(°)
20
25
Helix angle/(°)
6
0
Tooth width/mm
44.45/38.10
66.04/59.94
Shaft angle/(°)
122.5
50.2
Phase difference
10.889
24.542Figure 4
An example of flank deviation surface for a tooth.The dynamic equations of the reference design system are established according to the previous section, and ode45 order is adopted to solve the functions with MATLAB software. The curves of mesh force and mesh stiffness of the first stage in both paths are evaluated at given power and speed, and parts of the curves are shown in Figures5 and 6.Figure 5
Mesh force curves of the first stage in both paths at given power and speed.Figure 6
Mesh stiffness curves of the first stage in both paths at given power and speed.As Figure5 illustrates, FL, FR are the mesh forces of left and right paths of the first stage, and F0 is the static mesh force of the first stage. It can be seen from the figure that the mesh force of the left path differs significantly from the right one, and the evaluated load sharing coefficient kls has reached 1.268.In Figure6, kL, kR are the mesh stiffness of left and right paths of the first stage, and inside the dashed rectangle box is the errorless mesh stiffness curve over one mesh period calculated by ISO6336. It can be seen that the actual mesh stiffness of both paths decrease sometimes, compared to the errorless mesh stiffness. It means the gear tooth of both paths does not always fully contact over the whole theoretical contact line. Therefore, contact length and mesh stiffness cannot be predicted before the nonlinear model is solved. Besides, the phase difference of mesh stiffness can be found easily in the figure.The load sharing coefficientkls under different power and speed is evaluated as is illustrated in Figure 7. It can be seen that the “load sharing map” (curved surface of load sharing coefficient kls under different operating conditions) is complicated, which comes from the nonlinearity of the split torque transmission system. The evaluated load sharing coefficient kls varies from 1.058 to 2.136, with input power varying from 223.7 kW to 522 kW and speed varying from 5268 r/min to 12292 r/min. The root mean square (RMS) of load sharing coefficient kls is 1.391. From a global point of view, the load sharing coefficient kls increases with the speed increasing and power decreasing, which correspond with [6]. The mesh load of the gear pair is so little under the condition of high speed and light power that the elastic deformation is not large enough to offset the initial deviations at all the contact points. However, it is noteworthy that the law of load sharing and operating conditions proposed here is not strictly correct, and there might be some counterexamples. It is because the fact that the actual length of contact line depends much on gear deviation under light load, which increases the nonlinearity of system dynamics.Figure 7
“Load sharing map” of the reference design.


## 4. Mathematical Model of Optimization

### 4.1. Objective Function

According to the previous section, the load sharing coefficient of a split torque transmission system changes greatly with different operating conditions, and the law of changing is complicated. Therefore, in order to obtain better load sharing from a system point of view, multiple operating conditions have to be taken into consideration. Considering the average case, the first objective function is promoted by minimizing the root mean square of load sharing coefficient under a wide range of operating conditions (possible operating conditions): input power varying from 223.7 kW (60% ofPin) to 522 kW (140% of Pin) and input speed varying from 5268 r/min (60% of nin) to 12292 r/min (140% of nin). When designing a gear transmission, light weight and safety are always important design targets. Safety is always measured by safety factors of contact fatigue strength and bending fatigue strength [24]. Therefore, the second and the third objective functions can be promoted by minimizing the total system mass and maximizing the total safety factors.To sum up, the whole objective functions are expressed as(23)min⁡Y1=kls,RMS,min⁡Y2=∑i=16Mgi,min⁡Y3=-∑S=-SH1+SH2+SF1+SF2,where kls,RMS is the RMS of load sharing coefficient (L-S-C RMS) under the possible operating conditions, Mgi is the mass of gear i, SH1, SH2 are the safety factor of contact fatigue strength of first and second stages, and SF1, SF2 are the safety factor of bending fatigue strength of first and second stages. Here only the mass of gears is considered in the total mass of the system.The safety factors of gear pairs can be evaluated according to ISO6336; it will not be discussed in detail here. However, the calculation of gear mass is a problem; for the method to calculate gear mass is associated with its wheel structure. Generally there are three types of wheel structure: solid type, panel type, and spoke type. The reason to adopt different types of wheel structure is reducing weight: as the gear gets larger the more percentage of mass is removed from the wheel. Here the light weight coefficientη is introduced to measure the extent of light weight, and then the gear mass Mgi can be calculated by(24)Mgi=ηiMgi∗i=1,2,…,6,where Mgi∗=ρ(π/4)bidi2 is the solid mass of gear i, with ρ: the material density, bi: the tooth width of gear i, and di: the pitch diameter of gear i.The introduction of light weight coefficientη unifies the different methods to calculate gear mass under different wheel structures, and the difference of three types of wheel structure is presented by varying the value of η. The type of wheel structure is decided by the tip diameter da, so the value of η is directly related to the tip diameter η=ηda. According to wheel structure design criteria, 3D parameterized model of a spur gear is created in CATIA V5 system. A series of gear models are created by varying the tip diameter da in a wide range and the masses of them are measured in CATIA V5 system, and then the values of η for gears with different tip diameters can be calculated. In the process of optimization, the value of η for a gear is obtained by interpolating according to its tip diameter.The values of the three objective functions of the reference design are evaluated: the root mean square of load sharing coefficientkls,RMS is 1.391, the system mass is 40.910 kg, and the total safety factors is 8.966.


### 4.2. Designing Variables

There are a lot of designing parameters in a split torque transmission system; some of them are independent while others are not. Picking up appropriate parameters as the designing variables is the prerequisite for optimization design. The special arrangement of the split torque transmission leads to the special mounting condition: the proportioning of gear tooth has definite interrelation with the two shaft angles [2]. Once the proportioning of gear tooth and the two shaft angles are determined, the center distances of the two stages are determined at the same time, which are restricted to the center distance between input and output shaft. Therefore, the modules of the two stages can hardly be the standard value. In other words, in order to guarantee the correct arrangement of a split torque transmission system, the standard of modules has to be sacrificed.Based on the above considerations, the designing variables selected here includes the gear ratio of the first stagei1, the pinion tooth number of the first and second stages Z1, Z4, the helix angle of the first stage βb1, and the two shaft angles Φ1, Φ2, as expressed in (25):(25)X=i1,Z1,Z4,βb1,Φ1,Φ2T.Other parameters can be evaluated by(26)i2=i0i1,Z2=round⁡i1·Z1,Z6=round⁡i2·Z4,m1=2Hcos⁡βb1Z1+Z2·sin⁡Φ2/2sin⁡Φ1/2+Φ2/2,m2=2H1Z4+Z6·sin⁡Φ1/2sin⁡Φ1/2+Φ2/2,where i0 is the total gear ratio, i2 is the gear ratio of the second stage, Z2, Z6 are the gear tooth numbers of the first and second stages, m1, m2 are the module of the first and second stages, and H is the center distance between input and output shafts. The function round· here means round · to the nearest integer.


### 4.3. Constraints

#### 4.3.1. Boundary Constraints

The design variables meet the following constraints:(27)i1,min≤i1≤i1,max,Z1,min≤Z1≤Z1,max,Z4,min≤Z4≤Z4,max,βb1,min≤βb1≤βb1,max,Φ1,min≤Φ1≤Φ1,max,Φ2,min≤Φ2≤Φ2,max,where i1, Z1, Z4, βb1, Φ1, Φ2 with the subscripts min and max are the boundaries of design variables, which are determined empirically according to the initial design.


#### 4.3.2. Performance Constraints

(i) Contact and bending fatigue strengths should be below the allowable values:(28)σH1,2≤σH1,2,σF1,2≤σF1,2,where σH1,2 is the contact stress of the first and second stages, [σH1,2] is the allowable contact stress of the first and second stages, σF1,2 is the bending stress of the first and second stages, and [σF1,2] is the allowable bending stress of the first and second stages.(ii) First stage gears should not interfere with each other:(29)m1Z2+2ha∗<m1Z1+Z2cos⁡βb1·sin⁡Φ12,where ha∗ is the addendum factor of the first stage gear.(iii) Safety margin of each stage should be balanced:(30)14∑i=14ΔSi-ΔS¯2<ε,where (31)ΔSi=SH,i-SHmini=1,2SF,i-2-SFmini=3,4is the safety margin, with SHmin: the minimum safety factor for contact fatigue strength and SFmin: the minimum safety factor for bending fatigue strength; ΔS¯ is the mean value of safety margin ΔS; ε is the allowable upper limit of safety margin standard deviation.


## 4.1. Objective Function

According to the previous section, the load sharing coefficient of a split torque transmission system changes greatly with different operating conditions, and the law of changing is complicated. Therefore, in order to obtain better load sharing from a system point of view, multiple operating conditions have to be taken into consideration. Considering the average case, the first objective function is promoted by minimizing the root mean square of load sharing coefficient under a wide range of operating conditions (possible operating conditions): input power varying from 223.7 kW (60% ofPin) to 522 kW (140% of Pin) and input speed varying from 5268 r/min (60% of nin) to 12292 r/min (140% of nin). When designing a gear transmission, light weight and safety are always important design targets. Safety is always measured by safety factors of contact fatigue strength and bending fatigue strength [24]. Therefore, the second and the third objective functions can be promoted by minimizing the total system mass and maximizing the total safety factors.To sum up, the whole objective functions are expressed as(23)min⁡Y1=kls,RMS,min⁡Y2=∑i=16Mgi,min⁡Y3=-∑S=-SH1+SH2+SF1+SF2,where kls,RMS is the RMS of load sharing coefficient (L-S-C RMS) under the possible operating conditions, Mgi is the mass of gear i, SH1, SH2 are the safety factor of contact fatigue strength of first and second stages, and SF1, SF2 are the safety factor of bending fatigue strength of first and second stages. Here only the mass of gears is considered in the total mass of the system.The safety factors of gear pairs can be evaluated according to ISO6336; it will not be discussed in detail here. However, the calculation of gear mass is a problem; for the method to calculate gear mass is associated with its wheel structure. Generally there are three types of wheel structure: solid type, panel type, and spoke type. The reason to adopt different types of wheel structure is reducing weight: as the gear gets larger the more percentage of mass is removed from the wheel. Here the light weight coefficientη is introduced to measure the extent of light weight, and then the gear mass Mgi can be calculated by(24)Mgi=ηiMgi∗i=1,2,…,6,where Mgi∗=ρ(π/4)bidi2 is the solid mass of gear i, with ρ: the material density, bi: the tooth width of gear i, and di: the pitch diameter of gear i.The introduction of light weight coefficientη unifies the different methods to calculate gear mass under different wheel structures, and the difference of three types of wheel structure is presented by varying the value of η. The type of wheel structure is decided by the tip diameter da, so the value of η is directly related to the tip diameter η=ηda. According to wheel structure design criteria, 3D parameterized model of a spur gear is created in CATIA V5 system. A series of gear models are created by varying the tip diameter da in a wide range and the masses of them are measured in CATIA V5 system, and then the values of η for gears with different tip diameters can be calculated. In the process of optimization, the value of η for a gear is obtained by interpolating according to its tip diameter.The values of the three objective functions of the reference design are evaluated: the root mean square of load sharing coefficientkls,RMS is 1.391, the system mass is 40.910 kg, and the total safety factors is 8.966.


## 4.2. Designing Variables

There are a lot of designing parameters in a split torque transmission system; some of them are independent while others are not. Picking up appropriate parameters as the designing variables is the prerequisite for optimization design. The special arrangement of the split torque transmission leads to the special mounting condition: the proportioning of gear tooth has definite interrelation with the two shaft angles [2]. Once the proportioning of gear tooth and the two shaft angles are determined, the center distances of the two stages are determined at the same time, which are restricted to the center distance between input and output shaft. Therefore, the modules of the two stages can hardly be the standard value. In other words, in order to guarantee the correct arrangement of a split torque transmission system, the standard of modules has to be sacrificed.Based on the above considerations, the designing variables selected here includes the gear ratio of the first stagei1, the pinion tooth number of the first and second stages Z1, Z4, the helix angle of the first stage βb1, and the two shaft angles Φ1, Φ2, as expressed in (25):(25)X=i1,Z1,Z4,βb1,Φ1,Φ2T.Other parameters can be evaluated by(26)i2=i0i1,Z2=round⁡i1·Z1,Z6=round⁡i2·Z4,m1=2Hcos⁡βb1Z1+Z2·sin⁡Φ2/2sin⁡Φ1/2+Φ2/2,m2=2H1Z4+Z6·sin⁡Φ1/2sin⁡Φ1/2+Φ2/2,where i0 is the total gear ratio, i2 is the gear ratio of the second stage, Z2, Z6 are the gear tooth numbers of the first and second stages, m1, m2 are the module of the first and second stages, and H is the center distance between input and output shafts. The function round· here means round · to the nearest integer.


## 4.3. Constraints

### 4.3.1. Boundary Constraints

The design variables meet the following constraints:(27)i1,min≤i1≤i1,max,Z1,min≤Z1≤Z1,max,Z4,min≤Z4≤Z4,max,βb1,min≤βb1≤βb1,max,Φ1,min≤Φ1≤Φ1,max,Φ2,min≤Φ2≤Φ2,max,where i1, Z1, Z4, βb1, Φ1, Φ2 with the subscripts min and max are the boundaries of design variables, which are determined empirically according to the initial design.


### 4.3.2. Performance Constraints

(i) Contact and bending fatigue strengths should be below the allowable values:(28)σH1,2≤σH1,2,σF1,2≤σF1,2,where σH1,2 is the contact stress of the first and second stages, [σH1,2] is the allowable contact stress of the first and second stages, σF1,2 is the bending stress of the first and second stages, and [σF1,2] is the allowable bending stress of the first and second stages.(ii) First stage gears should not interfere with each other:(29)m1Z2+2ha∗<m1Z1+Z2cos⁡βb1·sin⁡Φ12,where ha∗ is the addendum factor of the first stage gear.(iii) Safety margin of each stage should be balanced:(30)14∑i=14ΔSi-ΔS¯2<ε,where (31)ΔSi=SH,i-SHmini=1,2SF,i-2-SFmini=3,4is the safety margin, with SHmin: the minimum safety factor for contact fatigue strength and SFmin: the minimum safety factor for bending fatigue strength; ΔS¯ is the mean value of safety margin ΔS; ε is the allowable upper limit of safety margin standard deviation.


## 4.3.1. Boundary Constraints

The design variables meet the following constraints:(27)i1,min≤i1≤i1,max,Z1,min≤Z1≤Z1,max,Z4,min≤Z4≤Z4,max,βb1,min≤βb1≤βb1,max,Φ1,min≤Φ1≤Φ1,max,Φ2,min≤Φ2≤Φ2,max,where i1, Z1, Z4, βb1, Φ1, Φ2 with the subscripts min and max are the boundaries of design variables, which are determined empirically according to the initial design.


## 4.3.2. Performance Constraints

(i) Contact and bending fatigue strengths should be below the allowable values:(28)σH1,2≤σH1,2,σF1,2≤σF1,2,where σH1,2 is the contact stress of the first and second stages, [σH1,2] is the allowable contact stress of the first and second stages, σF1,2 is the bending stress of the first and second stages, and [σF1,2] is the allowable bending stress of the first and second stages.(ii) First stage gears should not interfere with each other:(29)m1Z2+2ha∗<m1Z1+Z2cos⁡βb1·sin⁡Φ12,where ha∗ is the addendum factor of the first stage gear.(iii) Safety margin of each stage should be balanced:(30)14∑i=14ΔSi-ΔS¯2<ε,where (31)ΔSi=SH,i-SHmini=1,2SF,i-2-SFmini=3,4is the safety margin, with SHmin: the minimum safety factor for contact fatigue strength and SFmin: the minimum safety factor for bending fatigue strength; ΔS¯ is the mean value of safety margin ΔS; ε is the allowable upper limit of safety margin standard deviation.


## 5. Algorithm and Improvement

Classical optimization methods suggest converting the multiobjective optimization problem to a single-objective optimization problem by weighted sum of all the objectives, which lead to disadvantage of subjectivity when determining the weights of objectives. Whereas the improved nondominated sorting genetic algorithm (NSGA-II) proposed by Deb avoids this disadvantage. In this paper, NSGA-II is adopted to solve the proposed multiobjective optimization model.However, there are large numbers of nonlinear dynamic equations to be solved under multiple operating conditions when evaluating the fitness, and the solving time can be too long to accept. Therefore, an improved NSGA-II algorithm is being put forward to solve the problem of time consuming: prediction strategy is used in the fitness evaluation step so as to avoid the evaluation of load sharing property which is computationally very expensive. The key of the improvement is predicting instead of evaluating the real fitness.


### 5.1. The Fitness Prediction Strategy

As is known to all, NSGA-II is a population based evolutionary algorithm for multiobjective optimization problems and the population evolves with the generation increasing. In the improved NSGA-II algorithm, each individuali in the population has its fitness vector: (32)ftni=ftn1i,ftn2i,…,ftnNiand the fidelity Ri of the fitness vector, where N is the number of objectives. The value of the fitness vector can be evaluated using the fitness functions or predicted through the values of other individuals. If the fitness vector ftni is evaluated using the original fitness functions, the fidelity Ri=1; if ftni is estimated, the fidelity 0≤Ri<1.As shown in Figure8, for each individual i, specify its fitness sharing radius rs. The area in which the dimensionless Euclidean distance between individual i and any other one is no greater than the fitness sharing radius is called fitness sharing area for individual i, expressed as Ωi. Assume there are m other individuals in the fitness sharing area Ωi, which composed a collection S=s1,s2,…,sm. And the evaluation method of ftni is as follows.Figure 8
Fitness sharing area.Firstly, evaluate the fidelityRi of individual i:(33)Ri=∑j=1mωsj,i·Rsj,where sj is an individual including fitness sharing area Ωi, Rsj is the fidelity of sj, ωsj,i is the weight of sj donated by individual i. Let dimensionless Euclidean distance between individuals s1,s2,…,sm and individual i be d¯1i,d¯2i,…,d¯mi, respectively; then ωsj,i can be evaluated by(34)ωsj,i=exp⁡-γ·d¯ji∑k=1mexp⁡-γ·d¯ki,j=1,2,…,m,where γ is weight rescaling factor. The closer the individual is to individual i, the greater contribution of fidelity it makes.As depicted in Figure9, if fidelity Ri is greater than a given threshold R∗, then the fitness vector ftni can be predicted as(35)ftni=∑j=1mωsj,i·ftnsj.Else if Ri is less than the threshold R∗, then the fitness vector ftni should be evaluated using its original fitness functions.Figure 9
Fitness prediction model.To take full advantage of the historical population whose fitness has been evaluated, it is necessary to establish a database of historical populations of coordinates, fitness, and fidelities. As the population evolves, the database will expand the scale gradually. In order to reduce space complexity and the amount of computation, redundant data need to be eliminated after each generation. The concept of individual redundancy is introduced to determine whether the data is redundant, which is defined as(36)Iri=∑k=1nΔx¯ki,where Δx¯ki is the coordinate difference (absolute value) of individual i between the former individual and the latter one in the kth dimensional and n is the number of design parameters. If the redundancy value of an individual is less than a given threshold, the individual is knocked out.In addition, since not all of the individuals’ fitness vectors are evaluated using its original fitness functions, the predicted fitness vectors are not accurate. Thus, as the population evolves gradually, the fidelities of individuals with predicted fitness vectors should decline gradually. Assume that individuali is with the predicted fitness vectors, let the fidelity of individual i in generation t be Ri,t, and then the fidelity in generation t+1 can be updated as(37)Ri,t+1=β·Ri,t,where β is fidelity drain factor, with 0<β<1. As the population evolves, the fidelity drops below a given threshold R0, and the individual should also be removed from the database.


### 5.2. Algorithm Flow

Step 1.
Initialize historical population database, set the initial population blank, and set fitness vectors and fidelities to 0.Step 2.
Find the fitness sharing area for each individuali, and find the collection of individuals in the area from the database.Step 3.
Evaluate the fidelityRi of individual i, and determine whether Ri is greater than the threshold R∗. If Ri≥R∗, predict fitness vector of individual i according to (35); otherwise, evaluate the fitness vector using its original functions and set the fidelity Ri to 1.Step 4.
Add individuali to the database.Step 5.
Update the database as follows: (1) calculate redundancy for all individuals and eliminate all redundant individuals; (2) for all individuals with predicted fitness vectors, update its fidelity according to (37) and remove all the individuals with low fidelity.


### 5.3. Numerical Experiments

There are two purposes of conducting numerical experiments on the improved NSGA-II (iNSGA-II): (1) testing convergence of the algorithm, which tests whether the algorithm can correctly guide the evolution, so that the Pareto optimal solution (Pareto front) of the original problem can be obtained; (2) testing the effectiveness of the algorithm, namely, testing what extent the algorithm can reduce the amount of computation.Four benchmark problems are chosen from a number of significant past studies in multiobjective optimization area: Schaffer’s study (SCH) [25], Fonseca and Fleming’s study (FON) [26], Poloni’s study (POL) [27], and Kursawe’s study (KUR) [28]. The benchmark problems are described as follows.(1) SCH Problem (
n
=
1
)
(38)
min
⁡
f
1
x
=
x
2
,
min
⁡
f
2
x
=
x
-
2
2
,where the variable bound is -103,103 and the Pareto optimal front is convex.(2) FON Problem (
n
=
3
)
(39)
min
⁡
f
1
x
=
1
-
exp
⁡
-
∑
i
=
1
3
x
i
-
1
3
2
,
min
⁡
f
2
x
=
1
-
exp
⁡
-
∑
i
=
1
3
x
i
+
1
3
2
,where the variable bounds are -4,4 and the Pareto optimal front is nonconvex.(3) POL Problem (
n
=
2
)
(40)
min
⁡
f
1
x
=
1
+
A
1
-
B
1
2
+
A
2
-
B
2
2
,
min
⁡
f
2
x
=
x
1
+
3
2
+
x
2
+
1
2
,where(41)A1=0.5sin⁡1-2cos⁡1+sin⁡2-1.5cos⁡2,A2=1.5sin⁡1-cos⁡1+2sin⁡2-0.5cos⁡2,B1=0.5sin⁡x1-2cos⁡x1+sin⁡x2-1.5cos⁡x2,B2=1.5sin⁡x1-cos⁡x1+2sin⁡x2-0.5cos⁡x2,and the variable bounds are -π,π and the Pareto optimal front is nonconvex and disconnected.(4) KUR Problem (
n
=
3
)
(42)
min
⁡
f
1
x
=
∑
i
=
1
n
-
1
-
10
exp
⁡
-
0.2
x
i
2
+
x
i
+
1
2
,
min
⁡
f
2
x
=
∑
i
=
1
n
x
i
0.8
+
5
sin
⁡
x
i
3
,where the variable bounds are -5,5 and the Pareto optimal front is nonconvex and disconnected.The four benchmark problems are solved by iNSGA-II with MATLAB programming. Binary-coding, single-point crossover, and bitwise mutation are used in the algorithm. The algorithm parameters are settled as follows: population size is 100, evolution generation is 200, the crossover probability is 0.9, the mutation probability is 0.1, thresholdR∗=0.6, and fidelity drain factor β=0.9. Each problem is tested 20 times, respectively.Pareto optimal fronts of the four benchmark problems with iNSGA-II are illustrated in Figure10, where (a), (b), (c), and (d) represent SCH problem, FON problem, POL problem, and KUR problem, respectively. It can be seen that the improved NSGA-II algorithm (iNSGA-II) achieves Pareto fronts correctly in the four benchmark problems.Figure 10
Pareto optimal fronts of the four benchmark problems with iNSGA-II.
(a)
SCH (convex)
(b)
FON (nonconvex)
(c)
POL (nonconvex and disconnected)
(d)
KUR (nonconvex and disconnected)Percentage of real fitness vectors evaluated in each generation of the four benchmark problems are represented in Figure11, where (a), (b), (c), and (d) represent SCH problem, FON problem, POL problem, and KUR problem, respectively. Percentage of real fitness evaluated in each generation substantially stabilized at 30% to 50%. A certain percentage of the individuals’ fitness vectors is evaluated using the original fitness functions in each generation, so that the evolutionary direction can be guided correctly.Figure 11
Percentage of real fitness evaluated.
(a)
SCH
(b)
FON
(c)
POL
(d)
KURThe detailed testing results of effectiveness are illustrated in Table2, where “Max.Num” represents the maximum number of fitness vectors to be evaluated or predicted, “Eval.Num” represents the average number of evaluated fitness vectors, and “Percentage” represents the percentage of “Max.Num” and “Eval.Num.” It can be known from Table 2 that iNSGA-II can reduce much computation amount of real fitness compared with traditional NSGA-II. It means that, when the real fitness evaluation is computationally very expensive, using iNSGA-II can save approximately 2/3 of the computing time.Table 2
Effectiveness test results of the algorithm.
Problems
Max.Num
Eval.Num
Percentage
SCH
20000
7167.5
35.838%
FON
20000
7560.5
37.803%
POL
20000
7382.5
36.913%
KUR
20000
7212.05
36.060%


## 5.1. The Fitness Prediction Strategy

As is known to all, NSGA-II is a population based evolutionary algorithm for multiobjective optimization problems and the population evolves with the generation increasing. In the improved NSGA-II algorithm, each individuali in the population has its fitness vector: (32)ftni=ftn1i,ftn2i,…,ftnNiand the fidelity Ri of the fitness vector, where N is the number of objectives. The value of the fitness vector can be evaluated using the fitness functions or predicted through the values of other individuals. If the fitness vector ftni is evaluated using the original fitness functions, the fidelity Ri=1; if ftni is estimated, the fidelity 0≤Ri<1.As shown in Figure8, for each individual i, specify its fitness sharing radius rs. The area in which the dimensionless Euclidean distance between individual i and any other one is no greater than the fitness sharing radius is called fitness sharing area for individual i, expressed as Ωi. Assume there are m other individuals in the fitness sharing area Ωi, which composed a collection S=s1,s2,…,sm. And the evaluation method of ftni is as follows.Figure 8
Fitness sharing area.Firstly, evaluate the fidelityRi of individual i:(33)Ri=∑j=1mωsj,i·Rsj,where sj is an individual including fitness sharing area Ωi, Rsj is the fidelity of sj, ωsj,i is the weight of sj donated by individual i. Let dimensionless Euclidean distance between individuals s1,s2,…,sm and individual i be d¯1i,d¯2i,…,d¯mi, respectively; then ωsj,i can be evaluated by(34)ωsj,i=exp⁡-γ·d¯ji∑k=1mexp⁡-γ·d¯ki,j=1,2,…,m,where γ is weight rescaling factor. The closer the individual is to individual i, the greater contribution of fidelity it makes.As depicted in Figure9, if fidelity Ri is greater than a given threshold R∗, then the fitness vector ftni can be predicted as(35)ftni=∑j=1mωsj,i·ftnsj.Else if Ri is less than the threshold R∗, then the fitness vector ftni should be evaluated using its original fitness functions.Figure 9
Fitness prediction model.To take full advantage of the historical population whose fitness has been evaluated, it is necessary to establish a database of historical populations of coordinates, fitness, and fidelities. As the population evolves, the database will expand the scale gradually. In order to reduce space complexity and the amount of computation, redundant data need to be eliminated after each generation. The concept of individual redundancy is introduced to determine whether the data is redundant, which is defined as(36)Iri=∑k=1nΔx¯ki,where Δx¯ki is the coordinate difference (absolute value) of individual i between the former individual and the latter one in the kth dimensional and n is the number of design parameters. If the redundancy value of an individual is less than a given threshold, the individual is knocked out.In addition, since not all of the individuals’ fitness vectors are evaluated using its original fitness functions, the predicted fitness vectors are not accurate. Thus, as the population evolves gradually, the fidelities of individuals with predicted fitness vectors should decline gradually. Assume that individuali is with the predicted fitness vectors, let the fidelity of individual i in generation t be Ri,t, and then the fidelity in generation t+1 can be updated as(37)Ri,t+1=β·Ri,t,where β is fidelity drain factor, with 0<β<1. As the population evolves, the fidelity drops below a given threshold R0, and the individual should also be removed from the database.


## 5.2. Algorithm Flow

Step 1.
Initialize historical population database, set the initial population blank, and set fitness vectors and fidelities to 0.Step 2.
Find the fitness sharing area for each individuali, and find the collection of individuals in the area from the database.Step 3.
Evaluate the fidelityRi of individual i, and determine whether Ri is greater than the threshold R∗. If Ri≥R∗, predict fitness vector of individual i according to (35); otherwise, evaluate the fitness vector using its original functions and set the fidelity Ri to 1.Step 4.
Add individuali to the database.Step 5.
Update the database as follows: (1) calculate redundancy for all individuals and eliminate all redundant individuals; (2) for all individuals with predicted fitness vectors, update its fidelity according to (37) and remove all the individuals with low fidelity.


## 5.3. Numerical Experiments

There are two purposes of conducting numerical experiments on the improved NSGA-II (iNSGA-II): (1) testing convergence of the algorithm, which tests whether the algorithm can correctly guide the evolution, so that the Pareto optimal solution (Pareto front) of the original problem can be obtained; (2) testing the effectiveness of the algorithm, namely, testing what extent the algorithm can reduce the amount of computation.Four benchmark problems are chosen from a number of significant past studies in multiobjective optimization area: Schaffer’s study (SCH) [25], Fonseca and Fleming’s study (FON) [26], Poloni’s study (POL) [27], and Kursawe’s study (KUR) [28]. The benchmark problems are described as follows.(1) SCH Problem (
n
=
1
)
(38)
min
⁡
f
1
x
=
x
2
,
min
⁡
f
2
x
=
x
-
2
2
,where the variable bound is -103,103 and the Pareto optimal front is convex.(2) FON Problem (
n
=
3
)
(39)
min
⁡
f
1
x
=
1
-
exp
⁡
-
∑
i
=
1
3
x
i
-
1
3
2
,
min
⁡
f
2
x
=
1
-
exp
⁡
-
∑
i
=
1
3
x
i
+
1
3
2
,where the variable bounds are -4,4 and the Pareto optimal front is nonconvex.(3) POL Problem (
n
=
2
)
(40)
min
⁡
f
1
x
=
1
+
A
1
-
B
1
2
+
A
2
-
B
2
2
,
min
⁡
f
2
x
=
x
1
+
3
2
+
x
2
+
1
2
,where(41)A1=0.5sin⁡1-2cos⁡1+sin⁡2-1.5cos⁡2,A2=1.5sin⁡1-cos⁡1+2sin⁡2-0.5cos⁡2,B1=0.5sin⁡x1-2cos⁡x1+sin⁡x2-1.5cos⁡x2,B2=1.5sin⁡x1-cos⁡x1+2sin⁡x2-0.5cos⁡x2,and the variable bounds are -π,π and the Pareto optimal front is nonconvex and disconnected.(4) KUR Problem (
n
=
3
)
(42)
min
⁡
f
1
x
=
∑
i
=
1
n
-
1
-
10
exp
⁡
-
0.2
x
i
2
+
x
i
+
1
2
,
min
⁡
f
2
x
=
∑
i
=
1
n
x
i
0.8
+
5
sin
⁡
x
i
3
,where the variable bounds are -5,5 and the Pareto optimal front is nonconvex and disconnected.The four benchmark problems are solved by iNSGA-II with MATLAB programming. Binary-coding, single-point crossover, and bitwise mutation are used in the algorithm. The algorithm parameters are settled as follows: population size is 100, evolution generation is 200, the crossover probability is 0.9, the mutation probability is 0.1, thresholdR∗=0.6, and fidelity drain factor β=0.9. Each problem is tested 20 times, respectively.Pareto optimal fronts of the four benchmark problems with iNSGA-II are illustrated in Figure10, where (a), (b), (c), and (d) represent SCH problem, FON problem, POL problem, and KUR problem, respectively. It can be seen that the improved NSGA-II algorithm (iNSGA-II) achieves Pareto fronts correctly in the four benchmark problems.Figure 10
Pareto optimal fronts of the four benchmark problems with iNSGA-II.
(a)
SCH (convex)
(b)
FON (nonconvex)
(c)
POL (nonconvex and disconnected)
(d)
KUR (nonconvex and disconnected)Percentage of real fitness vectors evaluated in each generation of the four benchmark problems are represented in Figure11, where (a), (b), (c), and (d) represent SCH problem, FON problem, POL problem, and KUR problem, respectively. Percentage of real fitness evaluated in each generation substantially stabilized at 30% to 50%. A certain percentage of the individuals’ fitness vectors is evaluated using the original fitness functions in each generation, so that the evolutionary direction can be guided correctly.Figure 11
Percentage of real fitness evaluated.
(a)
SCH
(b)
FON
(c)
POL
(d)
KURThe detailed testing results of effectiveness are illustrated in Table2, where “Max.Num” represents the maximum number of fitness vectors to be evaluated or predicted, “Eval.Num” represents the average number of evaluated fitness vectors, and “Percentage” represents the percentage of “Max.Num” and “Eval.Num.” It can be known from Table 2 that iNSGA-II can reduce much computation amount of real fitness compared with traditional NSGA-II. It means that, when the real fitness evaluation is computationally very expensive, using iNSGA-II can save approximately 2/3 of the computing time.Table 2
Effectiveness test results of the algorithm.
Problems
Max.Num
Eval.Num
Percentage
SCH
20000
7167.5
35.838%
FON
20000
7560.5
37.803%
POL
20000
7382.5
36.913%
KUR
20000
7212.05
36.060%


## 6. Results

The improved NSGA-II algorithm for solving the multiobjective optimization model is realized by MATLAB programming. As is known to all, there are a set of optimal solutions (largely known as Pareto optimal solutions) in a multiobjective optimization problem, instead of a single optimal solution. The Pareto optimal solutions form a Pareto optimal front, which has the property that one solution in the Pareto optimal front cannot be said to be better than any of the others [21]. Therefore, when the Pareto optimal front is obtained, designers can choose any solution in it according to different requirements and favors. The Pareto optimal front for two objective problems is always a continuous curve, whereas for three objective problems it is a curved surface.Pareto optimal front (blue points) for the three objective problems of split torque transmission is illustrated in Figure12. In order to observe the Pareto optimal front conveniently, spatial curved surface fitting for the Pareto optimal solution points is applied. As shown in Figure 12, coordinate Y1 is the L-S-C RMS kls,RMS, coordinate Y2 is the total mass ∑M of the system, and coordinate Y3 is the opposite number of total safety factors ∑S.Figure 12
The Pareto optimal front.In order to choose an optimal design conveniently from the Pareto optimal front, the curved surface of the front is projected onto a 2D plane, as shown in Figure13. The x-coordinate is the total mass of the system ∑M, the y-coordinate is the opposite number of total safety factors ∑S, and the information of L-S-C RMS kls,RMS is expressed by color. As Figure 13 illustrates, the Pareto optimal front is divided into 10 belt areas and L-S-C RMS kls,RMS of each area decreases from area (a) to area (j). The distributions of values of the three objective functions in the Pareto optimal front are as follows: the L-S-C RMS kls,RMS varies from 1.166 to 1.709, the total mass of the system ∑M varies from 33.513 kg to 50.448 kg, and the total safety factors ∑S varies from 7.796 to 14.450. As is mentioned above, one solution in the Pareto optimal front cannot be said to be better than any of the others; therefore, any solution can be chosen according to designers’ requirements and favors. If the load sharing property is biased, then areas (i) and (j) will be suitable; if the lightweight is biased, then areas (e) and (f) will be selected; if the safety is biased, then areas (a) and (b) will be the answer. Besides, designer can also judge and weigh among the three objectives and then make an optimal choice.Figure 13
Projected Pareto optimal front.To promise that the optimal design is better than the referance design in all the three objectives, all the unqualified solutions are eliminated, and then the left solutions can be further divided based on RMS of load sharing coefficient. As Figure14 illustrates, there are still 8 optimal solutions which are better than the referance one in all the three objectives. Four alternative solutions 1#~4# are selected, with 1#: the solution with minimum L-S-C RMS, 2#: the solution with maximum total safety factor, 3#: the solution with minimum total mass of the system, and 4#: the optimal solution which weighs among the three objectives. The values of design parameters and objective functions of solutions 1#~4# are listed in Table 3.Table 3
The values of design parameters and objective functions of solutions 1#~4#.
Parameters
1#
2#
3#
4#
First stage
Tooth numberZ1/Z2
28/119
26/91
26/105
26/105
Modulem1/mm
1.538
1.861
1.683
1.694
Gear ratioi1
4.25
3.5
4.038
4.038
Helix angleβb1/(°)
18.900
14.296
18.196
17.402
Second stage
Tooth numberZ4/Z6
22/131
22/159
22/138
22/138
Modulem2/mm
3.270
2.981
3.095
3.119
Gear ratioi2
5.955
7.227
6.273
6.273
Shaft angleΦ1/Φ2 (°)
114.741/47.439
133.641/45.018
111.615/45.617
113.589/45.899
Phase differenceΔφ1/Δφ2
16.558/8.411
19.295/9.145
16.698/7.559
16.839/7.709
L-S-C RMSkls,RMS
1.287
1.381
1.317
1.318
Total mass∑M/kg
39.066
40.516
36.338
37.088
Total safety factor∑S
9.279
10.582
9.253
9.367Figure 14
Further divided optimal solutions.Solution 3# is selected as the final optimal design through balancing, with L-S-C RMS deceasing from 1.391 to 1.317, the system mass deceasing from 40.910 kg to 36.338 kg, and the total safety factors increasing from 8.966 to 9.253. Figure15 illustrates the “load sharing map” of solution 3# (the colorized one), and it can be seen that the nonlinearity still exists in the optimal design. The load sharing coefficient kls of solution 3# under possible operating conditions varies from 1.032 to 2.023, which under the normal work condition is 1.083. Compared with the “load sharing map” of the reference design (the gray translucent one), the optimal design achieves better load sharing property under wide range of operating conditions.Figure 15
Load sharing map of 3# and the reference design.


## 7. Conclusions

Gap and phase differences are two of the mainly causes of unequaled torques in a split torque transmission system. Various methods have been proposed to compensate for or minimize the effect of gap, promoting the load sharing property to a certain extent. In this paper, system design parameters were optimized to change the phase difference, thereby further improving the load sharing property. Therefore, a multiobjective optimization design of a split torque transmission system was conducted, with the promoting of load sharing property, lightweight, and safety considered as the objectives. The load sharing property, which was measured by load sharing coefficient, was evaluated under multiple operating conditions with dynamic analysis method. An improved NSGA-II algorithm was adopted to solve the multiobjective model and a satisfied optimal solution was picked up as the final optimal design from the Pareto optimal front.The main contributions of this paper include the following.(1)
A new method for promoting load sharing property of a split torque transmission system is proposed, which changes the phase difference by optimizing the system design parameters. This method is proved effective by the results with the reference design and can be used with other load sharing methods together.(2)
The concept of “load sharing map” is put forward to measure the load sharing property under multiple operating conditions, which allows us to find an optimized design over a wide range of operating conditions.(3)
An improved NSGA-II is proposed to overcome the problem of time consuming during the process of fitness evaluation, which promotes the effectiveness of optimization with nearly 2/3 of the total computations reduced.



---

*Source: 381010-2015-10-12.xml*