Load Sharing Multiobjective Optimization Design of a Split Torque Helicopter Transmission

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.


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.

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 Figure 1.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 , , and ) and 1 rotational degree (around axis ) each.Here, the angle  between axis  and the center line of the gear pair   →     is defined as direction angle.
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 Figure 1 illustrates, when a gear pair is nominally engaged, contact occurs at one certain point  * , which has the maximum deviation, namely, ( * ).Select any point  in the theoretical contact line different from  * , and name the flank deviations of  in the pinion and gear as   () and   (), 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  can be expressed as The gear model is shown in Figure 1, and pinion and gear are assimilated to rigid cylinders with four-degreesof-freedom connected by series of stiffness and damp.The generalized displacements can be expressed as where   ,   ,   are translational displacements of the pinion along axes , , , respectively.  is rotational displacement of the pinion around axis ;   ,   ,   are translational displacements of the gear along axes , , , respectively.  is rotational displacement of the gear around axis .
Since generalized displacements are small quantities, it can therefore be represented by infinitesimal translations and rotations.The perturbation q on pinion and gear gives a normal approach () on the base plane relative to rigid body positions: where  is a projective vector depending on gear geometry which projects displacements q onto the base plane;   is the base helix angle (the value is 0 to spur gears);  is the pressure angle;   ,   are the base radius of pinion and gear.
Contact and deflection in  occur only if the normal approach () is larger than the initial deviation (), and in such case the deflection Δ() can be evaluated as where () = ( * ) − (), ( * ), () is the deviation in  * , .Else if the normal approach () is less than the initial deviation (), then the deflection Δ() = 0.The elemental mesh force transmitted from the pinion onto the gear at one point of contact  is where () is the mesh stiffness at point  per unit of contact length which is calculated according to ISO6336 and  is the elemental contact length.The total mesh force   can be deducted by integrating over the time-varying and deflection dependent contact length:

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.
As shown in Figure 2, 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  and , with  to the right of .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 Figure 3, a right-hand Cartesian coordinate system is established such that the -axis is coincident with the output gear shaft, the positive -axis extends from the input gear center to the output pinion center, and the input gear drives clockwise.The first-stage pinion, gear (), and gear () are marked with gears 1∼3; the second-stage pinion (), pinion (), and gear are marked with gears 4∼6.The input shaft, compound shaft (), compound shaft (), and output shaft are marked with axes 1∼4, and the center of which is marked with  1 ,  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 points ,  (the length of arc  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 where Δ 1,2 is the phase difference of the first and second stages;  ,   are the length of arcs , ;  1 ,  2 are the circular pitch of the first and second stages;  1 ,  2 are the transverse module of the first and second stages;   is the pitch diameter of gear ;   is the tooth number of gear .
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 where the subscripts 1∼6 correspond to gears 1∼6.

Equations of Dynamics. The general form of dynamic equations of a gear transmission system is
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 where K  is the mesh stiffness term, K  is the supporting stiffness term, K  is the coupling stiffness term, F 0 is the constant torques term, and F() 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.

Mass Matrix M.
Mass matrix M is a 24 × 24 diagonal matrix, which is expressed as where   ,  = 1 ∼ 6 is the mass of gears 1∼6,   ,  = 1 ∼ 6 is the rotational inertia around -axis.

Stiffness Matrix K.
The stiffness matrix K consists of mesh stiffness K  , supporting stiffness K  , and coupling stiffness K  , and it is also a 24 × 24 symmetric matrix: (i) Meshing Stiffness Matrix K  .The mesh stiffness matrix of the whole system K  comes from mesh stiffness of all the gear pairs; it can therefore be obtained by assembling the entire mesh stiffness submatrix together: where K , = ∫    () ⋅      is the nonlinear and timedependent mesh stiffness submatrix of gear pair , , with   : the projective vector of gear pair  and : is assemble matrix, with R 1 = R 2 = diag(1, 1, 1, 1).
(ii) Supporting Stiffness Matrix K  .The supporting stiffness matrix K  comes from all the gears' supporting stiffness in all the translational motion freedoms: where   ,   ,   ,  = 1 ∼ 6 is the supporting stiffness of gear  in , ,  directions.
(iii) Coupling Stiffness Matrix K  .The coupling stiffness matrix K  of the whole system comes from the two compound shafts, and it can be obtained by assembling the two coupling stiffness submatrix together: where is the coupling stiffness submatrix of compound gear , , with  , ,  , ,  , : the torsional stiffness, bending stiffness, axial stiffness between the compound gears  and .

Damp Matrix C.
Here Rayleigh's damping is adopted: with , : two constants to be adjusted from experimental results and experience.

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   is expressed as where   ,   are the average mesh forces of the first stage in left and right paths;  0 is the static mesh force of the first stage.The larger load sharing coefficient   , the worth load sharing properties.
It is mentioned in Section 1 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  in is 372.85 kW, input speed  in is 8780 r/min, center distance  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.
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 Figures 5 and 6.
As Figure 5 illustrates,   ,   are the mesh forces of left and right paths of the first stage, and  0 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   has reached 1.268.
In Figure 6,   ,   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 coefficient   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   under different operating conditions) is complicated, which comes from the nonlinearity of the split torque transmission system.The evaluated load sharing coefficient   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   is 1.391.From a global point of view, the load sharing coefficient   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.

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% of  in ) to 522 kW (140% of  in ) and input speed varying from 5268 r/min (60% of  in ) to 12292 r/min (140% of  in ).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 min  1 =  ,RMS , where  ,RMS is the RMS of load sharing coefficient (L-S-C RMS) under the possible operating conditions,   is the mass of gear ,  1 ,  2 are the safety factor of contact fatigue strength of first and second stages, and  1 ,  2 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   can be calculated by where  *  = (/4)    2 is the solid mass of gear , with : the material density,   : the tooth width of gear , and   : the pitch diameter of gear .

Mathematical Problems in Engineering
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   , so the value of  is directly related to the tip diameter  = (  ).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   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 coefficient  ,RMS is 1.391, the system mass is 40.910kg, and the total safety factors is 8.966.

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 stage  1 , the pinion tooth number of the first and second stages  1 ,  4 , the helix angle of the first stage  1 , and the two shaft angles Φ 1 , Φ 2 , as expressed in (25): Other parameters can be evaluated by where  0 is the total gear ratio,  2 is the gear ratio of the second stage, where  1 ,  1 ,  4 ,  1 , Φ 1 , Φ 2 with the subscripts min and max are the boundaries of design variables, which are determined empirically according to the initial design.

Performance Constraints. (i) Contact and bending fatigue strengths should be below the allowable values:
where  1,2 is the contact stress of the first and second stages, [ 1,2 ] is the allowable contact stress of the first and second stages,  1,2 is the bending stress of the first and second stages, and [ 1,2 ] is the allowable bending stress of the first and second stages.
(ii) First stage gears should not interfere with each other: where ℎ *  is the addendum factor of the first stage gear.(iii) Safety margin of each stage should be balanced: where is the safety margin, with  min : the minimum safety factor for contact fatigue strength and  min : the minimum safety factor for bending fatigue strength; Δ is the mean value of safety margin Δ; [] is the allowable upper limit of safety margin standard deviation.

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.

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 individual  in the population has its fitness vector: and the fidelity () of the fitness vector, where  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 f tn () is evaluated using the original fitness functions, the fidelity () = 1; if f tn () is estimated, the fidelity 0 ≤ () < 1.
As shown in Figure 8, for each individual , specify its fitness sharing radius   .The area in which the dimensionless Euclidean distance between individual  and any other one is no greater than the fitness sharing radius is called fitness sharing area for individual , expressed as Ω().Assume there are  other individuals in the fitness sharing area Ω(), which composed a collection  = { 1 ,  2 , . . .,   }.And the evaluation method of f tn () is as follows.Firstly, evaluate the fidelity () of individual : where where  is weight rescaling factor.The closer the individual is to individual , the greater contribution of fidelity it makes.As depicted in Figure 9, if fidelity () is greater than a given threshold  * , then the fitness vector f tn () can be predicted as Else if () is less than the threshold  * , then the fitness vector f tn () should be evaluated using its original fitness functions.
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 where Δ  () is the coordinate difference (absolute value) of individual  between the former individual and the latter one in the th dimensional and  is the number of design Mathematical Problems in Engineering 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 individual  is with the predicted fitness vectors, let the fidelity of individual  in generation  be (, ), and then the fidelity in generation  + 1 can be updated as where  is fidelity drain factor, with 0 <  < 1.As the population evolves, the fidelity drops below a given threshold  0 , and the individual should also be removed from the database.

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 individual , and find the collection of individuals in the area from the database.
Step 3. Evaluate the fidelity () of individual , and determine whether () is greater than the threshold  * .If () ≥  * , predict fitness vector of individual  according to (35); otherwise, evaluate the fitness vector using its original functions and set the fidelity () to 1.
Step 4. Add individual  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.

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.
(3) POL Problem ( = 2) where and the variable bounds are [−, ] and the Pareto optimal front is nonconvex and disconnected.
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, threshold  * = 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 Figure 10, 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.
Percentage of real fitness vectors evaluated in each generation of the four benchmark problems are represented in Figure 11, 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.
The detailed testing results of effectiveness are illustrated in Table 2, 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.

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 Figure 12.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  1 is the L-S-C RMS  ,RMS , coordinate  2 is the total mass ∑  of the system, and coordinate  3 is the opposite number of total safety factors ∑ .
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 Figure 13.The -coordinate is the total mass of the system ∑ , the coordinate is the opposite number of total safety factors ∑ , and the information of L-S-C RMS  ,RMS is expressed by color.As Figure 13 illustrates, the Pareto optimal front is divided into 10 belt areas and L-S-C RMS  ,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  ,RMS varies from 1.166 to 1.709, the total mass of the system ∑  varies from 33.513 kg to 50.448 kg, and the total safety factors ∑  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.
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 Figure 14 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. 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.Figure 15 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   of solution 3# under possible operating conditions varies from 1.032 to The optimal design The reference design 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.

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.

Figure 1 :
Figure 1: Model of gear contact.

Figure 2 :
Figure 2: Full arrangement of a split torque transmission system.

Figure 3 :
Figure 3: Coordinates of a split torque transmission system.

Figure 4 :
Figure 4: An example of flank deviation surface for a tooth.

Figure 5 :
Figure 5: Mesh force curves of the first stage in both paths at given power and speed.

Figure 6 :
Figure 6: Mesh stiffness curves of the first stage in both paths at given power and speed. 400

Figure 7 :
Figure 7: "Load sharing map" of the reference design.

Figure 10 :
Figure 10: Pareto optimal fronts of the four benchmark problems with iNSGA-II.

Figure 11 :
Figure 11: Percentage of real fitness evaluated.

Figure 15 :
Figure 15: Load sharing map of 3# and the reference design.

Table 1 :
The geometry parameters of the reference design.
2,  6 are the gear tooth numbers of the first and second stages,  1 ,  2 are the module of the first and second stages, and  is the center distance between input and output shafts.The function round(⋅) here means round ⋅ to the nearest integer.

Table 2 :
Effectiveness test results of the algorithm.

Table 3 :
The values of design parameters and objective functions of solutions 1#∼4#.