Design and Optimization of PSD Housing Using a MIGA-NLPQL Hybrid Strategy Based on a Surrogate Model

A new power-split device (PSD) is designed as the core component of the multi-power coupling system in a hybrid electric vehicle. It is important to consider the influence of the factors of weight and stiffness on performance of PSD when designing PSD housing. In this paper, the overall arrangement of the power distribution system and PSD housing are redesigned. The finite element analysis is conducted to test PSD housing stiffness. According to the results of the analysis, this study adopts a hybrid optimization strategy based on surrogate models to obtain the least weight under the constraint of PSD housing stiffness. The surrogate models are established using a responsive surface method based on the data obtained by the optimal Latin hypercube design (OLHD). The hybrid optimization strategy combines a multi-island genetic algorithm (MIGA) with a nonlinear programming quadratic line search (NLPQL), which ensures obtaining optimal design parameters for PSD housing. In comparison with optimization using a single MIGA, this hybrid optimization strategy is more efficient and feasible for optimizing housing.


INTRODUCTION
The hybrid electric vehicle (HEV), an eco-friendly and energy-saving vehicle, has been a very popular in vehicle design and manufacturing [1] to [4].The power-split device (PSD) is a core component of a multi-power coupling system in a series-parallel HEV [5], because it can achieve energy coupling and conversion among the engine, motor, and generator under different working modes.
Based on the principle of differential, Yu et al. [6] proposed a PSD by which the multiple power coupling and decoupling can be achieved.In order to avoid the difficulty in the arrangement and installation of transmission gears in the PSD housing, the design of a PSD is explored in this paper.
In a traditional drive vehicle, the transmission is very compact, and the space is insufficient.However, new drive source and drive components are added in the hybrid electric drive vehicle, so the minimum optimization of the volume and quality is more important.The design and optimization of this kind of PSD are challenging.In particular, lightweight design has attracted intensive effort due to its great contribution to cost, material volume and time savings in engineering design [7] and [8].
With the field of complex engineering, numerous single optimization algorithms, such as the Karush-Kuhn-Tucker method [9], non-gradient optimization [10], genetic algorithm [11], and particle swarm optimization [12] are proposed to optimize structures.However, a single optimization algorithm does not always achieve a global solution [13] and can be time-consuming.Hybrid optimization algorithms, such as the combination of the global search (genetic algorithm or evolution strategies) and the local search (descent method) [14], the SASP method (hybridization of simulated annealing and the descent method) [15], the simulated annealing and the local proximal bundle method [16], the combination of the particle swarm optimization and the genetic algorithm [17], are of significant interest for the rapid speed of convergence and robustness in seeking globally optimal solutions [18].
In order to obtain a globally optimal solution and reduce computing time, this study adopts a hybrid optimization strategy based on the combination of a multi-island genetic algorithm (MIGA; a global optimizer), and a nonlinear programming quadratic line search (NLPQL; a local optimizer).In the MIGA, each group is divided into several subgroups called "islands".The selection, crossover, and mutation operations are performed in the subgroups, and the immigration operation is periodically performed among different targeted islands.The MIGA is a pseudo-parallel genetic algorithm, which can both avoid a local optimal solution as far as possible and accelerate convergence [19].When solving constrained nonlinear mathematical problems, the NLPQL algorithm shows stability, rapid convergence, and the capacity to seek globally optimal solutions [20].As the NLPQL can quickly determine the local optimal solution near the starting point, it can reduce the computing time during the optimization process.
This study proposes a new PSD and focuses on the sufficient level of stiffness of PSD housing, which ensures the minimum meshing misalignment between planet gear and half-axle gear.In such a case, a hybrid optimization strategy combining the MIGA and NLPQL based on a surrogate model is adopted.In this optimization, finite element analysis (FEA) is used to investigate the PSD housing stiffness.The surrogate models of weight and stiffness of the PSD housing are established.The hybrid strategy demonstrates higher efficiency and a stronger ability to find better solutions compared with a single MIGA.

Transmission System with a PSD
An integrated power distribution system including a PSD is shown in Fig. 1.The main parts of this system are the engine, the engine clutch, the gear reduction unit, the input clutch, the PSD, the motor, the gear acceleration unit for the generator, the generator clutch and the generator.The torque generated by the engine is introduced by the input axle and is passed to the PSD through the gear reduction unit.Then, according to the principle of the differential mechanism, the torque drives three planet gears and two half-axle gears.The input clutch connects with the left-housing, so that the parts of PSD housing cannot rotate.In this way, the bevel planet gear in PSD only rotates on own axis.The left side of the motor connects with the half-axis gear, and the right side connects with the rear axle.The electricity, which is produced by the generator, is stored in the battery.When the motor is working, the battery supplies electricity to it.The generator is connected to the generator clutch, and its operation can be stopped when the generator clutch disengages.

Structure of PSD
The PSD and the gear transmission train are shown in Fig. 2. The thrust bearing is applied at the shafts of each bevel planet gear to decrease the sliding friction between the bevel planet gears and the PSD housing; furthermore, the uneven distribution of the loading on the gear tooth is also improved.The left and right-housing is connected to the mid-housing of the PSD, which is composed of an outer ring and an inner ring.The ladder hole, which is opened on the inner ring, accommodates the roller bearing as well as makes the left spline axle through it.Hence, the driving torque transmitted to left and right half-axle gears becomes an output from one side of PSD.

Forces and Boundary Conditions of PSD Housing
Bevel planet gears and half-axle gears are the core transmission components of the PSD.The PSD housing deflection, which is introduced by the insufficient stiffness, usually leads to the meshing misalignment of the gears.The meshing misalignment increases the transmission error and transmission noise, and decreases transmission efficiency and gear life.In this paper, FEA is adopted to evaluate the PSD housing stiffness.
Before the finite element models are created, the loads applied to the PSD housing are computed.F al and F ar are the axial forces of the half-axle gears, which are applied to the left-and right-housing by half-axle gears.F tx (F t1 , F t2 , F t3 ) is the resultant force of the circumferential forces of the bevel planet gears, which is applied to the PSD housing.F ax (F a1 , F a2 , F a3 ) is the resultant force of the axial force created by the bevel planet gear meshing transmission and the centrifugal force created by the bevel planet gear revolution, which is also applied to the PSD housing.F al , F ar , F tx and F ax are defined as: The boundary conditions should be determined before a simulation is conducted.The engine transmits the input torque M (100 Nm) to the lefthousing, and the torque is distributed to the two halfaxle gears through the differential planetary gear train.In static load analysis, a hypothesis is proposed that the half-axle gears are fixed, so the torque M is applied to the left-housingof the PSD.F ax , F al and F ar can be converted into the pressure P acting on the corresponding plane, the loading results and the boundary conditions of the PSD housing are shown in Fig. 4.   Deflections of the inner and outer rings are shown in Fig. 7.The inner ring turns a certain angle relative to the outer ring, which leads the planet gears to deflect β angle.The outer ring expands b along the radial direction, which leads the planet gear to deflect along the axial direction.All of these events lead to the meshing misalignment among the bevel planet gears and the half-axle gears.The insufficient stiffness of the PSD housing affects the accuracy of gear meshing, which is considered to be a key constraint condition in designing the structural parameters of the PSD.Thus, this study aims to secure the PSD housing stiffness for decreasing the meshing misalignment when optimizing the weight of the PSD housing.

OPTIMIZATION BASED ON SURROGATE MODEL
The surrogate model method uses a simple mathematic model to replace a complex structural optimization problem, which can reduce computing costs and simplify problems.This paper adopts it in the design optimization of PSD housing.First, the design variables, the objective functions and the constraint conditions are determined; then, the design of the experiment based on the OLHD is implemented, and the surrogate models are created with response surface method according to the experimental data.

Design Variables
Based on the FEA results, the structure sizes of the inner ring, the outer ring and the radial plate have significant effects on the stiffness and weight of the PSD housing.Therefore, their sizes are selected as the design variables, as shown in Fig. 8.

Objective Function
To meet the lightweight design requirement, minimizing the weight of the PSD housing is considered the objective function, as shown in Eq. ( 5): where W is the weight of the PSD housing and X is the design variable matrix related to W.

Constraint Conditions
(1) Constraints of Design Variables: The ranges of design variables are determined under the geometric non-interference of the PSD housing structure, as shown in Table 1.

Table 1. Ranges of design variables
Design variables (2) Stiffness Constraint of the PSD Housing: According to the results of the FEA, the insufficient stiffness of the PSD housing will make the bevel gears deflect along the axle, which will lead to the stress concentration.The tooth surface is shaped like a drum in order to alleviate the above problems.In this study, the stiffness index is determined by comparing the changes of the contact stress along with the skew angles, which is between standard and modified tooth profiles.The value of the drum-gear profile [21] is presented in Eq. ( 6): where Δg [μm] is the drum size of the tooth profile, b [mm] is the gear width, f g = A(0.1b+10),A is the constant related to the accuracy magnitude of the gear as shown in Table 2.This paper adopts the accuracy magnitude 7, accordingly, A = 2.5.The gear parameters are used in Eq. ( 6) to obtain Δg = 14.2 μm.An axial modification is adopted [22], i.e. the theoretical involute of the two end surfaces of the modification gear rotates from the central axis to the gear tooth.The rotation angle is determined by the size of the drum-shape.Finally, the endpoints of the two theoretical involutes are connected into an arc.The arc radius is equal, which is used as the scanning trajectory of the tooth surface.With the variation of the skew angle β, the maximum contract stress of the modified and unmodified tooth profile is shown in Fig. 9.As shown in Fig. 9, when the axial skew angle is smaller, the contact stress of the modified gear is heavier.However, with the increasing of the skew angle, the maximum contact stress of the modified gear increases dramatically.When the skew angle is more than 0.04°, the maximum contact stress of the unmodified gear exceeds that of the modified gear.When the skew angle is 0.1°, the load distribution along the width of the gear tooth is very uneven, and the gear meshing line is shorter; the maximum contact stress increases by 85% compared with the normal meshing.
The skew angle β can be effectively reduced by improving the PSD housing stiffness.The uneven load distribution of the gear tooth can also be improved.Therefore, in the optimization process, the PSD housing stiffness is selected as the constraint condition to ensure that the uneven load distribution of the gear tooth, which is mainly caused by the insufficient housing stiffness, is minimized.
According to the above analysis, in the PSD housing optimization, the skew angle of the bevel gear is restricted to less than about 0.04°.The relationship between the skew angle of the bevel gear axis and the deflection of axis endpoint is presented in Table 3.By fitting calculation based on Table 3, the maximum deflection of the PSD housing is no more than about 0.02 mm.Thus, the constraint can be defined in the following form: 3.2 Establishment of the Surrogate Model

Design of experiment based on OLHD
The surrogate models are created based on the design matrix.In this study, the design matrix is generated using the design of experiment (DOE) and the real responses of the sample points are calculated by the FEA.The typical DOE includes the orthogonal design, uniform design, full factorial design, OLHD, etc.The OLHD improves the uniformity of the LHD by spacefilling and balance [23].Thus, OLHD is selected as the DOE in this study.

Surrogate Model
The surrogate model uses a simple mathematical model to replace a complex relation between the design variables and corresponding response during structural optimization process.At the same time, the calculation cost can be reduced effectively, and the results are very close to the true values [24] and [25].RSM uses different order polynomials to express the relationship between the design variables and their responses.One of the most used RSM is the second order polynomial model, which has the advantages of low computational cost, better approximation effect, being easy to solve, etc.Therefore, in this paper, the second order polynomial RSM is adopted to establish the surrogate models.The mathematical model of the second-order polynomial RSM can be expressed in Eqs. ( 8) and ( 9): where y is the output variable, x i is the design variable, n is the number of design variables, and A is the undetermined coefficient vector, α i , α ii and α ij are the regression coefficient.
The required sample points of RSM for establishing surrogate models are (n+1)×(n+2)/2, in which n is the number of design variables.In this optimization process, it has 8 design variables, so 50 sets of the sample points obtained using OLHD are used to establish the surrogate models of the PSD housing weight and maximum deflection.In the meantime, the other 20 sets of sample points, which can be used to test the accuracy of the surrogate models, are also obtained in the same manner.Based on the 50 sets of sample points, the surrogate models are established using the second-order polynomial RSM, and the coefficients of weight and deflection are presented in Table 4.
According to the sample points obtained using OLHD, the influence degree of the design variables on output W and U are obtained based on the statistics principle, as shown in Fig. 10a and b, respectively.As the surrogate model is an approximate equation between design variables and response function, conducting an error analysis of a surrogate model is necessary before replacing the true model.The accuracy of the surrogate model is usually measured via the multiple correlation coefficient R² in engineering application, which is defined as Eq. ( 10).The value of R² is between 0 and 1.The accuracy of the surrogate model is higher when R² is close to 1, and the accuracy is acceptable when the value of R² exceeds 0.9 in engineering applications [26].
where S R is the regression sum of squares, S y y . S E is the residual sum of squares,

S y y
. S T is the total deviation sum of squares, S T = S E + S R .
Using the error analysis, the R² values of the surrogate models of the housing weight and the maximum deflection equal 0.99995 and 0.96758, respectively, which indicates that the surrogate models can be used to predict performance.

Optimization Algorithm and Hybrid Optimization Strategy
A hybrid optimization strategy that combines the MIGA for global optimization and the NLPQL for gradient optimization is used to obtain the optimal results.Brief explanations of the MIGA and NLPQL algorithms are given, and the hybrid strategy is then explained.

Optimization Algorithm
A genetic algorithm (GA) is an effective method of global search optimization, which has been widely used in engineering optimization problems [27].The MIGA is an enriched algorithm based on GA.However, compared with a traditional GA, the MIGA is more efficient in finding the global optimum and more suitable for problems that have difficulty in obtaining gradient information.The MIGA is an exploration optimization method that uses "selection", "crossover" and "mutation" mechanisms to obtain the optimal design.The MIGA avoids obtaining a local optimal solution, thus can avoid the precocious phenomenon [28].
The NLPQL is used for gradient optimization to reduce the computing cost in the process of optimization.When the NLPQL is used to solve nonlinear mathematical problems, it has superior stability and fast convergence, and can easily obtain the global optimal solution [29].

Hybrid Optimization Strategy
The MIGA-NLPQL hybrid optimization strategy is performed as follows: firstly, the MIGA is used to determine the target area of the extreme value in the design space.Next, the NLPQL is used for the accurate optimization in target area defined by the MIGA to obtain the optimal design results.The flowchart of the hybrid optimization based on surrogate models is shown in Fig. 11.First, the advantages of the MIGA are brought into full play in rapidly roughly locating the sensitive target area in the overall design space, so that the low efficiency of the global optimization algorithm in detail can be avoided.Then, the global optimal results, which are obtained from MIGA optimization, are used in the NLPQL optimization module.Thus, the advantages of the gradient optimization algorithm of the NLPQL are brought into full play on local optimization, and the optimal solution can be found accurately.Misguided results, which are caused by seeking optimal results in highly nonlinear or discrete design spaces using the gradient optimization algorithm, can be avoided.

Optimization Results
In the optimization process, eight structural sizes are selected as design variables that belong to a multidimensional design, which is difficult to search in order to find the optimal value.To solve this problem, the MIGA-NLPQL hybrid optimization strategy is adopted in this study.In comparison with the optimization results of the MIGA-NLPQL and the single MIGA, the advantages of the hybrid optimization strategy are demonstrated.The iteration processes of the MIGA-NLPQL and MIGA are shown in Figs. 12 and 13, respectively.
The optimization process is time-consuming, so the number of optimization steps can determine its efficiency.Fig. 12 shows that the total iteration of the MIGA is 2000 steps.The objective function begins convergence at about 400 steps, and the global optimal solution is obtained at 1662 steps.Figs.13a and b show the iteration process of hybrid optimization; the iteration steps total about 1360.In the hybrid optimization, the MIGA calculates about 1000 steps and the NLPQL calculates about 360 steps.In Fig. 13a, a preliminary optimal solution is found at about 770 steps by the MIGA.Although the convergence of the objective function is not ideal, it shows a convergent trend.In this hybrid optimization strategy, the MIGA is only used to find the preliminary global optimal solution, so the convergence does not affect the NLPQL for further optimization.Fig. 13b shows the optimization process of the NLPQL.The objective function begins to converge when iterating at about 200 steps and the global optimal solution is finally obtained at 352 steps around the preliminary optimal solution found by the MIGA.Therefore, the total iteration steps of the hybrid optimization strategy are 1122, which is less than the 1662 steps of a single MIGA.In comparison with the iteration process of the MIGA, the computing cost is obviously reduced; therefore, the MIGA-NLPQL is more efficient than the single MIGA algorithm.The results of the optimization are shown in Table 5.The second row, "Initial", stands for the initial design without optimization; the third row, "MIGA", stand for the results based on the MIGA technique; the fourth row, "MIGA-NLPQL", stands for the results received from the hybrid strategy, "MIGA-NLPQL"."Rounded results" stands for the results after rounding according to the hybrid strategy of the MIGA-NLPQL.
Comparing the optimization results of the single MIGA and MIGA-NLPQL with the original design, it is shown that the weight of PSD housing is reduced by 26% after optimization using the single MIGA and by 26.6% after optimization using the MIGA-NLPQL.Furthermore, the computing time of the MIGA-NLPQL 1360 steps is less than the computing time of the MIGA 2000 steps.Therefore, the MIGA-NLPQL strategy can perform more efficiently to reduce calculation time and has a stronger ability to obtain optimal results than a single MIGA.The rounded results of the MIGA-NLPQL are selected as the final design parameters.Table 6 shows the geometric dimensioning and the performances of the original and the final designs of PSD housing.
The values of the objective functions U and W obtained by the surrogate models are 0.0095 and 1.1686, respectively.The optimal design is verified by FEA as shown in Fig. 14.The error between the surrogate model and the FEA of U is 13.6%, and that of W is 10.4%; which the error values of U and W are within the acceptable range, it can be verified that the surrogate model can be used to replace the performance of PSD.As shown in Fig. 14, the maximum displacement of FEA results is 0.011 mm, which satisfies the constrict condition of deflection.In  comparison with previous optimization, the maximum displacement is reduced by 56%.The maximum stress of the optimal PSD is 93.2 MPa, which occurs at the holes of end covers, and the stress of the joint between the radial plate and inner ring of the midhousing is improved.The above analysis shows that the optimization of PSD is reasonable, and the performances of PSD are improved.

CONCLUSIONS
This paper proposes a new PSD design scheme used in HEVs to improve the support form of bevel planet gears.To decrease the meshing misalignment of the internal gears caused by the lack of PSD housing stiffness, improving the is imperative.However, increasing stiffness and reducing weight are often contradictory.Therefore, the optimization design of PSD housing is performed, where the minimum weight is selected as the optimization objective and the stiffness is selected as the constraint.To enhance the optimization efficiency and obtain global optimal results simultaneously, the hybrid optimization strategy of the MIGA-NLPQL, based on a surrogate model, is adopted in this paper.In this way, the weight of PSD housing is reduced by 26.6% under the stiffness constraint.The design optimization shows that the MIGA-NLPQL hybrid optimization strategy has strong optimization ability and rapid convergence speed, as proven by the comparison between the single MIGA and MIGA-NLPQL results.

Fig. 1 .
Fig. 1.Integrated power distribution system of parallel-series HEV with PSD

Fig. 3 .
Fig. 3. Forces applied to PSD housing force of the bevel planet gear, m is the weight of the bevel planet gear, n is the rotational speed of the revolution, r is the reference radius of the half-axle gear, α is the pressure angle of the bevel planet gear, δ is the reference cone angle of the bevel planet gear, d m is the reference diameter of the bevel planet gear at the tooth width midpoint.

Fig. 4 .
Fig. 4. Loads and boundary conditions of PSD housing

2. 2
Results of FEAThe loads and boundary conditions are applied to the finite element models, and FEA is conducted.The results of Von Mises equivalent stress are shown in Fig.5.The maximum stress, approximately 92 MPa, occurs on the joint between the radial plate and inner ring of the mid-housing.Displacements are shown in Fig.6in the contour and vector forms.The maximum displacement is 0.025 mm.

Fig. 9 .
Fig. 9. Maximum contact stress of modified and unmodified tooth profile

Fig. 10 .
Influence degree of design variables on W and U; a) Influence degree of design variables on W, b) Influence degree of design variables on U

Fig. 11 .
Fig. 11.Flowchart of hybrid optimization based on surrogate models

Table 2 .
Constant A related to the accuracy

Table 3 .
Relationship between the skew angle of the bevel gear axis and the deflection of axis endpoint

Table 4 .
Coefficient of surrogate model of weight and maximum deflection

Table 5 .
Optimization results of MIGA-NLPQ and MIGA

Table 6 .
Comparison of the original and the final designs