Minimum-Fuel Planar Earth-to-Mars Low-Thrust Trajectories Using Bang-Bang Control

: Recent advance in electric propulsion systems have demonstrated that these engines can be used for for long-duration interplanetary voyages. Constant speciﬁc impulse engine described as a thrust-limited engine is an example of this type of engine, processing the ability to operate at a constant level of impulse. The determination of minimum-fuel, planar heliocentric Earth-to-Mars low-thrust trajectories of spacecraft using a constant speciﬁc impulse is discussed considering the ﬁrst-order necessary conditions derived from Lawden’s primer vector theory. The minimum-fuel low-thrust Earth-to-Mars optimization problem is then solved in two-dimensional, heliocentric frame using both indirect and direct methods. In the indirect method, two-point-boundary-value problems are derived to solve boundary value problems for ordinary differential equations. In the direct method, a general-purpose optimal control software called GPOPS-II is adopted to solve these optimal control problems. Numerical examples using two different optimization methods are presented to demonstrate the characteristics of minimum-fuel planar low-thrust trajectories with on-off-on thrust sequences at three chosen ﬂight times and available maximum powers. The results are useful for broad trajectory search in the preliminary phase of mission designs.


Introduction
Electric propulsion systems have demonstrated that low-thrust engines have the capability to be used for long-duration travels by the planetary and interplanetary space missions. Electric propulsion has been used by NASA's Deep space 1 [1] and Dawn [2], ESA's SMART-1 [3], and JAXA's Hayabusa and Hayabusa 2 [4]. Electric propulsion systems have demonstrated that lowthrust engines can be used for long-duration travels by the planetary and interplanetary space missions [5]. Low-thrust electric propulsion spacecraft is known to have a greater payload capability than conventional chemical propulsion spacecraft. Low-thrust trajectory optimization generally accompanies determining the control variables, which may include thrust magnitude, and direction and parameter, and the corresponding trajectories while minimizing a given performance index (propellant consumption or time-of-flight) and satisfying boundary conditions (departure and arrival orbits), mid-point conditions, and path constraints.
Low-thrust trajectory optimization problems can be generally solved by either indirect or direct method [6]. The optimization of low-thrust trajectories have been mathematically formulated as an Optimal Control Problem (OCP) [7]. Indirect methods solve the optimal control problem by obtaining the solution to the corresponding two-point boundary value problem (TPBVP) which results from the calculus of variations. The application of the Lagrange multipliers (costates) doubles the number of differential equations that have to be integrated along the trajectory. However, the solution to the TPBVP is very sensitive to initial guess for costate variables which do have any intuitive physical meanings. To overcome the disadvantage in the indirect methods, heuristic and/or evolutionary techniques have been developed [8,9]. In contrast, direct methods solve the optimal control problem by converting it into a nonlinear programming (NLP) problem with various transaction schemes applied to either states and controls, or both states and controls. Since the control variables are explicitly parameterized, a good initial control guess for direct methods can be easily produced. In direct methods, the modifications of the performance index, equality constraints, state and control inequality constraints can easily made for different problem formulations while they should result in a new derivation of TPBVP in indirect methods. This paper presents indirect and direct methods for obtaining fuel-optimal planar Earth-to-Mars trajectories using CSI engines where the differential equations of motion are numerically integrated and the continuous-time control is parameterized. Thus, the purpose of this paper is to develop an indirect and direct optimization methods to solve minimum-fuel planar low-thrust Earth-to-Mars trajectories in the heliocentric frame with bang-bang control structure. In addition, the solutions of this problem using the indirect and direct methods will be validated by comparing them at the different flight time and maximum powers of the spacecraft and analyzing the optimization results. A performance index representing propellant consumed for CSI engines is formulated to minimize its total propellant consumption, resulting in on-off-on thrusting sequence by Lawden's primer vector theory [10,11]. Minimizing the performance index is equivalent to maximizing the final mass of the spacecraft. The formulation of the problem treats the spacecraft mass as a state variable, thus updating the spacecraft mass to the optimal trajectory design. For both optimization methods, the equations of the motion are normalized for the fundamental distance, velocity, mass and time to streamline the numerical computations. In the the indirect method, the analytical formulations of the problem are presented to set up TPBVP using the primer vector theory, and the necessary conditions for an optimal solution are discussed. In the direct method, the bounds of the states, the flight time, and the control including the maximum available power, the equality constraint and the boundary conditions are explicitly specified. The optimal control problem is then converted to the parameter optimization problem that can be solved by nonlinear programming (NLP). As the NLP solver, a general-purpose optimal control software called GPOPS-II [12] is adopted to solve optimal control problem using variable-order Gaussian quadrature collocation methods where the continuous-time optimal control problem is transcribed to a large sparse nonlinear programming problem.
Numerical studies show that the fuel-optimal, low-thrust heliocentric Earth-to-Mars trajectories for the specific arrival time and maximum power are obtained with different thrust magnitudes to find "on-off-on" thrusting sequence. The primer vector theory is employed to analyze the control structure by monitoring the variations of the switching functions. Finally, the fuel-optimal planar Earth-to-Mars trajectories at different flight time and maximum power are validated.

Problem Formulation
In this section, the governing equations of motion are given and the scaled equations of motion are derived. The spacecraft is propelled by electric propulsion.

Equations of Motion
This trajectory optimization problem is modeled in a two-dimensional, heliocentric (suncentered) polar coordinate system to make the NLP problem both efficient and robust instead of Cartesian coordinates which is the simplest but most disadvantageous choice [13]. All motions are assumed to be confined to the ecliptic or fundamental plane. Figure 1 illustrates the geometry of this coordinate system along with the steering angle.
where r is the heliocentric radius of the spacecraft, θ is the phase angle, u and ν are the radial and transverse velocity components, respectively, and α is the steering or thrust orientation angle. The steering angle is measured relative to the instantaneous tangential direction and is positive in the clockwise direction to the thrust vector. The two-dimensional equations of motion for the spacecraft in the polar coordinate system [13,14] are given bẏ where m is time varying spacecraft mass, u 1 and u 2 are the unit vector components in the thrust vector, µ is the solar gravitational constant, which equals to 1.32712441933 × 10 11 , T represents the thrust magnitude and c = I sp g 0 is the exhaust velocity. I sp and g 0 = 9.80665 m/s 2 are the specific impulse and the standard gravitational acceleration at sea level, respectively. The radial and tangential acceleration components due to thrust are defined ad The boundary conditions and the inequality should be specified for a complete optimization problem.
Here, the initial boundary conditions can be formulated as while the terminal boundary conditions are given by The first equation in Eq. (5) states that the radial velocity should be zero and the second equation is a boundary condition is that forces the final velocity to be equal to the local circular velocity at the final radial distance. The equality path constraint is given by

Scaled Equations of Motion
We scale the state variables to remain close to unity to streamline numerical computations. All heliocentric distances are normalized with respect to the AU (Astronomical Unit) which is equal to 149597870.691 km. Likewise, all velocity values are normalized with respect to the "local circular velocity", υ θ 0 = µ/r 0 at the heliocentric distance of the Earth's circular orbit, r 0 . Using the initial orbit radius and initial velocity as reference values, and the normalizing the time by the final time, we define the new state variables as: Applying the chain rule to change variables results in the following new non-dimensional equations of motion: For r : For υ r : where a = 2ηP mc t f υ θ 0 .
For υ θ : For θ : Form : In summary, newly scaled equations of the motion are: The scaled initial boundary conditions can be formulated as while the scaled terminal boundary conditions are given by

Low-Thrust Trajectory Optimization
The trajectory is controlled by the thrust magnitude and direction. The thrust level is obtained by selecting exhaust velocity and the engine input power P (thruster efficiency η is assumed constant throughout), taking the operational constraints int account. The input power is always limited by the availability of the solar power (P ≤ P max ) and the exhaust velocity is constant. The payload is the performance index to be maximized. For the minimum-fuel problem with CSI engines, the performance index can be established as It is obvious from Eq. (15) that minimizing J is equivalent to maximizing the final mass m f or the payload.
where λ υ r υ θ = [λ υ r , λ υ θ ] T and u = [u 1 , u 2 ] T . The optimal control law is derived to minimize the Hamiltonian. Thus, the direction of the thrust is opposite to the adjoint vector λ υ r υ θ and the primer vector p(t) is defined, so that λ T υ r υ θ u = −λ υ r υ θ = −p(t).
for which where p(t) = −λ υ r υ θ is the primer vector and λ υ r υ θ = λ 2 υ r + λ 2 υ θ . Thus, the unit vector u is given by The terms in the bracket of the Hamiltonian in Eq. (17) is written as H using a in (9) During a constant I sp operation where c is constant, the input power, or equivalently, is the only control. Equation (21) is usually rewritten as and the sign of the switching function S F determines when the thruster is turned on or off. The choice of the input power, P, that minimizes the Hamiltonian in Eq. (21) is then given by the bang-bang control law: The maximum available power is adopted when S F < 0, whereas the engine is switched off when S F < 0 to minimize H , according to the PMP. In addition, the thrust magnitude for 0 ≤ T ≤ T max will also depend on the algebraic sign of the switching function S F . The Euler-Lagrange equations yield: From the known initial conditions in Eq. (12) with the known initial time t 0 = 0 s and the known final conditions in Eq. (13) with the known final time t f , we can get 10 boundary conditions. We apply the transversality condition to pose a well-defined TBPVP.
Expanding the nonzero terms in Eq. (26) while noting that we have a Mayer form of the performance index, where φ = J = −m f : Using the the terminal boundary conditions, we get the following boundary conditions: All differential equations and boundary conditions give the well-defined TPBVP:

Direct Method
The same minimum-fuel problem is solve with the direct method. Minimize the performance index in Eq. (15) subject to the dynamics constraints which is the scaled equation of motion in Eq. (11) The continuous-time state variables y and control variables u are given, respectively as The equality path constraint function and boundary condition function are given as The right-hand side function of the dynamics, the path constraint function, and the boundary functions are written, respectively, as

a(y(t), u(t),t) =
Finally, the lower and upper bounds on the path constraints and boundary conditions are all zero. Finally, the lower and upper bounds on the path constraints and boundary conditions are all zero. Because the first five boundary conditions, (b 1 , ..., b 7 ), are simple bounds on the initial and final continuous-time state, they will be enforced in the NLP as simple bounds on the NLP variables corresponding to the initial and terminal state. The 7th boundary condition, b 7 , on the other hand, is a nonlinear function of the terminal state and, thus, will be enforced in the NLP as a nonlinear constraint.

Numerical Results
In this section, the mission scenarios of the planar Earth-to-Mars orbit transfer at five different input powers and the flight time are given to validate minimum-fuel low-thrust trajectory optimization with the constant specific impulse engines. Both the indirect and direct methods are used to solve the same problem and their solutions are compared. The initial and final orbital parameters are displayed in Table 1. It is observed that the final position in the mission orbit (Mars orbit) is not specified since the purpose of the simulation is to reach the orbit not the planet. As the forces acting on the spacecraft, the Sun's gravity and the thrust produced by the engines are considered. The spacecraft is assumed to have initial mass 1500 kg and the CSI engines with the I sp = 3300 s and the constant thruster efficiency, η = 0.7. The used input powers and the corresponding flight time are listed in Table 2. The input powers and flight time conditions of the scenario 1, 2 and 3 are listed in Table 2. The scenario 1 has 19 kw of input power and 240 days of flight time, the scenario 1 has 7.5 kw of input power and 365 days of flight time, and the scenario 3 has 3.6 kw of input power and 2 years of flight time. In these scenarios, the indirect method solutions are obtained by solving the well-defined TPBVP in Eq. (29) with bvp4c function in MATLAB while the direct method solutions are obtained with GPOPS-II in the Section 3.2. In addition, SNOPT [15] is used as a NLP solver in GPOPS-II.  Figures 2, 3, and 4 present the results of a transfer from Earth to Mars orbit in the heliocentric reference frame in astronomical units (AU) using the indirect and direct methods for the scenarios 1, 2 and 3, respectively. The dashed inner circle represents the Earth's orbit while the dashed outer circle represents the Mars' orbit. The curves represent the spacecraft trajectory while the arrows represent the thrust unit vectors. Since the indirect method uses 10000 mesh points, the thrust unit vectors in the Earth-to-Mars trajectories are not obviously indistinguishable. On the other hand, the direct method obvious thrust unit vectors along the trajectories because the number of nodes adjusted in the the optimization calculation. However, the Earth-to-Mars trajectories using the indirect method and direct methods are matching with no remarkable difference. These numerically matching results are verified in the control acceleration and the thrust angle as shown in Figs. 5, 6 and 7. In these results, the control accelerations by the thrust on-off-on thrusting sequences are obtained by the ban-bang control law (23) for the CSI engines. The thrust on-off-on thrusting sequences are also analyzed using the primer vector theory, which appear in the next figures. The control accelerations and thrust angle (α) are also almost matching except the thrust off phases. It is due to some numerical difference in the indirect and direct methods. However, the numerical differences in the results are not remarkable enough to result in different optimization results. Figures 8 and 9 present time histories of thrust and acceleration magnitudes, and time histories of the propellant mass and total mass using the indirect and direct methods for scenario 1. Numerically matching results are also obtained for scenario 1. Table 3 shows the ∆v and the propellant masses for scenario 1, 2 and 3, respectively using the indirect and direct methods. The numerical integration of the control acceleration is computed as ∆v using   Figures 10, 11 and 12 present the switching functions and thrust profiles for scenario 1, 2 and 3, respectively. They show on-off-on thrusting sequences corresponding to the bang-bang control law in (23). The thrust magnitude switches its limiting values of 0 (null-thrust arc) and T max (maximum-thrust arc) each time S F (t) passes through 0 according to Eq. (22). The left side of Figs. 10, 11 and 12 are obtained in the indirect method while the right side of them are obtained in GPOPS-II. In general, the costate vectors are computed in indirect method. However, GPOP-II can provide the computed costate vectors in the optimization algorithm. With these costate vectors, the switching function S F and the thrust magnitude (T ) are also obtained in the right side of Figs. 10, 11 and 12. The switching functions and thrust profiles using the indirect and direct methods for scenario 1, 2 and 3 present matching results. Thus, minimum-fuel planar Earth-to Mars low-thrust trajectories using bang-bang control for the CSI engines have been validated by both the indirect method and direct method. In addition, it has been demonstrated that ∆v and the propellant mass m p can be saved in longer flight time and smaller input power through three selected mission scenarios.

Conclusion
In this paper, minimum-fuel planar heliocentric Earth-to-Mars low-thrust trajectories using bang-bang control for constant specific impulse engines of spacecraft have been studied. For the problem formulation, the derivation of the scaled equations of motion is described in detail. A well-defined TPBVP is derived for the indirect method while GPOPS-II is emoyed to solve the problem. Using the three selected mission scenarios with different flight time and input powers, minimum-fuel planar low-thrust trajectories with on-off-on thrust sequences are validated with both the indirect and direct methods. Numerically matching optimization results are obtained by both methods with the same on-off-on thrust sequence. The results are useful for broad trajectory search using the CSI engine in the preliminary phase of mission designs.