Optimization of Geostationary Orbit Transfers via Combined Chemical–Electric Propulsion

: For geostationary orbit transfers, a long duration is required using electric propulsion and a large propellant mass is needed with chemical propulsion. Hybrid transfers can achieve a balance between the fuel consumption and transfer time. In this paper, a trajectory optimization method is proposed for time-ﬁxed minimum-fuel orbital transfer with combined chemical–electric propulsion. The necessary conditions and transversality conditions related to impulsive burns are derived theoretically with Pontryagin’s maximum principle. The long-duration geostationary orbit transfer is a many-revolution transfer, and is solved with the homotopic approach from the short-duration transfer problem. The variation in fuel consumption with transfer time is nearly linear, and the variation in the magnitude of impulsive burn is exponential. A simple model is presented for the estimation of fuel consumption and magnitude of impulsive burn with given transfer time, speciﬁc impulse of propulsion system and low-thrust magnitude.


Introduction
Most commercial communications satellites, some weather satellites and navigation satellites are placed in geostationary orbit (GEO). The transfers to GEO are generally accomplished by a chemical propulsion (CP) system. The satellites are injected into geostationary transfer orbits (GTOs). Then, they are boosted into GEO using CP in several hours or several days. This transfer is fast but requires a large amount of fuel due to its low specific impulse. In the 1960s, electric propulsion (EP) was introduced as a practical alternative to CP, being able to reduce fuel consumption at the cost of a longer mission time [1]. There have been many works about the optimization of all-electric orbit transfers [2][3][4]. Most of the optimization methods can be generally categorized into direct methods, indirect methods and hybrid methods, which combine the two. In the direct method, the optimization problem is converted into the parameter optimization problem and then nonlinear programming (NLP) can be applied to obtain the solution, which requires a large amount of computation [2,5]. As for the indirect method, the optimization problem is usually solved through a two-point boundary-value problem (TPBVP) or multi-point boundary-value problem (MPBVP), which is difficult to solve due to the small convergence radius [6,7]. Hybrid methods are combinations, which have the advantages of both methods, where the time histories of costate variables are directly parameterized and the optimal solution is obtained through the NLP method [8,9]. Apart from these trajectory optimization methods, the heuristic control law has also been applied to achieve low-thrust orbit transfer. Although the obtained solutions are usually not optimal, the heuristic control law has the advantages of simple formulation and quick calculation speed [4,10,11]. Moreover, the shape-based method was developed to determine the analytical solutions to the low-thrust transfer trajectory [12]. Boeing's 702SP satellite was the first all-electric propulsion satellite launched into GEO [13]. Although EP is fuel saving, the transfer time is also essential because of the economic benefits. Therefore, all-electric satellites are not applicable in some GEO transfer missions when the mission time is limited. Combined chemical-electric multimode propulsion has been the subject of most recent interest [14]. This type of propulsion has shown benefits for commercial spacecraft, especially for orbit-raising missions. It has been found that the advanced onboard propulsion systems can be used to perform both the north-south station keeping and part of the orbit transfer to GEO [15]. This work shows that the use of EP for portions of the transfer allows a greater net mass of spacecraft by extending the transfer time, in comparison with chemical-only transfer. A hybrid orbital transfer from low Earth orbit (LEO) to GEO was investigated by Mailhe and Heister [16]. The hybrid transfer was proven to be effective in reducing the trip time and radiation damage in the Van Allen belts, and dramatically reducing the vehicle gross weight. Both high-low and low-high-low thrusting strategies are studied in their investigation, and the low-high-low thrusting strategy was proven to be more efficient. Oh et al. developed a simple analytic multistage model for combined chemical-electric orbit-raising missions [17]. They applied the analytic model to estimate the performance of the transfer, along with a graphical-analytic method for a more accurate estimation. Jenkin performed trade studies in a case study of GEO transfer using combined chemical-electric propulsion [18]. In his work, trajectory optimization was used to maximize the GEO insertion mass. His trade studies determined trends among various missions and system parameters, such as the elliptical orbit for the start of the EP phase, the input power for the EP system and the specific impulse of the system. Kluever considered the mission sequence beginning with injection into an elliptical orbit with an arbitrary apogee altitude, an apogee burn using the onboard chemical stage to raise perigee and/or reduce inclination and a low-thrust transfer to geostationary orbit [19]. Kluever also proposed a purely analytical algorithm that can determine spacecraft mass requirements for a desired electric propulsion system and desired transfer time [20].
The main difficulty of hybrid transfer is how to determine the time, magnitude and direction of impulsive burns to minimize the fuel consumption of the whole transfer. Previous works have discussed the relations between the EP transfer stage and starting orbit parameters, such as apogee, perigee altitudes and inclination. With these relations, the delta-v of the chemical stage for a desired EP stage can be determined with Hohmann transfer or some other methods. In this paper, the optimizations of EP stage and CP stage are combined together. Different from previous works, optimal control theory is employed here to determine the timing, magnitude and direction of the impulsive burns along with the low-thrust control law. The impulsive burns can lead to discontinuities in state variables at those moments. The necessary conditions and transversality conditions related to impulsive burns are derived theoretically with Pontryagin's maximum principle (PMP). The time-fixed minimum-fuel hybrid transfer can be turned into a MPBVP with PMP. This optimization method is similar to the indirect method in the optimization of low-thrust transfer. Moreover, as the long-duration transfer from GTO to GEO is many-revolution, the homotopic approach is taken to solve long-duration transfer with the solution of shortduration transfer. Moreover, a model to estimate the performance of the hybrid transfer is developed here for the determination of the parameters of the propulsion system to meet the mission's requirements, such as the mission time and fuel consumption. The variations in the fuel consumption and magnitude of impulsive burn are studied to build the interpolation model. With this model, the performance of the hybrid transfer, such as the fuel consumption, magnitude of impulsive burn and transportation rate, can be estimated quickly, while the errors between the estimated values through this model and the accurate solutions through the indirect method are relatively small. There has been some trajectory optimization software developed by the direct method, indirect method or some other methods [21]. For example, GMAT (General Mission Analysis Tool) is a space mission design software for computing multigravity-assisted interplanetary and Earth orbit transfers in accurate dynamical models that implement direct methods combined with gradient-based solvers [22]. In the tool SEPTOP (Solar Electric Propulsion Trajectory Optimization Program), indirect single shooting has been implemented [23]. However, it would be difficult to use the existing software directly for the optimization of the timing, magnitude and direction of the impulsive burns. This work can be used to provide the necessary information for this software to further design the transfer orbit.
In this investigation, the optimization of GEO transfers via combined chemical-electric propulsion is based on the following assumptions: (i) no effect of Earth shadow eclipses is considered, (ii) the chemical burns are modeled as impulsive burns, (iii) switch time between the work of chemical propulsion and electric propulsion is ignored, (iv) electric propulsion is turned on for the whole trip (continuous thrust).
This paper is organized as follows: in Section 2, the method for solving time-fixed minimum-fuel transfer using combined chemical-electric propulsion is introduced. The optimization result and estimation model of hybrid transfer from GTO to GEO with one impulsive burn is presented in Section 3.

Dynamic Model
In order to avoid singularities, the orbit of a spacecraft is described using modified equinoctial elements (MEE) [24]: where a is the semi-major axis, e is eccentricity, i is inclination, Ω is the longitude of the ascending node, ω is the argument of perigee, and ν is a true anomaly.
In this study, the direction of thrust is described in the RSW (where the unit vectors R and S in the orbital plane are along the local radial and transverse direction, and W is along the angular momentum direction) reference frame, where the x axis is oriented along the position vector, z along the angular momentum and y perpendicular to the radius vector in the direction of motion, and r RSW , v RSW are used to denote position and velocity vectors in the RSW reference frame. For the preliminary design, it is assumed that the spacecraft is only subject to the central force of Earth's gravity and the thrust of the electric propulsion system. The differential equations of motion are given as [24]: where s L = sin L, c L = cos L, µ is the gravitational constant of Earth, m is the mass of the spacecraft, T max is the maximum of the thrust, I sp is the specific impulse of EP, g 0 is the gravitational acceleration of Earth on the sea level, u ∈ [0, 1] is the thrust ratio, and the unit vector α denotes the direction of thrust, w = 1 + f c L + gs L and s 2 = 1 + h 2 + k 2 .
In this study, all the quantities are nondimensionalized by the Earth's radius R E , Earth's gravitational constant µ and the spacecraft's initial mass m 0 . Let x = [p, f , g, h, k, L], A and B, respectively, denote a 6 × 1 vector and 6 × 3 matrix on the right-hand side of Equation (2). The differential equations of motion can be rewritten in a compact form:ẋ Actually, matrix B is the derivative of MEE with respect to the velocity state in the RSW reference frame, i.e., B = ∂x ∂v RSW . For low-thrust transfer, there are boundary conditions as follows: where t 0 and t f are, respectively, the initial and final times. The terminal boundary conditions can also be denoted with MEE: In this problem, the inner constraints regarding impulsive burns should be considered. Assuming that the number of impulsive burns is n, the inner constraints are addressed as where s [r; v], r and v are position and velocity vectors in the Earth-centered inertial (ECI) reference frame; ∆s i [0 3×1 ; ∆v i ], ∆v i denotes the change in velocity vector in the ECI reference frame conducted by the impulsive burn; ∆v i = ||∆v i ||, t − i and t + i are, respectively, time just before and after the i-th impulsive burn; r and v represent the position and velocity vectors of the spacecraft in the ECI reference frame, and I spc is the specific impulse of CP.

Optimal Control Theory
The purpose of this problem is to minimize the fuel consumption. The mass of the spacecraft after the i-th impulsive burn is derived as (normalized by m 0 ) Define the n + 1-th impulsive burn at t f as ∆v n+1 , and assume ∆v n+1 = 0. The final mass of the spacecraft can be expressed as Therefore, the performance index is As there are constraints ψ i and φ i , Lagrange multipliers χ i and κ i are applied to transform it into an unconstrained problem. With the assumption that the electric propulsion system turns on for the whole trip, the performance index of the unconstrained problem is given as For convenience, define The Hamiltonian can be built as where λ λ p ; λ f ; λ g ; λ h ; λ k ; λ L is the costate vector related to MEE, λ m is the costate related to mass, and the superscript i in H (i) and L (i) means the expressions when t i < t < t i+1 (i = 0, 1, 2, . . . , n). According to PMP, the optimal thrust direction and magnitude, which minimize the Hamiltonian, are determined by where SF (i) is the switching function, which holds the form Substituting the optimal control law shown by Equation (17) into the Hamiltonian function, the costate differential equations, which are also called the Euler-Lagrange equations, can be obtained as [25] where the derivatives ∂A ∂x and ∂B ∂x are shown in Appendix A. Because the final mass m t f and L t f are not constrained, the transversality conditions show that the final costates λ m t f and λ L t f should satisfy Now, the necessary conditions related to impulsive burns can be derived from Equation (13). The first-order variation in performance index is [25] where . . . denotes the terms related to Euler-Lagrange equation (20) and control variables u and α, which have been introduced before. It is obvious that ∂Φ ∂t i = 0. Then, let us analyze the derivatives ∂Φ dt. Applying Equations (9) and (14), the following relations can be deduced: The expression of matrix ∂s ∂x is shown in Appendix A. In the remaining part, A sx is used to denote matrix ∂s ∂x . Substituting Equations (23)- (27) into Equation (22), λ(t), λ m (t) and ∆v i should satisfy The Lagrange multipliers χ i can be determined with λ t − i and A sx t − i according to Equation (28), and κ i can be determined using Equation (30). Costates λ t + i and λ m t + i are updated using Equations (29) and (31). From Equation (32), the direction and magnitude of impulsive burn ∆v i satisfy Equation (34) shows the direction of the impulsive burn in the ECI reference frame. Lagrange multiplier χ i is derived from Equation (28) as Combining Equations (17), (34) and (36), it can be deduced that where α ECI is the direction of low thrust in the ECI reference frame. Equation (37) proves that the direction of the impulsive burn stays the same as the direction of low thrust just before and after the impulsive burn in the ECI reference frame. The differential equations shown by Equations (4), (5) and (20) are numerically integrated using the classic eighth-order Runge-Kutta integrator with a seventh-order RKF7(8) for automatic step-size control [26]. The unknowns include seven initial values of costates λ(t 0 ), λ m (t 0 ), the timing and magnitude of impulsive burn t i , ∆v i (i = 1, 2, . . . , n). Numerical Lagrange multipliers χ i , κ i (i = 1, 2, . . . , n) are obtained using Equations (28) and (30). The direction of impulsive burn can be determined with Equation (34). At the same time, there are also equations consisting of the 5-D Equation (8), the 2-D Equation (21), transversality condition Equation (33) and first-order necessary condition Equation (35). States  (29) and (31). Hence, the solution of the MPBVP can be obtained using the numerical method for solving nonlinear equations.

Steps of Solving Hybrid Transfer Problem
Similar to the all-electric transfer problem, the hybrid transfer problem is difficult to solve when there are many revolutions. For short-duration transfer, the initial guess of the solution can be obtained by the global optimization method improved cooperative evolutionary algorithm (ICEA) [27]. However, as the solution of MPBVP is sensitive to the initial guess, it is difficult to solve the long-duration transfer problem with the initial guess obtained by global optimization. To deal with this issue, homotopy techniques, also called numerical continuation, as introduced in the all-electric transfer problem [28,29], are employed. In this work, the long-duration transfer problem is solved through transfer time continuation and thrust continuation.
The steps of solving the hybrid transfer problem can be concluded as the following.

1.
Apply the ICEA algorithm to search the approximate solutions of the initial costate vector [λ(t 0 ); λ m (t 0 )], the timing and magnitude of impulsive burn [t i ; ∆v i ] (i = 1, 2, . . . , n) and numerical Lagrange multipliers χ j to the high-thrust short-duration transfer prob-

2.
Take the approximate solutions of step 1 as the initial guess to provide the accurate solution of the high-thrust short-duration transfer problem Γ T

5.
When t (j) f = t f , output the desired solutions of the hybrid transfer problem Γ T max , t f . The time continuation process and thrust continuation process are accomplished separately in the algorithm above. The two steps can also be combined together by reducing the thrust and increasing the transfer time at the same time.

Results
Orbit transfer from GTO to GEO via combined chemical-electric propulsion is discussed in this section. The gravitational acceleration of Earth on the sea level is g 0 = 9.806 m/s 2 , the specific impulse of chemical propulsion is I spc = 300 s, and the initial mass of the spacecraft is m 0 = 800 kg. The initial orbital elements of the spacecraft are listed in Table 1.

Solution of Hybrid Minimum-Fuel Hybrid Transfer
Firstly, hybrid transfer orbits with thrust magnitude T max = 200 mN and low-thrust specific impulse I sp = 3000 s are optimized with various transfer times. MinPack-1 [30] and GNU Scientific Library (GSL) are used to solve the MPBVP. In this way, the convergence and efficiency of this algorithm can be discussed with different solvers. In GSL, the hybrid solver, discrete Newton solver and Broyden solver are considered to solve the MPBVP [31][32][33]. Moreover, the hybrid solver in GSL is based on MinPack. The computations are executed on a desktop personal computer with a CPU of 1.8GHz, and Visual Studio 2019 is used.
Considering the homotopic approach from t (j) f days, and the maximum iteration number of the solvers in GSL is 100. When it does not converge, set h t f = 0.8h t f . Part of the performance of the algorithm with the four solvers is listed in Table 2. According to Table 2, this homotopic approach converges with all of the four solvers. The fuel consumption of the solution with four solvers is the same. However, the efficiency of solving this algorithm with MinPack-1, the hybrid solver and the Broyden solver in GSL is slightly higher than with the discrete Newton solver in GSL. Therefore, it is suggested to use Minpack-1, the hybrid solver and the Broyden solver in GSL to solve this problem. In the following works, MinPack-1 is used to solve this MPBVP. In this case, chemical-only transfer orbit and all-electric transfer orbit are optimized first for comparison. The minimum time transfer problem with electric propulsion is solved by the indirect method [34] and thrust continuation. The fuel consumption of chemical-only transfer with two impulsive burns is ∆m c = 366.86 kg. Moreover, the magnitude of the two burns is, respectively, 1.805 km/s and 0.22 m/s. The minimum transfer time from GTO to GEO with all-electric propulsion is 115.942 days, and the fuel consumption is 68.10 kg.
For hybrid transfer with two impulsive burns with a transfer time of several days, the solution shows that the first impulsive burn is the same as that of hybrid transfer with one impulsive burn, and the magnitude of the second impulsive burn is 0. Moreover, the positions where the two impulsive burns are conducted are the same as those of chemical-only transfer. The first impulsive burn is conducted at the apogee to raise the perigee and lower the inclination, which is also the descending node. For GTO orbits, the apogee or the descending node is usually close to the GEO orbit (usually with an apogee of 42,164 km), which means that the delta-v of the second impulsive burn is almost 0 for chemical-only transfer. For hybrid transfer from GTO to GEO, the very small delta-v of the second impulsive burn can be completely replaced with the delta-v of EP in a short transfer time. As the magnitude of the second impulsive burn of hybrid transfer with two impulsive burns is 0, only hybrid transfer with one impulsive burn is considered for this transfer mission.
A series of hybrid transfer problems with transfer time varying from several days to 115 days are optimized. There are multiple local optimal solutions for the hybrid transfer problem. In the time continuation process, the timing of impulsive burn t 1 varies in a small range. When beginning from solutions with different t 1 , the homotopy approach will lead to different solutions. In this work, two series of solutions with different values of t 1 are studied. The variation in fuel consumption ∆m with transfer time t f and variation in the magnitude of impulsive burn ∆v 1 with t f are demonstrated in Figure 1. The timing of impulsive burn t 1 of the first series of solutions is 1.10 days, and t 1 of the second series of solutions is 0.66 days. According to Figure 1, when the transfer time is short enough, the hybrid transfer is close to chemical-only transfer, which requires fuel consumption of 366.86 kg. When the transfer time is close to the minimum transfer time with allelectric propulsion, the performance of hybrid transfer orbit is close to that of all-electric transfer orbit. The local optimal solutions with t 1 = 1.10 days perform better than those with t 1 = 0.66 days. Therefore, only solutions with t 1 = 1.10 days are discussed in the following section. As is shown in Figure 1, there is a linear relation between fuel consumption ∆m and transfer time t f , and also an exponential relation between ∆v 1 and t f . The curve-fitting result with t 1 = 1.10 days is presented in Figure 2. The curves are fitted with linear models, which hold the forms Combining Equations (12), (38) and (39), it can be derived that Equation (40) describes the relation between coefficients a 0 , a 1 and b 0 , b 1 . By MATLAB tool curve-fitting, the linear coefficients are solved: where the units of a 0 and a 1 are, respectively, 1 and 1/day, and the units of b 0 and b 1 are, respectively, kg and kg/day. The goodness of the two models is guaranteed, with an R 2 statistical value of 0.9999. For a set of data y j j = 1, 2, . . . , n y , the estimated value isŷ j , and the average value isȳ. The definition of R 2 is [35] where S = ∑ n y j=1 y j −ŷ j 2 is the sum of the squares of the error, and S yy = ∑ wheret EP represents the estimated minimum time of all-electric transfer. Because the value oft EP is independent of the term related to chemical propulsion, we have Several previous investigations yielded the empirical conclusion that the product of lowthrust magnitude T max and minimum time of all-electric transfer t EP was nearly constant [6,28,36]. There is also a model suggested for the estimation oft EP [35,37] The values of the coefficients vary with the initial and terminal orbits. However, for transfer missions from GTO to GEO and from LEO to GEO, the values of d 1 and d 2 are both near -1, which indicates that the empirical conclusion about the constancy of the product t EP × T max is practical in some cases. Combining Equations (44) and (46), we have the interpolation model in this work: When the coefficients c 1 , c 2 , d 1 , d 2 and d 3 are determined, the fuel consumption ∆m and impulsive burn ∆v 1 can be estimated with transfer time t f , low-thrust magnitude T max , specific impulse of CP I spc and specific impulse of EP I sp .
The transportation rate T R , which is the payload mass advantage gained by using EP divided by the low-thrust transit time [17], can be deduced from Equation (39): Because ∆m c − b 0 = 2.44 kg, the first term on the right-hand side of Equation (48) can be neglected when t f is large enough, and the transportation rate is near constant −b 1 . The actual transportation rate, deduced transportation rate in Equation (48) and simplified quadratic model are presented in Figure 3. For short-duration transfer with several days, the deduced model cannot be applied to estimate the transportation rate. This is because the error is enlarged when the term ∆m c − ∆m is divided by a small value of t f . For hybrid transfer over 6 days, the deduced model is applicable for the estimation of the transportation rate. Transportation rate, kg/day Now, hybrid transfer orbits with transfer time of 40 days, 80 days and 115 days are shown among the transfer orbits with t 1 = 1.10 days. Moreover, the time-optimal transfer orbit with all-electric propulsion is also shown for comparison. The time histories of the semi-major axis, eccentricity, inclination and mass are presented in Figure 4. As is shown in Figure 4, the transfer orbit with t f = 115 days is close to the all-electric transfer orbit. For hybrid transfer orbits with transfer time of 40 days and 80 days, time histories of thrust direction angles in the RSW reference frame are illustrated in Figure 5. The trends of thrust direction and the transfer orbit with transfer time of 80 days are shown in Figure 6. As is presented in Figure 6, the impulsive burn is conducted at the ascending/descending node, which is also the apogee, to raise the perigee and lower the inclination.

Performance of Hybrid Transfer with Various Low-Thrust Maximum Magnitude and Specific Impulse
As is mentioned above, the linear coefficients in Equations (38) and (39) are related to low-thrust maximum magnitude T max and specific impulse I sp . A group of values of T max and I sp are considered to determine the coefficients in Equation (46). The hybrid transfer problems are solved when T max and I sp are set with different values. When T max takes various values, the variation in fuel consumption ∆m with transfer time t f and the variation in the magnitude of the impulsive burn ∆v 1 with t f are demonstrated in Figure 7. When I sp takes various values, the variation in fuel consumption ∆m with transfer time t f and the variation in the magnitude of impulsive burn ∆v 1 with t f are demonstrated in Figure 8. The data are fitted with the linear model presented in Equations (38) and (39). Then, the coefficients a 0 and a 1 can be obtained on different occasions. With the model described in Equation (47) for curve fitting, the coefficients of the interpolation model are shown in Table 3. The goodness of the fit is guaranteed, with an R 2 statistical value of 1.0000. The interpolation surface is illustrated in Figure 9. In these cases, the coefficient a 0 varies between 0.5413 and 0.5453, and its average value is 0.5441. Therefore, the coefficient a 0 is determined as a 0 = 0.5441. Then, the coefficient a 1 can be obtained using Equation (47) (47). Moreover, the value of d 2 is 0.98333, which is close to −1. This suggests that the product (1 − a 0 )/a 1 × T max is nearly constant when I sp is fixed, which has only a small range of variation. This simplified conclusion can be applied for a fast estimation of a 1 with low precision.  Once the coefficients are determined, the model can be applied for estimating the fuel consumption ∆m and impulsive burn ∆v 1 . The hybrid transfer problems with another group of values of T max and I sp are solved for validation. When I sp = 2500 s and I sp = 3500 s, considering different values of T max , the variation in fuel consumption ∆m with transfer time t f , and the variation in the magnitude of impulsive burn ∆v 1 with t f are, respectively, illustrated in Figures 10 and 11. The values of RMSE of the model for estimation are presented in Figure 12. Applying this model for estimation, the RMSE of the fuel consumption is no more than 2 kg, and the RMSE of the magnitude of impulsive burn is approximately 10 −3 km/s.     Among these cases, the case where T max = 100 mN, I sp = 3500 s is discussed as an example. The variation in relative error (∆m − ∆m)/∆m with transfer time t f and the variation in error ∆v 1 − ∆v 1 with transfer time t f are presented in Figure 13. It can be seen that the estimated fuel consumption error is less than 5%, and the estimated error of the impulsive burn's magnitude is within 0.015 km/s. The estimation of the transportation rate is shown in Figure 14. Starting from t f = 10 days, the estimated transportation rates are essentially the same as the real transportation rates. The large error with a small value of t f has been explained before. The relative error of the transportation rate is presented in Figure 15, where the error with t f less than 10 days is ignored. It can be seen from Figure 15 that the transportation error is mainly within 2%. Transportation rate, kg/day

Discussion and Conclusions
In this paper, a trajectory optimization method is proposed for time-fixed minimumfuel orbital transfer with combined chemical-electric propulsion. The first-order necessary conditions for impulsive burn constraints are derived analytically. The long-duration geostationary orbit transfer is many-revolution, and is solved through the homotopic approach from the short-duration transfer problem. For hybrid transfer missions from GTO to GEO, the solution of hybrid transfer with two impulsive burns shows that the magnitude of the second impulsive burn is 0, which indicates that only hybrid transfer with one impulsive burn should be considered.
There are multiple local optimal solutions of this problem. Some of the local optimal solutions have different timings of impulsive burn t 1 . When the transfer time is close to 0, the performance of hybrid transfer is close to that of chemical-only transfer. Moreover, when the transfer time reaches the minimum all-electric transfer time, the performance of the hybrid transfer is close to that of all-electric transfer. The variation in fuel consumption with transfer time is nearly linear, and the variation in the magnitude of the impulsive burn is exponential. The models in Equations (38) and (39) are applied to describe the two relations. According to the linear coefficients a 0 and b 0 , the magnitude of impulsive burn ∆v 1 and fuel consumption at t f = 0 are close to those of chemical-only transfer. The coefficients a 0 and b 0 are constant when the specific impulse of chemical propulsion I spc is fixed. The coefficient a 1 can be expressed with a 0 and the estimated minimum time of all-electric transfert EP . The value oft EP is estimated with an interpolation model. The coefficients of the interpolation model are related to the initial and terminal orbits. Once the coefficients are determined, the model can be applied for estimating the fuel consumption and magnitude of the impulsive burn with a given transfer time, specific impulse of the propulsion system and low-thrust magnitude. A group of hybrid transfer cases with various low-thrust magnitude T max and specific impulse of EP I sp are optimized to offer data for the determination of the coefficients. Then, a new group of hybrid transfer cases are optimized for the validation of the model. The initial mass of the spacecraft is 800 kg, and the specific impulse of CP is fixed with I spc = 300 s. For the cases for validation, the RMSE of the fuel consumption is within 2 kg, and the RMSE of magnitude of the impulsive burn is approximately 10 −3 km/s. For hybrid transfer with T max = 100 mN and I sp = 3500 s, the fuel consumption prediction error is less than 5%. The error of the impulsive burn's magnitude is within 0.015 km/s. Finally, the transportation error is mainly within 2%.
Because this paper is based on some assumptions, the solution of this method is not accurate enough to design a precise transfer trajectory. The results and model can be used to provide an initial solution for existing software or programs such as GMAT, which can be used to further design a precise transfer trajectory.

Conflicts of Interest:
The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.