Research on operating domain optimization of power split hybrid electric vehicle based on global bifurcation and chaos threshold

The power-split hybrid electric vehicle (PS-HEV) has multiple working modes to maintain high operation efficiency according to different conditions. The main modes involved in the vehicle driving process are pure electric mode and the hybrid driving mode. Because the electromechanical coupling problem is involved in the above two working modes, the transmission system exhibits strong non-linear characteristics. If the operation range of the engine and motor are unreasonable, the rotor system will vibrate and become instability. In this paper, the non-linear dynamic equations of the electromechanical coupling of the transmission system are established for electric driving mode and hybrid driving mode. The closed-homoclinic phase trajectory equation at the center point of the disturbance-free Hamilton system is determined. The chaotic thresholds for the pure electric and hybrid driving modes are derived through the Melnikov’s method to obtain the optimal working domain of the engine and motor. Finally, numerical simulation analysis is conducted to verify the feasibility of the work domain optimization scheme. Simulation results show that the proposed engine and motor working area optimization scheme can effectively avoid the homoclinic bifurcation in the PS-HEV during the driving process and prevent the vehicle from entering the chaotic state.


Introduction
The configuration of the hybrid vehicle powertrain differs significantly from that of the traditional automobile. For the hybrid vehicle, there is electromechanical coupling between the motor and the load, as well as the motor and the engine during the power transmission process, which demonstrates strong nonlinear dynamics. One of the key challenges of hybrid electric vehicles is the electromechanical coupling due to multiple energy sources. 1,2 If the operation range of the engine and motor are unreasonable, the rotor system will vibrate and become instability, even resulting in the shaft fracture of the power system. The nonlinear electromechanical coupling characteristics directly affect the performance, work efficiency, and safety of 1 School of Automotive and Traffic Engineering, Jiangsu University, Zhenjiang, China 2 Automotive Engineering Research Institute, Jiangsu University, Zhenjiang, China hybrid electric vehicles. Therefore, research on the dynamic characteristics and instability mechanism of hybrid systems has important theoretical significance on the system design and application.
With the development of the non-linear dynamics theory, non-linear characteristics of non-linear systems, such as vibration, instability, periodic solutions, bifurcations, and chaos, have become research hotspots. Hu et al. proposed an irregular internal excitation (engine excitation and motor excitation) and an external excitation (road excitation) would cause torsional vibrations of the transmission system. Besides, a simplified two-mass nonlinear dynamic model of a hybrid vehicle power system was established, and then applied to the prediction of the non-linear dynamics of the torsional instability range of the hybrid electric vehicle (HEV) powertrain. The analysis results were used to optimize the hybrid control strategy in order to improve the torsional stability of the hybrid vehicle power system under three typical driving modes. 3 The research on the torsional vibration of a hybrid vehicle power system mainly focused on the influence of the engine excitation, natural frequency analysis, and the corresponding modal analysis. Further, Hu et al. proposed a non-linear dynamic model of a hybrid power transmission system was established considering the multi-frequency excitation, and then the model was transformed into a fast-slow model with two-time scales. The De Moivre equation was introduced to unify the slow variables into a single variable, and the dynamic equations at different excitation frequencies and amplitudes were derived. Also, the bifurcation theory was used to study the bifurcation and stability analysis of the equilibrium point. Through the numerical analysis, the effects of the excitation frequency and amplitude on a dynamic behavior were studied using the equilibrium point curve, phase change phase diagram, and time history curve. The simulation results showed that the bifurcation would cause the jumping of the system trajectory and the corresponding blasting oscillation. With the continuous power and speed improvements, the electromechanical coupling effect became more and more prominent. 4 Zhou et al. studied the nonlinear dynamic characteristics of a system under the internal and external excitations were studied. The effects of the speed and load fluctuations on the system dynamic response were analyzed. The results showed that the increase in the excitation amplitude would cause the system to enter the chaotic motion through frequency doubling, and the speed change could effectively control system vibrations. 5 Liu et al. studied the dynamic characteristics of the gear transmission system under a multifrequency excitation based on the Hopf bifurcation theory and revealed the unstable mathematical mechanism. The largest Lyapunov index revealed the influence of the system parameters, including the friction coefficient and modal damping, on the dynamic response from the table period to the unstable, chaotic motion. 6 Zhao et al. studied in-depth research on the torsional vibration of wind turbine transmission systems was conducted. By using the maximum Lyapunov exponent and Poincare interface, the effects of the internal and external incentives on the system dynamic response were studied. The results showed that the external excitation has the greatest influence on the system dynamic characteristics, and changes in the excitation amplitude led to the system chaotic motion. In addition, it was also found that various internal stimuli, such as stiffness and damping, could significantly change the system dynamic response. 7-9 Bai et al. proposed a more comprehensive system model, including a nonlinear permeance network model (PNM), of a motor with a lateral torsion dynamic planetary gear system was established. The non-linear PNM considered the winding distribution of the motor, slotting of the stator and rotor, and non-linearity of the material. The lateral torsional dynamic planetary gear system considered the time-varying gear meshing stiffness, and the timevarying meshing and damping in the planetary gear system were studied in the proposed model. The simulation results were compared with the results predicted by the electromechanical model of a dynamic motor model with constant inductance. The results showed that the electromechanical coupling effect could cause additional severe vibration of the gear. In addition, it was found that external conditions, especially voltage transients, could have a significant impact on the dynamic characteristics of electromechanical systems. 10 Although the electromechanical coupling characteristics of the hybrid electric vehicle have been studied by many scholars, most of the studies focus on the influence of motor parameters, motor torque fluctuations, or engine torque fluctuations on the electromechanical coupling characteristics, or have only considered the electromechanical coupling process between the motor and the load. Namely, different working modes have been rarely included in studying the electromechanical coupling process between the motor and engine, and the motor and load. The chaotic phenomenon of the electromechanical coupling system was not considered during the mode switching coordination process. The main innovation of this paper is to analyze the system chaotic and bifurcation characteristics in pure electric mode and hybrid drive mode by establishing the electromechanical coupling torsional vibration model, so as to provide the basis for the torque distribution scheme during the mode switching of the PS-HEV. In the chaos analysis of the electromechanical coupling characteristics, the effects of different excitations, such as the engine, motor and the load, on the system are fully considered.
The remainder of this paper is organized as follows. Section I analyzes the configuration of the HEV powertrain and the corresponding working mode. Section II establishes the motor torque, engine torque models, and electromechanical coupling torsional vibration models of the power system in the pure electric and hybrid driving modes. Section IV analyzes the bifurcation and chaos thresholds of the two working modes, and Section IV determines the optimization scheme of the engine and motor operation area based on the chaos threshold. Section V conducts the numerical simulation analysis. Section VI concludes the paper.
The main contributions of this paper can be summarized as follows.
(1) The nonlinear dynamic equations of the power split system for the pure electric and hybrid driving modes are established to demonstrate the system electromechanical coupling characteristics. The equilibrium point of the disturbance-free Hamilton system is determined and analyzed. The closeness of the disturbancefree Hamilton system at the center Homoclinic orbit phase trajectory equation is obtained. (2) By adding the disturbance term to the Hamilton system, the chaos thresholds for the pure electric and hybrid driving modes are determined with Melnikov method to obtain the optimization scheme of the working domains of the engine and motor.

PS-HEV powertrain structure
The schematic diagram of the PS-HEV powertrain system studied in this paper is shown in Figure 1. The schematic diagram of the PS-HEV powertrain 11 studied in this paper is shown in Figure 2, where it can be seen that the transmission system consists of the engine, motor MG1, motor MG2, brakes CB1, CB2, and two planetary rows PG1 and PG2. The engine is connected to carrier C1, MG1 is connected to the sun gear S1, and MG2 is connected to the sun gear S2. In addition, the ring gear R1, carrier C2, and the output shaft are connected, and the ring gear R2 is fixed. Since the torque of MG1 can regulate the speed of the engine, the transmission in the power system is eliminated, leading to more compact configuration. The brakes CB1 and CB2 can control the engine and motor MG1, making them rotate freely or be fixed, and can also reduce the torque ripple of the engine and motor MG1. By controlling the working status of each power source and brake, the hybrid vehicle can switch between different working modes. According to the working status of the power source and the brake, eight working modes of the hybrid vehicle can be realized: electric driving mode, compound pure electric driving mode with MG1 and MG2, engine starting mode, hybrid driving mode, charging while standstill mode, compound braking mode, mechanical braking mode and stopping mode. Among them, the switching process from electric driving mode to hybrid driving mode involves various intermediate states such as engine startup and speed regulation. Compared with other driving mode switching, it is more complicated, as shown in Table 1.
The electromechanical coupling of the transmission system in the pure electric mode includes mainly the coupling between motor MG2 and load. On the other hand, the electromechanical coupling of the transmission system in the hybrid driving mode includes mainly the coupling between motor MG2, engine, and load. Therefore, the nonlinear dynamics of an electromechanical coupling system that are studied in this paper are as follows: (1) the nonlinear torsional vibrations caused by the electromechanical coupling between motor MG2 and load; (2) the nonlinear torsional vibrations caused by the electromechanical coupling between motor MG2, engine, and load.
Electromechanical coupling torsional vibration model of power system

Electromagnetic torque model of permanent magnet synchronous motor
In order to simplify the analysis, the following assumptions are made 12,13 : (1) The winding current changes with the time according to the cosine law, and the windings are star-connected; (2) the rotor is cylindrical, and the air gap is uniform; (3) the rotating air gap magnetic field ignores the effects of higher harmonics; and (4) the motor loss is neglected.
When the coil has the alternating current that changes with the cosine law, and the influence of harmonics is not considered, according to the theory of electrical machinery, the magnetomotive force of the three-phase symmetrical stator winding is expressed as 12,14 : where F sm denotes the fundamental amplitude of the stator magnetomotive force, and F sm =1.35NI m k N1 /p; N represents the number of per phase windings, I m denotes the current, pre presents the pole pairs of a PMSG; k N1 is the coefficient of the fundamental wave magnetization, and v 1 is electrical angular frequency of the motor. The magnetomotive force of the rotor is given by: where F rm represents the fundamental amplitude of the magnetomotive rotor force, and F rm = B r h m sin a p p 2 À Á =p=m 0 ; B r denotes the material remanence of the permanent magnet, h m is the permanent magnet thickness, m 0 represents the air permeability coefficient, and m 0 = 4p 3 10 À7 ; and a p denotes the arc coefficient of the permanent magnet.
The air-gap magnetic field energy of a permanent magnet synchronous motor is given by: where R denotes an average radius of the air gap, and l denotes the effective length of a rotor.
By considering the torsional vibration angle of the PMSG, the rotor synthesis fundamental wave magnetomotive force can be expressed as: If the rotor torsional angle is denoted as u, then the electromagnetic torque is expressed as: where F m = ppRlLF sm F rm , R denotes the inner radius of a stator, l represents the effective length of a rotor, and L is the air gap permeability of a PMSG.
The simplified equation of the electromagnetic torque for the PMSG can be obtained by Taylor approximation.
where k 0 = F m cosu, k 1 = pF m sinu, k 2 = p 2 F m cosu=2, k 3 = p 3 F m sinu=6, and k 0 , k 1 , k 2 , and k 3 represent the electromagnetic parameters. When the rotor torsional vibration angle is equal to zero, that is, u = 0, the motor is in a stable running state, and

Engine shaft torque model
The unbalanced torque the engine is mainly composed of the reciprocating moment of inertia T j and the combustion gas moment T g . 15,16 The variation law of the combustion gas torque is complicated, which can be expressed by Fourier series theory as follows: where T 0 denotes an average driving torque, and v 2 denotes the angular velocity of the shafting.
Since the fourth-and higher-order unbalanced moments account for only 2.5% of that of the secondorder, only the second-order unbalanced torque is considered and given by: Where m j represents the reciprocating equivalent mass of the connecting rod and piston group,r denotes the crank radius, D is the piston diameter, and l represents the connecting rod ratio is expressed as l = r=l.
The cosine component of the combustion gas moment T g is small compared with its sine component. When the torsional vibration angle is considered, the simplified torque model of the engine is given by: where T E0 denotes a constant part of the engine torque.

Torsional vibration model
Gear pair dynamics model. In order to study the electromechanical coupling characteristics of the PS-HEV powertrain, the dynamical model of the planetary gear is introduced and shown in Figure 3.
The general equation of the pure torsional dynamics of the planetary gears dynamics model that is presented in Figure 3 is as follows. 17 For the sun wheel: where u s denotes the rotational angle of the sun gear, J s represents the rotational inertia of the system equivalent to the sun gear, R s denotes the radius of the base circle for the sun gear, T 1 denotes the driving torque; C spi denotes the meshing damping between the sun gear and the ith planetary gear, k spi denotes the meshing stiffness between the sun gear and the ith planetary gear, i is the number of planetary gears, and d spi is the relative displacement between the sun gear to the ith planet gear. For planetary wheels: where u pi denotes the rotation angle of the planet gear, J p denotes the rotational inertia of the system equivalent to the planetary wheel, R P denotes the radius of the base circle for the planetary wheel radius, C rpi denotes the meshing damping between the ring gear and the ith planetary gear, k rpi denotes the meshing stiffness between the ring gear and the ith planetary gear, and d rpi is the relative displacement of the ring gear to the ith planetary gear.
For the planet carrier: where u C denotes the rotation angle of the planet carrier, J 0 C denotes the rotational of inertia of the system equivalent to the planet carrier, and J ; m pi is the mass of the ith planet wheel, and R C denotes the radius of the base circle of the planet carrier.
Also, d spi and d rpi can be respectively expressed as: where R p denotes the radius of the base circle of the planetary gear, u p denotes the rotation angle of the planetary gear, and R r denotes the radius of the base circle of the ring gear.
Nonlinear model of electromechanical coupling rotor in electric driving mode. According to the torsional vibration model of the gear pair given by equations (10)-(12) the electromechanical coupled torsional vibration model at the pure electric mode is established, 18 and it is given by: where T e0 denotes the constant part of the electromagnetic torque, T L0 denotes the constant part of the load torque, DT e denotes the electromagnetic torque disturbance, and DT L is the load torque disturbance.
According to the torsional vibration equation of the gear pair, equation (14) can be substituted into equations (10)-(12), we have: New variables are set to x = R s u s À u C ð Þ+ R p u p À u C À Á . Then, the nonlinear torsional vibration model of the hybrid powertrain is given by: Also, D represents the ring gear disturbance item, and Nonlinear model of electromechanical coupling rotor in hybrid driving mode. According to the torsional vibration model of the gear pair given by equations (10)-(12) the electromechanical coupled torsional vibration model at the hybrid driving mode is established, 3 and it is expressed as: where T E0 denotes the constant part of the engine shaft torque, and DT E denotes the engine torque disturbance.
According to the torsional vibration equation of the gear pair, equation (17) can be substituted into equations (10)- (12) which results in the following relations: New variables are set to x = R s u s À u C ð Þ+ R p u p À u C À Á . Then, the nonlinear torsional vibration model of the hybrid powertrain is given by: where m = Work area optimization based on global bifurcation and chaos threshold

Equilibrium point analysis of undisturbed Hamilton system
The excitations of the engine and load input, and the engine damping terms are regarded as disturbance.
Equations (16) and (19) are respectively rewritten as the form of the state equations 19,20 : where e denotes the small parameter, and 0\e\1; also, the value of eD is small enough and can be ignored here.
The potential function of the undisturbed system has an important on studying different nonlinear dynamic behaviors. Therefore, it is necessary to study the function of the undisturbed system in conducting electromechanical coupled nonlinear dynamics. When e = 0, the damping and disturbance terms are excluded, and equations (20) and (21) can be transformed into disturbance-free system, which can be expressed as follows: It should be noted that by discarding the damping and harmonic excitation terms, the transmission system can be converted into a Hamiltonian system. The two working modes have the same Hamilton form, and the properties of their equilibrium points are also the same. The Hamiltonian function is given by: And the asymmetric potential energy is expressed as: Let y = _ y = 0 in equation (24), we obtain the equilibrium points P 0 0, 0 ð Þ, P 1 x 1 , 0 ð Þ, P 2 x 2 , 0 ð Þ. Where: Then, the equilibrium point and its stability are analyzed according to the conditions of the system parameters g, b, and 1 À k ð Þ. Based on the theory of electrical engineering, the properties of a fixed point u s , 0 ð Þcan be determined as follows.
(1) If V 00 u s ð Þ\0, where V u s ð Þ is the maximum, the fixed point is the saddle point. There exit the trajectory to intersect the homoclinic orbit.
There exit heteroclinic orbit in the response of an unperturbed system.
There exit the trajectory to intersect the homoclinic orbit.
There exit heteroclinic orbit in the response of an unperturbed system.
Based on the parameters given in Table 2, b = 0:03, g = 0:25 and k = 1:5 can be calculated. Then, the equilibrium point of the disturbance-free Hamiltonian system satisfies the condition C1, that is, P 0 0, 0 ð Þ is the saddle point, and P 1 x 1 , 0 ð Þ and P 2 x 2 , 0 ð Þ are the center points. The potential function and phase trajectory of the undisturbed Hamilton system are shown in Figure 4. It can be seen that the undisturbed Hamilton system has a closed orbit at the center point and a homoclinic orbit. When the Hamilton system is disturbed the external interference, the homoclinic orbit is split by disturbance; thus, the system produces homoclinic bifurcation, which causes chaos.
In this case, P 0 0, 0 ð Þ is the saddle point, P 1 x 1 , 0 ð Þand P 2 x 2 , 0 ð Þ are the center points. Then, the Hamiltonian of the system is given by: When H 0, 0 ð Þ= 0, there are two homoclinic orbits q 0 t ð Þ È É connecting the saddle point 0, 0 ð Þ; thus, the homoclinic orbits satisfy the following conditions: Thus, the equations of x and y can be obtained, that is, the phase trajectory equations of the homoclinic orbit of the Hamilton system, and they are as follows: Chaos threshold determination According to the above analysis, the undisturbed Hamilton system transformed by the PS-HEV powertrain has the homoclinic orbits, and the PS-HEV transmission system has uncertain disturbances. When the disturbance terms is added to the Hamilton system, the closed homoclinic orbit is destroyed, which will lead to the same-direction unfolding bifurcation The existence of global homoclinic bifurcation often indicates the chaotic motion. Therefore, Melnikov method can be used to estimate the parameter space that the homoclinic bifurcation may be transformed into the chaos. In this way, the threshold of chaotic motion of the system can be determined. The essence of the Melnikov method is to calculate the distance between the stable and unstable manifolds after the disturbance of the homoclinic orbit of the undisturbed Hamilton system. Determine whether the stable and unstable orbits meet the conditions of crosssection intersection. Based on the given threshold conditions for the chaotic vibration of electromechanical coupling to determine whether there is chaos in the system. 21,22 Since the disturbance terms for the pure electric and hybrid driving modes are different, the Melnikov  analysis needs to be performed separately for these two working modes.
Melnikov analysis of electric driving mode. equation (20) can be transformed into: where X = x y ! , p 1 x, y ð Þ= y, and q 1 x, y, t ð Þ= 0, According to Melnikov's principle, the distance between stable and unstable manifolds can be expressed by the Melnikov equation, 23,24 as follows: where X h = x 6 , y 6 ð Þrepresents the homoclinic orbit of the Hamilton system, and p^q = p 1 q 2 À p 2 q 1 .
By substituting equation (28) into equation (29), the Melnikov integral is obtained as: where Þdt. The first integral I 1 can be directly calculated using the mathematical software, whereas the second integral I 2 is difficult to solve. Consequently, the numerical integration method is used in the paper to derive the integral I 2 . and the chaotic threshold curve are drawn by using the numerical integration method. Based on the threshold curve, the system can be qualitatively and quantitatively analyzed.
According to the Melnikov theory, ifM t 0 ð Þ= 0, The condition of producing chaos with a fixed frequency is expressed as: Assuming f is the control parameter and m = 0:08, the Melnikov threshold curve is shown in Figure 5. The three-point f -valuesin Figure 5 are f = 1:4, f = 1:6, and f = 1:8. v 1 = 610 corresponds to f = 1.58 on the threshold curve, and the corresponding amplitude of the load excitation is 120 N Á M. In Figure 5, it can be seen that the change in f determines whether the system appears chaotic motion. When f \1:58, the load excitation is smaller than 120 N Á M, the system is in the stable state with periodic vibration, and there is no chaos; on the other hand, when f ø 1:58, the load excitation is F 1 ø 120 N Á M, the system passes through the stable region and then enters into the chaotic area, and the chaos phenomenon occurs.
Melnikov analysis of hybrid driving mode. Assume v 1 = 2v 2 , and express equation (21) as: where X = x y ! , p 1 x, y ð Þ= y, and q 1 x, y, t ð Þ= 0, Þ . According to Melnikov's principle, the distance between stable and unstable manifolds can be expressed by the Melnikov equation as follows: where X h = x 6 , y 6 ð Þrepresents the homoclinic orbit of the Hamilton system, and p^q = p 1 q 2 À p 2 q 1 .
By substituting equation (29) into equation (30), the Melnikov integral is given by: Where Similarly, the first integral I 1 can be directly calculated using the mathematical software, whereas the second and third integral I 2 and I 3 are difficult to solve, so in this work, the numerical integration method is used to solve them. Although the corresponding analytical formula cannot be obtained, the chaotic threshold curve can be drawn using the numerical integration method. According to the threshold curve, the system can also be qualitatively and quantitatively analyzed.
According to the Melnikov theory, ifM t 0 ð Þ= 0, The condition of producing chaos with fixed frequency is expressed as: À1\ mI 1 fI 2 + f 2 I 3 \1 Assume f 1 is a control parameter and m = 0:08, then the Melnikov threshold curve is as shown in Figure 6.
The three-point f -values in Figure 6 are f 1 = 1:5, f 1 = 1:7, and f 1 = 1:8; where v 1 = 715 corresponds to f 1 = 1.62 on the threshold curve, and the corresponding load excitation amplitude is 125 N Á M. In Figure 6, it can be seen that the change in f 1 determines whether the system appears chaotic motion. When f 1 \1:62, the load excitation is F 1 \125 N Á M, the system is in a stable state of periodic vibration, and there is no chaos; on the other hand, when f 1 ø 1:62, the load excitation F 1 ø 125 N Á M, the system passes through the stable region, and then enters into the chaotic area; thus, the chaos phenomenon arises.
Assume f 2 is a control parameter and m = 0:08, then the Melnikov threshold curve is as shown in Figure 7. The blue line represents the relationship between the engine excitation amplitude and the angular velocity v 1 , and the red line represents the stable boundary. When v 1 \887 or v 1 .1079, f 2 is in the stable region at this time, and the system runs stably. When 887\v 1 \1079, the system runs in an unstable area, and chaos occurs.

Engine and motor working area optimization scheme
According to the Melnikov analysis of the PS-HEV powertrain in pure electric and hybrid driving modes, the optimization scheme of the engine and motor working domain based on the global bifurcation and chaos threshold is designed. It is assumed that the load excitation is 120 N Á M and that it always smaller than 120 N Á M during the vehicle running process.
In the pure electric mode, the speed of motor MG2 and the corresponding working torque are 5830r=min v = 610rad=s ð Þand 74 N Á M to avoid chaos. If the speed is too small, and the system will cross the chaos threshold curve, the system will enter the chaotic region. Therefore, in the pure electric mode, the torque of motor MG2 should be controlled to be not more than 74N Á M, that is T MG2 ł 74 N Á M.
In the hybrid driving mode, the speed of motor MG2 and the corresponding working torque needed to ensure that the chaos does not appear are 6800r=min v = 715rad=s ð Þ and T MG2 = 63 N Á M. Therefore, under the condition that the required torque of the vehicle is satisfied, the torque of the control motor MG2 is T MG2 ł 63 N Á M. If the torque T MG2 provided by motor MG2 is greater than 63 N Á M, the system will enter the chaotic region, and chaos will occur. At the same time, according to the established control strategy, the engine rotational speed is 1500r=min ł n E ł 6000r=min, and the angular velocity is 157rad=s ł v 2 ł 628rad=s. It has been analyzed that 887rad=s\v 1 \1079rad=s, that is, when 443:5rad=s ł v 2 ł 539:5rad=s, the system runs in an unstable area. Therefore, the optimal angular speed range of the engine to ensure that the system does not appear chaotic is 157rad=s ł v 2 ł 443:5rad=s, and the corresponding optimal speed range is 1500r=min ł n E ł 4237r=min.
In summary, when the vehicle is operating in the pure electric mode, MG2 provides the driving torque to drive the vehicle. The torque of MG2 is smaller than 74N Á M. Under this condition, if the torque provided by MG2 is insufficient, the mode should be changed and switched to the hybrid driving mode. When the vehicle is operating in the hybrid driving mode, MG2 and the engine provide the driving torque simultaneously, the torque of MG2 is T MG2 \63N Á M, and the remaining driving torque is provided by the engine. The optimal speed is n E = 1500r=min. Combined with the fuel economy operation area of the engine, the  speed range of the control engine is 1500r=min\n E \4237r=min.

Numerical simulation verification
In order to verify the accuracy of the results obtained by the above-mentioned Melnikov analysis, Runge-Kutta method is used to verify the torsional vibration equations (20) and (21) of the transmission system. The calculated time response history curve (excluding the transient instability part at the beginning of the simulation), phase trajectory diagram, Poincare section diagram and amplitude spectrum are analyzed. The purpose of the numerical simulation is to verify the effectiveness of the proposed work area optimization scheme.

Numerical analysis of electric driving mode
When T MG2 = 65N Á M, the time history, phase plan, Poincare map, and Fourier spectra obtained by the numerical analysis are shown in Figure 8. It can be seen that the time response history curve has an orderly cycle, the phase trajectory is a circular circle, the Poincare map has only one point, and only one energy pulse appears in the Fourier spectra. This indicates that the transmission system is in a periodic motion state, and the maximum dynamic tangential displacement of the system is relatively small. At this time, the system has not lost its stability, and no bifurcation and chaotic motion have occurred.
WhenT MG2 = 85N Á M, the time history, phase plan, Poincare map, and Fourier spectra obtained by the numerical analysis of the transmission system are shown in Figure 9. As shown in Figure 9, the time history curve starts to be chaotic and disordered. The phase trajectory has multiple intersecting loops. There are multiple points on the Poincare map. The Fourier spectra has regions of energy concentration. This indicates that the timedomain response of the transmission system is non-periodic, and the maximum dynamic tangential displacement of the system is relatively large, which indicates that the system loses stability, a global homoclinic bifurcation occurs, and the system is in chaos.
The numerical simulation analysis results show that in pure electric mode, the torque of MG2 was T MG2 ł 74N Á M, and it could ensure the system stability without bifurcation and chaotic motion. Simulation results are consistent with the results obtained by the Melnikov analysis method.

Numerical analysis of hybrid driving mode
When T MG2 = 55N Á M, v 2 = 400rad=s, the time history, phase plan, Poincare map, and Fourier spectra obtained by the numerical analysis of the transmission system are shown in Figure 10. It can be seen that the time response history curve has an orderly cycle, the phase trajectory is a circular circle, the Poincare map has only one point, and only one energy pulse appears in the Fourier spectra. This indicates that the transmission system is in a periodic motion state, and the maximum dynamic tangential displacement of the system is relatively small. At this time, the system has not lost its stability, and no bifurcation and chaotic motion have occurred.
When T MG2 = 70N Á M, v 2 = 500rad=s, the time history, phase plan, Poincare map, and Fourier spectra obtained by the numerical analysis of the transmission system are shown in Figure 11. As shown in Figure 11, the time history curve starts to be chaotic and disordered. The phase trajectory has multiple intersecting loops. There are multiple points on the Poincare map. The Fourier spectra has regions of energy concentration. This indicates that the time-domain response of the transmission system is non-periodic, and the maximum dynamic tangential displacement of the system is relatively large, which indicates that the system loses stability, a global homoclinic bifurcation occurs, and the system is in chaos.
The numerical simulation analysis results showed that in the hybrid driving mode, the torque of motor MG2 was T MG2 ł 63N Á M, and it was controlled. Because the engine was involved in the driving process, in order to ensure the stability of the transmission system, no bifurcation and chaotic motion occurred, and the engine could work in a low fuel consumption area, achieving good fuel economy, while controlling the engine output torque to ensure meeting the driving requirements, and the engine speed range was 1500r=min\n E \4237r=min. These results are consistent with those obtained by the Melnikov analysis method.

Conclusion
In this paper, a PS-HEV is studied, and its two working modes that are involved in the mode switching process are analyzed. Also, the bifurcation and chaos characteristics of the hybrid vehicle transmission system are studied in different working modes, and the electric driving mode and hybrid driving mode are analyzed by Melnikov's method. In addition, based on the chaos threshold, the working domain optimization scheme of the engine and motor is obtained. Finally, numerical simulation analysis is used to verify the feasibility of the work area optimization scheme. chaotic motion. The homoclinic can be estimated by Melnikov's method, and the bifurcation may be transformed into a parameter space of chaos so that the threshold of the chaotic motion of the torsional vibration of electromechanical coupling can be determined. In the pure electric mode, in order to ensure the stability of the transmission system without bifurcation and chaotic motion, the torque of motor MG2 should be T MG2 ł 74N Á M; and, in  the hybrid drive mode, it should be T MG2 ł 63N Á M. Because the engine is involved in the driving process, no bifurcation and chaotic motion occurred, and the engine could work in the low fuel consumption area, achieving better fuel economy. By controlling the engine output torque to ensure the driving requirements, the engine speed range was 1500r=min\n E \4237r=min.
However, there are still some shortcomings that need to be addressed. For instance, a simplified construction of the transmission system mode should be obtained, and the internal motor parameters, the number of pole pairs of the motor, the magnetic circuit saturation coefficient, and other electromagnetic characteristics should also be considered.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.