Multi-objective low-thrust spacecraft trajectory design using reachability analysis

One of the fundamental problems in spacecraft trajectory design is ﬁnding the optimal transfer trajectory that minimizes the propellant consumption and transfer time simultaneously. We formulate this as a multi-objective optimal control (MOC) problem that involves optimizing over the initial or ﬁnal state, subject to state constraints. Drawing on recent developments in reachability analysis subject to state constraints, we show that the proposed MOC problem can be stated as an optimization problem subject to a constraint that involves the sub-level set of the viscosity solution of a quasi-variational inequality. We then generalize this approach to account for more general optimal control problems in Bolza form. We relate these problems to the Pareto front of the developed multi-objective programs. The proposed approach is demonstrated on two low-thrust orbital transfer problems around a rotating asteroid. © 2022 The Author(s). Published by Elsevier Ltd on behalf of European Control Association. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Since the Galileo mission in 1991 we have seen a steady increase in proposed missions to asteroids and comets, as they might hold the key to many scientific questions including the origins of life on earth [3]. The Dawn mission to Vesta and Ceres proved the viability of low-thrust electric propulsion for asteroid exploration [1,32], and it is expected that many upcoming missions will rely on similar low-thrust propulsion. While there has been a significant study of interplanetary transfer trajectories using low-thrust propulsion, comparatively little research has been conducted on the trajectory design in the vicinity of asteroids. We investigate a spacecraft trajectory design problem around an asteroid, where the objective is to use minimal amounts of propellant to raise an orbit while keeping flight times as short as possible. This is a multi-objective optimal control (MOC) problem, whereby one seeks to find the optimal way a dynamical system can perform a certain task, while minimizing or maximizing a set of, usually contradictory and incommensurable, objective functions [30]. Conventional optimization techniques for spacecraft trajectory design often fall into two categories, indirect methods, based on the calculus of variations, and direct methods, whereby the optimal control problem is reformulated as a nonlinear program. Direct approaches rely on parametrization and while a candidate solution is found, there are no guarantees on the optimality of the solution. Indirect methods, meanwhile rely on necessary analytic conditions for optimality using Lagrange multipliers. Yet while optimality of the obtained solutions may be guaranteed, indirect approaches, such as multiple shooting methods, rely on a good initial approximation of the optimal trajectory [39,41]. A third approach is dynamic programming, whereby the optimality conditions are formulated in continuous time based on the Hamilton-Jacobi-Bellman (HJB) equation [39], however, it is hampered by the so-called curse of dimensionality. Despite this, unlike direct approaches, optimality is guaranteed, and unlike indirect approaches, the solution does not rely on an initial approximation of the optimal trajectory. We expand on this third approach by taking advantage of recent developments in reachability analysis.
Reachability analysis aims to find the set of points from which a target can be reached within a given time, subject to constraints. It forms a fundamental part of the dynamics and control literature and has been used extensively for controller synthesis of complex systems [4,24,29]. In recent years we have seen considerable research being conducted into computing reachable sets using Hamilton-Jacobi (HJ) reachability analysis, whereby the reachable set is derived from the viscosity solution of a HJB equation accounting also for the presence of state constraints. Such a HJB framework is presented in [11,25] with more general value problems bypassing previous regularity issues presented in [2]. In [17] an extension of the HJB framework to time-varying targets and constraints is considered and in [15] the approach is extended to multi-objective control problems. HJ reachability has also been successfully applied to various aerospace applications including air traffic control [26], the climbing problem of multi-stage launchers [8], payload optimization [9], as well as most recently to the complete model of the ascent problem of multi-stage launchers [7]. One of the advantages of using HJ reachability is that the optimal trajectory can easily be constructed once the reachable set has been computed. This makes HJ reachability attractive for problems that require computing trajectories for various different initial states.
For the spacecraft trajectory design problem considered in this paper, there are two possible formulations for minimizing the burnt propellant. The first assumes that the initial mass is a free optimization variable. This approach is common during mission design where the total required fuel budget is being calculated. To this end, we formulate the spacecraft trajectory design problem as a MOC problem and show that it can be equivalently stated as an optimization problem subject to a constraint that involves the sub-level set of a certain value function. The latter is shown to be the unique continuous viscosity solution of a quasi-variational inequality that involves a HJB equation. Such value functions have been defined in [11,25] to account for the presence of state constraints. This formulation allows characterizing the Pareto front of the formulated MOC problem and also facilitates its computation by means of available numerical tools.
The second formulation of the spacecraft trajectory design problem assumes a fixed initial mass and the objective is to maximize the remaining mass after completing a given orbital maneuver. This approach is more common when the maneuver needs to be added to a given mission and the available fuel is nonnegotiable. This formulation had been previously investigated by the authors in [40]. To solve this second formulation of the spacecraft trajectory design problem, we draw on research from [15] and extend our formulation of the MOC problem to introduce an auxiliary state, allowing us to account for arbitrary problems in Bolza form, and, together with appropriate normalization and approximations allowing for a reduction of the state space, greatly improving on the method presented in [40].
Thus our contributions can be summarized as 1. the formulation of an efficient constrained MOC problem for low-thrust spacecraft trajectory design that optimizes only over the set of admissible initial states and transition times, 2. the reduction of the state space through the use of appropriate approximations, 3. the expansion of the MOC problem to allow for a generalization of the proposed methodology for arbitrary multi-objective problems in Bolza form.
This paper is organized into six sections. Section 2 contains details regarding the derivations of the spacecraft dynamics as well as the definitions of the constraints pertaining to its behavior. In Section 3 the optimal control problem is formulated while Section 4 describes how the set of admissible initial states is derived from the viscosity solution of a quasi-variational inequality. Section 5 is dedicated to the numerical computation and case study of an orbital transfer around a rotating asteroid. Finally, Section 6 provides concluding remarks and directions for future work.

Spacecraft equations of motion
We begin by modeling the dynamics of the spacecraft. The spacecraft thrust is defined in spherical coordinates as where α(t) ∈ [−π, π] is the incidence angle, δ(t) ∈ [− π 2 , π 2 ] is the sideslip angle and T(t) ∈ [0, T max ] is the variable thrust, with T max denoting the maximal allowable thrust. The Cartesian transformation of the thrust vector is denoted by u x (t), u y (t) and u z (t), respectively. The compact set U = [−π, π]×[− π 2 , π 2 ]× [0, T max ] is the set of possible control input values while u ∈ U ad denotes the control policy and U ad denotes the set of admissible policies which is the set of Lebesgue time measurable functions from [−∞, 0] to U. Note that time here is considered to be non-positive to facilitate the reachability problem exposition in Section 3. Boldface notation is used to denote time varying functions such as trajectories and policies, while non-boldface notation is used to denote scalars and vectors. The equations of motion of the spacecraft around a rotating body can be expressed in 3-dimensional Euclidean space as a second-order ordinary differential equation (see eg., [21]) where R(t) is the radius vector from the asteroid's center of mass to the particle, the first and second time derivatives of R(t) are with respect to the body-fixed coordinate system, U(R(t)) is the gravitational potential of the asteroid and Ω is the rotational angular velocity vector of the asteroid relative to inertial space. The term 2Ω(t) × dR(t) dt describes the Coriolis forces, Ω(t) × (Ω(t) × R(t)), the centrifugal forces and dΩ(t) dt × R(t) the Euler forces. We consider an asteroid rotating uniformly with constant magnitude ω around the z-axis. Therefore, the Euler forces can be neglected and we can express the rotation vector as Ω := ωe z , where e z is the unit vector along the z-axis. Following [18], the radius vector and its derivatives are given by The Coriolis and centrifugal forces (the first two terms in (2)) acting on the spacecraft are thus To model the current position, velocity, and available propellant, we define the state vector where ∆m ∈ R ≥0 denotes the available propellant. The total spacecraft mass can be expressed as m(t) = m 0 + ∆m(t), where m 0 denotes the dry mass of the spacecraft. Following our derivations from (2), we can formulate the dynamics of the spacecraft,ṙ = f (r, u), as where v exhaust ∈ R ≥0 is the exhaust velocity used to express the depletion of mass as propellant is burned, U x , U y and U z are the derivatives of the gravitational potential in the direction of the unit vectors e x , e y and e z , respectively, and where for brevity we neglect the time dependence by denoting r = r(t), v x = v x (t), and similarly for the other states.

State constraints
Since the dynamics of the spacecraft were derived for orbits in the vicinity of the asteroid, we need to enforce state constraints on x, y, z. We naturally also need to ensure that we bound the amount of propellant available. Assuming that the burnout mass of the spacecraft is the same as the dry mass, we set m min := 0 and m max := m propellant and impose m min ≤ ∆m ≤ m max .
Due to particles ejected from the asteroid, we do not want to fall below a circular orbit with radius ρ := x 2 + y 2 + z 2 of approximately ρ min = 1 km. Furthermore, in order for the two-body problem under discussion to be valid and the influence of other bodies in the solar system to be negligible, we need to stay within the sphere of influence (SOI) of the asteroid. The SOI can be approximated as in [37] by ρ SOI ≈ a M1 M2 2 5 , where a is the semi-major axis of the asteroid's orbit around the sun (1.5907 · 10 8 km), M 1 is the Mass of the asteroid (1.4091 · 10 12 kg) and M 2 is the mass of the sun (1.9890 · 10 30 kg). Therefore, the sphere of influence of the asteroid is approximately ρ max = ρ SOI ≈ 8.74 km. The set of states that satisfy the aforementioned restrictions is given by The target orbit that we would like to transfer to is denoted by the closed target set C ⊂ K.
The initial orbit that we start at is denoted by the closed initial set I ⊂ K. Note that the initial and the target orbit restrict only the position and the velocity, but allow the mass to take any admissible value within [m min , m max ].
While Cartesian coordinates are useful for modeling the behavior of an object around a rotating body, since we restrict all admissible states to lie within the set K, which constrains the radius ρ, it is more efficient to recast our problem in spherical coordinates. To this end, define a ρ , a θ , a ψ as the transformations of The tangential velocity in the x-y plane, v t , and its perpendicular counterpart, v ⊥ , can then be defined as follows: Then we can restate the system dynamics in spherical coordinates as With a slight abuse of notation, we also redefine I, C, and K in spherical coordinates. Finally, we impose the following assumptions on the spacecraft dynamics.
Due to norm equivalence, the choice of norm is irrelevant and not further discussed. Using Assumptions 2.1 and 2.2, for any control policy u ∈ U ad , any initial state r 0 ∈ K and transfer time t f > 0, the system admits a unique, absolutely continuous solution on [−t f , 0] (see [34]).

Multi-objective optimal control problem
Having defined the system dynamics, we are now in a position to discuss how to find trajectories that start on an initial orbit, I, and take the spacecraft to some final orbit, C. Additionally, the objective is to keep the flight time and required propellant as small as possible. Thus, the multi-objective optimal control problem can be formulated as a minimization problem whereby the first goal is to minimize the required propellant, ∆m, and the second is to minimize the required time for the orbit change, i.e., the transfer time, denoted by t f . The trajectory, r, which is the solution of (11), belongs to the Sobolev space . The set of trajectory-control pairs on [−t f , 0] starting at r 0 with transfer time t f is denoted as: Note that as in [14] we adopt the convention that 0 denotes the terminal time hence the transfer time t f denotes the time duration. Under Assumption 2.1 and by Filippov's Theorem [23, pg.121], we can conclude, that Π r0,t f is compact.
Remark 3.1. For similar applications as the one discussed in this paper, where Assumption 2.1 might not hold, we refer to [7] where convexification of the dynamics is considered in order to ensure that the set of absolutely continuous solutions of the problem is closed.
The set of admissible (in the sense of satisfying the state constraints) trajectorycontrol pairs on [−t f , 0] starting at r 0 with transfer time t f is denoted as: Finally, the set of admissible initial state and transfer time pairs is denoted as For a given initial state r 0 ∈ R 7 and transfer time t f ∈ [0, +∞), we can define the cost functions as J 1 (r 0 , t f ) := ∆m and J 2 (r 0 , t f ) := t f , where ∆m is the 7-th element of the state vector r 0 . The 2-dimensional objective function J : We are now in a position to formulate the multi-objective optimal control problem under study as minimize

Pareto optimality
The solution of (13) in general does not consist of a single isolated point, but rather a set of optimal compromises between the objectives J 1 and J 2 [27].
where a vector a is considered less than b (denoted a < b) if for every element a i and b i the relation a i < b i holds. The relations ≤, ≥, > are defined in an analogous way. Following Definition 3.2, a solution (r 0 , t f ) is considered Pareto optimal if it is not possible to improve all its performance metrics J 1 (r 0 , t f ), J 2 (r 0 , t f ) simultaneously. The set of Pareto optimal solutions is called the Pareto set P S , while its image is the Pareto front P F . Therefore, the solution of (13), i.e., the set of minimizing (r 0 , t f ) pairs, is the desired Pareto set, while the cost function corresponding to the minimizing (r 0 , t f ) pairs is the Pareto front.
To allow for mission designers to determine a compromise between minimizing required propellant and transfer times, we wish to compute the Pareto front. However, the unconventional constraint in (13) ensuring a solution (r 0 , t f ) is feasible, prevents us from solving (13) with standard MOC problem solvers. Therefore, we will next discuss how we can recast the constraint (r 0 , t f ) ∈ Γ to a standard nonlinear inequality constraint, which will allow us to compute the Pareto front by means of conventional MOC problem solvers.

Solution to multi-objective optimal control problems
To find an equivalent formulation for the constraint (r 0 , t f ) ∈ Γ, let g(r) and ν(r) be two Lipschitz functions (with Lipschitz constants L g and L ν , respectively) chosen such that g(r) ≤0 ⇐⇒ r ∈ K, ν(r) ≤0 ⇐⇒ r ∈ C. This can be achieved by choosing g(r) and ν(r) as the signed distance to the set K and C, respectively.
Next, we consider the value function ω: where a b denotes max(a, b). We are now in a position to use the value function to decide if, for a given initial state and transfer time, there exists a corresponding admissible trajectory. Thus we can introduce an equivalent formulation of (13).
Theorem 4.1 implies that the Pareto front can be computed from the solution of (15). To achieve this we discuss how to compute ω.

Value function computation
To begin to discuss how ω can be obtained, we introduce the Hamiltonian H : where q ∈ R 7 is the costate vector.
Since the Dynamic Programming Principle [6] holds the proof of Theorem 4.2 follows standard arguments for viscosity solutions, as shown in [2,25]. Note that the infimum should be understood to be over the restriction of Π r0,t f over In order to solve the quasi-variational inequality in Theorem 4.2, we employ a finite differences scheme. As in [11,25], a consequence following from Theorem 4.2 is the Lipschitz continuity of the value function. The proof is similar to the more general proof of Proposition 4.9 that is introduced in the sequel. Proposition 4.3 allows us to make statements about the discrete-continuous error estimate, for which we refer to Theorem 5 in [11], as well as the convergence of the value function, for which we refer to Proposition 6 in [16].
The Hamiltonian admits an explicit form. To this end, consider the term Then we can write the Hamiltonian as As T is always positive, the thrust angles can be optimized separately; see Appendix for more details. To this end, the derivation of α * and δ * follows a similar procedure as in [7] and is listed in the Appendix. After applying the optimal thrust angles, the Hamiltonian becomes affine in T and we are able to find the optimal thrust magnitude. In particular, H(r, q) becomes = 0, is negligible for the consideration of the optimal Hamiltonian. For the optimal control computation during the trajectory calculation, we have numerically investigated the occurrence of singular arcs and found no instances in which √ q 2 4 +q 2 5 +q 2 6 m0+∆m = − q7 v exhaust over an extended interval, thus we do not further consider the singular control case.
Finally, applying T * and rewriting the minimum as the maximum of the negation of the associated function, the Hamiltonian takes the following analytic form H(r, q) = −C(r, q) + max q 7 T max v exhaust + T max m 0 + ∆m q 2 4 + q 2 5 + q 2 6 , 0 .

Extension to problems in Bolza form
We will now generalize our approach to problems where the objective functions do not rely only on the initial state and are written in Bolza form. To achieve this we need to introduce auxiliary states. As in [2,15], we show how problems in Bolza form are reformulated into Mayer form, and then show how problems in Mayer form are solved in a similar fashion as in Section 3. Consider the p-dimensional objective function defined as: where J t denotes the terminal cost and J r denotes the running cost. We impose the following assumptions, as in [15].  Remark 4.7. To ease notation, we omit the dependents on the neighborhood for the Lipschitz constants and instead assume the existence of a global Lipschitz constant L t and L r , respectively.
Next, we define an auxiliary state z ∈ R p as: where z 0 becomes an optimization parameter and z ∈ W 1,1 (R p ). The auxiliary state captures the cumulative running cost and thus is treated as an additional state. In the same manner we previously ensured a trajectory, r, stayed within the set K, we bound z and ensure that the integrated running cost, added with the terminal cost, stays below some value z 0 . To capture all possible trajectories, we introduce the set and make the following assumption.
Assumption 4.8. For every r ∈ R 7 ,the set is a compact convex subset of R 7 × R p .
We now introduce the auxiliary value function ϑ: g(r(s)) , (23) where i x i denotes the maximum element of the vector x. As with ω, the term max s∈[−t f ,0] g(r(s)) and ν(r(0)) ensures that any trajectory r remains in K and terminates in C. The additional term J i t (r(0)) − z i (−t f ) ensures that the integrated running cost, J r , combined with the terminal cost, J t , never grows larger than z 0 . Thus, in addition to ensuring that (r, u) are admissible trajectory control pairs, the sub-zero level set of ϑ bounds the terminal and integrated running cost. Therefore, We are now in a position to introduce the generalized multi-objective optimal control problem for objective functions in Bolza form: where z 0 represents an upper bound for the term J Bolza (r, u, t f ), without explicit knowledge of r or u. The generalized value function can again be obtained as the unique continuous viscosity solution of a quasi-variational inequality where the Hamiltonian is defined as Proposition 4.9. The value function ϑ is Lipschitz continuous.
The proof can be found in the Appendix. Proposition 4.9 can be used to show that a numerical solution of (26) (in the viscosity sense) can always be determined. Under Assumptions 2.2, 4.5, 4.6 and 4.8, by Filippov's Theorem [23], the problem (23) admits an optimal solution, which implies the existence of an admissible r and z [15]. This yields the following relationship due to (21)

Numerical Approximation and Results
We will now discuss how the value functions can be obtained numerically, prior to discussing how the spacecraft trajectory design problem is solved. Following Proposition 4.3, a numerical solution to (26) can be found. To this end, we employ the Level Set Methods toolbox of [28]. For the computation of ω, we use a Lax-Friedrich Hamiltonian where p + and p − are the right and left derivatives computed using an appropriate fifth-order weighted essentially non-oscillatory (WENO) scheme. The Lax-Friedrich Hamiltonian consists of an analytic expression of the Hamiltonian (derived previously), as well as a dissipation term, which is scaled by the dissipation coefficients α k . The dissipation coefficients α k needs to satisfy Since we need α k to serve as an upper bound, we consider the control input that maximizes the Hamiltonian: For a further discussion of the Lax-Friedrich Hamiltonian and WENO scheme, we refer to [31], while for a discussion of the convergence of ω and the derivation of a necessary Courant-Friedrichs-Lewy condition, we refer to [11,19,28].

Implementation
To illustrate the theoretical results of the previous sections, we consider a spacecraft on an initial near circular orbit around asteroid Castalia 4769. The goal is to compute an efficient transfer trajectory that raises the spacecraft to

Scale
Values Distance Initial radius ρ 0 Velocity 2 m/s Time Maximum thrust T max a stable orbit at an altitude of 6117.5 m above the asteroid. For the derivation of a stable orbit around Castalia 4769 we refer to [22] and references therein. The gravity of Castalia 4769 was modeled by means of a spherical harmonic expansion as discussed in [20,35,36]. Even though the proposed theoretical framework allows us to tackle problems of any state dimension, the available numerical tools and computational power limit us to only study the lowerdimensional planar case of the application. We, therefore, omit the states ψ and v ⊥ to consider only To avoid ill-conditioning when solving the HJB equation, the state vector is normalized using the constants introduced in Table 1. This results in the following dynamics: where c = T max ρ 0 /(m 0 V 2 0 ) is a normalization constant. ρ 0 , m 0 , and v 0 denote the initial radius, mass, and velocity of the initial orbit, respectively. The optimal control policy and trajectory (r, u) ∈ Π K,C r0,t f can be constructed efficiently using the numerical approximation of ω. For a given N ∈ N we consider the time step h = 1 N and a uniform grid of [−t f , 0] with spacing s k = k N . Let us define the state {r k } N k=0 and control {u k } N −1 k=0 for the numerical approximation of the optimal trajectory and control policy. Setting r 0 as the initial orbit, we proceed by iteratively computing the control value u k (r k ) ∈ arg min u∈U ω(r k + hf (r k , u), s k ) g(r k ).
For a given ω(r k , s k ), this is done by numerically taking the partial derivatives along each grid direction to estimate the costate vector, q, and then determining the optimal control values as the minimizer of the Hamiltonian, H. After u k is determined we compute r k+1 using the Matlab ode113 function, a variableorder Adams-Bashforth-Moulton method of order 1 to 13 [38], and increment k. For the implementation, we discretized the interval [−t f , 0] using N = 4000 grid points. As shown in Figure 1, when considering orbits further than 4 km away from the surface of the asteroid, the variation of the gravitational acceleration along θ becomes negligible. It is, therefore, possible to approximate the gravitational terms in spherical coordinates as Using this approximation makes a ρ and a t independent of θ. This allows us to omit a grid dimension while numerically solving the quasi-variational inequality in Theorem 4.2, greatly reducing the computational cost. Thus the final set of states used for the computation of the value function is During the calculation of the optimal trajectory, θ can easily be reconstructed by forward integrating the dynamics at each time step s k , i.e.,

Simulation results
The spacecraft is modeled with 750 kg of dry mass, 600 mN of maximum thrust, and an exhaust velocity of 40 km/s. Using an initial orbit with radius 5.1 km and tangential velocity of −2.4 m/s, we are able to compute the numerical approximation of ω using the spatial grid described in Table 2 in combination with a temporal grid using N = 1000 grid points. The propagation of the zero level set over time is shown in Figure 2 and 3. Since the temporal grid for the calculation of the value function is more coarse than that used for the trajectory calculation, we need to interpolate the value function while computing the final trajectory. The final computed trajectory, for an initial propellant mass of 24.89 g and transfer time 2757 s is shown in Figure 4. The asteroid rendering for Figure 4 was computed as in [22].
The accuracy of the final orbit is within 41 meters of the target orbit. Calculating ω took 9 hours using a 3 GHz 8-Core Intel Core i7-9700 processor running Matlab with an extension of the Level Set Methods toolbox [13]. The run time and accuracy can be significantly improved upon when using optimized code  such as [10,12]. To show how the accuracy of the solution, the memory usage as well as the CPU time varies with the grid size, we recompute the value function and the trajectory with coarser spatial grids, as presented in Table  3. We use only single-precision arrays to store the value function, yet utilized double-precision arrays for all numerical calculations on the value function. The subscript final denotes the values of the final point of the computed trajectory, i.e. r(0), while target refers to the target orbit used to define C and subsequently initialize the computation of the value function. Once ω is computed, it is incorporated into (15), which is solved using Matlab's paretosearch function.
Solving the MOC problem took 120 seconds and the resulting Pareto front is shown in Figure 5. A comparison of the thrust magnitude of the smoothed optimal control policy is shown in Figure 6. As can be seen, the time-optimal solution uses near    Figure 4 is derived using the point marked by the black diamond, while the time-optimal solution used in Figure 6 uses the point denoted by the magenta pentagram.  continuous thrust to reach the target, at the cost of using a large amount of fuel. The control policy of the trajectory presented in Figure 4 meanwhile, has noticeable cruising phases where no fuel is consumed. As expected, the thrust magnitude of both policies follows a bang-bang structure.
To illustrate the results of problems in Bolza form, we consider the case of optimizing the remaining propellant in oppose to the initial propellant. Therefore, let us fix the initial propellant to 100 g. We use a similar setup as in 5, with the addition of the auxiliary state, z, defined as a terminal cost. We consider the uniform spaced grid over r and z, defined in Table 4. Storing the value function using single precision requires 15.2 GB of memory. Since the objective is to maximize the remaining propellant, the optimization problem needs to minimize −∆m. Therefore, for a given final state r f ∈ R 7 and transfer time t f ∈ [0, +∞), we define the cost functions as J 1 (r f , t f ) := −∆m and J 2 (r f , t f ) := t f , where ∆m denotes the 7-th element of the state vector r f (the mass in our case). The 2-dimensional objective function J : R 7 × [0, +∞) → R 2 can then be written as  The resulting Pareto front is shown in Figure 7. Calculating the reachable set took approximately 32 hours of CPU time. Using the reachable set, calculating the Pareto front took approximately 258 seconds. As expected, the Pareto front looks similar to that of Figure 5, yet not identical due to numerical inaccuracy and change in the initial mass, resulting in modified dynamics.
Using an initial propellant mass of 100 g, z of −75 g, and transfer time of 2698 s, the transfer orbit is computed in the same way as before. As expected, the transfer trajectory is similar to the trajectory shown in Figure 4, and the accuracy of the transfer orbit is within 66 meters of the target orbit. For comparison to the first approach, the accuracy of the second approach is shown in Table 5. From comparing the Pareto front, it can be seen that both methods have a minimum transfer time of around 2500 seconds as well as a minimum propellant requirement of just over 20 g.

Conclusion
We have presented a novel method of using the value function of a quasivariational inequality to compute the decision space of multi-objective optimization problems. The feasibility and effectiveness of the proposed approach was demonstrated by applying it to the problem of low-thrust trajectory design. The approach is applicable to arbitrary multi-objective optimization problems where the control variable is required to lie within a reachable set.
Since the Hamiltonian becomes affine with respect to the thrust magnitude, once the thrust angles have been fixed, future research aims to exploit this fact by integrating classification based approaches [33]. Furthermore, utilizing approximations of the reachable set as in [5], as well as decomposing the reachable sets as in [14], seems promising.
The proof for ω is similar to that of ϑ and will therefore be omitted in the interest of space.