Design of 3-D trajectory sequences for multiple asteroid flyby missions

Prospecting is a necessary pre-requisite of future asteroid mining ventures. It is generally assumed that inspection for the purpose of identifying the asteroid composition can be effectively accomplished from a distance through remote sensing. To be carried out in a timely and economically viable way, prospecting is best performed by devising trajectories such that a single spacecraft manages to fly by as many asteroid as possible, yet seeking to minimize a cost function that we assume to be coincident with propellant consumption. In this paper, we present a method to identify trajectory sequences to multiple asteroids. We restrict our analysis to Near-Earth Asteroids (NEAs), i.e.,, those with perihelion at less than 1.3 AU from the Sun, focusing on Apollo class NEAs only. The volume of space where encounter seeking takes place is a torus-shaped 3-D region in the proximity of the ecliptic. Under the assumption of using impulsive maneuvers to connect ballistic coast arcs, we show that a deterministic building blocks approach is successful in finding a significant number of multi-flyby mission profiles with the desired characteristics. Using this scheme, it is possible to envisage realistic asteroid prospecting missions using a single launch to deploy a number of small spacecraft, with tens—or possibly hundreds—of asteroids visited in a few years.

problem to find the best combination. Another strategy is the automatic planning and scheduling algorithm introduced by [15] and named Physarum Algorithm (PA), after the unicellular organism Physarum Polycephalum. This method is very reminiscent of B&B algorithms, since it is inspired by the "intelligent" behavior shown by Physarum Polycephalum, but the bounding phase is heuristic in this case. Finally, [16,17] analyze a multiple rendezvous mission to near-Earth asteroids with a solar sail propulsion system.
The MTOP methods and results are very dependent on the related mission concept and assumptions. Indeed, a "standard" form of both problem and resolution does not exist in literature. Most previous studies focus on sequence elements rather than on a mission-level performance figure of merit (such as, e.g., propellant consumption), utilizing heuristic or hybrid algorithms with widely different assumptions and, sometimes, with little or no discussion about trajectory determination. Heuristic methods may also suffer from early convergence, thus not ensuring that the optimal global solution will be actually identified.
In a previous work by our team [18], a 2-D algorithm for multi-target exploration missions at NEAs is formulated, attempting to create a deterministic direct method for asteroid target sequence selection as a function of the mission parameters (e.g., propellant consumption, spacecraft-Sun distance, etc.). A MATLAB ® software package was implemented and successfully tested that provides a useful tool for preliminary analysis of multi-target missions. The reader is referred to [18] for a detailed discussion of the fundamentals and motivations underlying the structure of the algorithm. The present work builds on the method illustrated in [18] to develop a substantially enhanced version of the procedure, with extension of the algorithm to full 3-D trajectories and with detailed modeling of specific phases of the trajectories, such as escape from a parking orbit and asteroid fly-by.
The method here presented does not aim at solving a MTOP, nor does it perform a proper mathematical optimization, other than by selecting among discrete sets of alternative possibilities at various phases in the procedure. This direct procedure, however, allows the mission planner to quickly identify practical, viable trajectories allowing for multiple asteroid flyby's under realistic constraints. In this paper we introduce the algorithm, present several series of tests and discuss the relevant results.

Main assumptions
The algorithm is designed as a tool to find an efficient asteroid visit sequence for a spacecraft equipped with a propulsion system. The goal of the algorithm we seek to build is to determine a trajectory that encounters the largest possible number of asteroids with the lowest possible propellant consumption.
It is understood and accepted that the resulting trajectories are, generally speaking, not the optimal ones in a strict mathematical sense; they are however selected so to meet the needs of the mission planner during preliminary mission design. For this purpose, our earlier 2-D work [18] defined a "mission area" constrained by a minimum and a maximum distance from the Sun (d min and d max ) where the spacecraft's trajectory takes place during a period of time defined by a starting and a final date (t 0 and t end ). That choice was motivated by the fact that most of the asteroids pass through the ecliptic plane at a distance from the Sun between 1 and 1.2 AU, in the so-called high-density transit zone (HDTZ). The spacecraft was bound to move in a 2-D trajectory on the ecliptic plane, trying to intercept the asteroids during their transit. Therefore each asteroid could only be passed by at the nodes of its orbit, thus providing limited opportunities of visiting it.
However, most Apollo asteroids have low orbit inclinations (< 10 deg), so that they stay relatively close to the ecliptic plane during large portions of their orbit. It is therefore reasonable to extend the mission area to a threedimensional volume of space in the vicinity of the plane of the ecliptic. In this work the mission area is thus assumed to be a torus that provides the spacecraft a track on the ecliptic plane from which to depart, or towards which to aim for approach, to reach asteroids. In particular, as described later in Sect. 3, the Sun-centered torus area can be defined by means of Sunspacecraft's maximum and minimum distances, respectively d min and d max (Eq. 5). This feature allows to manage the size of solar system's belt to be explored and to consider the planning of various simultaneous prospecting missions in different areas.
Since we are interested in a multi-target mission, orbital maneuvers will be needed to guide the spacecraft in chasing the asteroids. To minimize the delay between successive asteroid visits, we assume to use impulsive maneuvers only, such as can be provided by a restartable liquid rocket engine. The quantity of propellant needed can be defined in terms of the total velocity change v tot . Since the spacecraft will have to perform multiple maneuvers, it is also reasonable to set a maximum consumption per maneuver v max , so to preserve enough propellant for the final phases of the mission.
The asteroid orbital elements {a, e, i, , ω} with respect to the heliocentric-ecliptic reference frame are assumed known and constant over time. In this scenario, perturbations due to close encounters with planets are neglected. This issue must be taken into account when using the method for the preliminary design of long duration missions; in such case, a parallel check of the selected asteroid trajectories should be performed. The Earth's orbital elements {a ⊕ , e ⊕ , i ⊕ , ⊕ , ω ⊕ } are also known and constant. The initial position of the spacecraft, from which the search for target asteroids begins, is assumed at the Earth.  The algorithm provides the possibility to take into account the escape phase from any initial parking orbit around the Earth. This feature has been introduced for those missions where it is impossible to escape from the Earth by means of the launcher vehicle only. The escape maneuver is assumed impulsive and requires the user to input the parking orbit elements a par , e par , i par , par , ω par with respect to a geocentric-equatorial reference frame. Earth gravity assists, which might indeed be feasible given the near-Earth orbits, are not taken into consideration at this stage.
In summary, the main inputs to the algorithm are as follows: The assumptions discussed so far are summarized in Table 1.

Asteroid flyby analysis
As knowledge of size and mass of asteroids is very scarce, the gravitational effect on spacecraft trajectory during the flyby phase is mostly impossible to calculate. It is therefore interesting to assess whether such effect can actually be neglected. To this end, we take into consideration asteroid 4179 Toutatis as a representative example ( Table 2).
We seek to demonstrate that, for a typical flyby, the gravitational effect of asteroid Toutatis (here assumed to be of spherical shape) on the spacecraft trajectory is negligible. Given the estimated size-frequency distribution of observed NEA's reported in Fig. 1, a similar result will hold for most of Apollo asteroids.
Regardless of the fact that flyby occurs upstream or downstream of the asteroid velocity vector, let v S and v AST be respectively the spacecraft and the asteroid velocities in the heliocentric reference frame (Fig. 2). The spacecraft velocity with respect to the asteroid is, therefore, The quantity v ∞ is the velocity at the border of asteroid's sphere of influence (SOI). Using the superscripts (1) and (2) respectively for entry and exit conditions from the SOI, the module of velocity variation resulting on the spacecraft is: As mass of the spacecraft does not change during the flyby, mechanical energy is conserved. Therefore, the module of the velocity vector relative to the asteroid at the boundary of the sphere of influence is the same both for the incoming leg of the trajectory and for the outgoing leg, with the flyby maneuver resulting in a change of the velocity vector direction only. Let's indicate δ for the flyby turning angle (see Fig. 2); for the conservation of specific mechanical energy we have:  Taking advantage of the properties of a hyperbolic trajectory, it is possible to determine eccentricity and semi-major axis of the flyby trajectory: The maximum v fb is obtained with the largest value of turning angle δ, that is when the pericenter distance of the hyperbolic trajectory is equal to the asteroid radius: Finally, obtaining the sine term from Eq. 2 and replacing it in Eq. 1, we have Equation 4 gives the spacecraft velocity variation as a function of the hyperbolic excess velocity at the border of SOI. To estimate a upper boundary to such change of velocity, let's consider the following function: The denominator of f is always positive, so the function is positive in the domain considered. The derivative of f is The function derivative is positive for 0 < v ∞ ≤ 1/ √ k, hence there is a global maximum for: Using Eqs. 2 and 3: Furthermore, as the hyperbolic excess velocity increases we get: The function f is represented in Fig. 3. As the relative velocity of the spacecraft increases, the gravitational influence of the asteroid vanishes and the trajectory tends to a straight line. Therefore, it is possible to neglect the flyby phase when the initial spacecraft relative velocity v ∞ is much larger than about 10 −2 km/s, as will happen in all practical cases. The flyby phase at the asteroids will, therefore, not be considered in further analysis, with no detriment to the quality of the trajectory determination.

Preliminary phase
Since asteroids will be reachable only within the torus-shaped mission area (Fig. 4), it is necessary to calculate when they enter or leave this domain. To this end, let us consider a parametric definition of the torus surface: where R 1 > 0 is the distance between the torus center and the tube center, while R 2 > 0 is the tube radius.
The torus geometric parameters are related to the maximum and minimum input distances through the following equations: With reference to Fig. 4, let r A = [x A , y A , z A ] be a generic asteroid position vector and r A = x 2 A + y 2 A + z 2 A its modulus. The asteroid is on the surface of the torus when the O z (x, y) Torus cross section passing through the asteroid position with main geometrical parameters indicated asteroid-tube center distance d A is A the Sun-asteroid distance projected on the ecliptic plane. The asteroid is inside the torus when: An asteroid is defined as potentially observable, when Therefore, the asteroid will be selected only if it is already inside the mission area at the time t 0 , or if it is initially external to the torus but its entry is expected at time t < t end ; on the contrary, it will instead be discarded if is initially outside the torus but it will enter at time t > t end , or if it will never enter. Since the torus is not a simply connected domain, it is reasonable to assume that more than one orbital arc exists for which the asteroid is within the torus (usually there are two at maximum). Thus, there are time intervals in which the asteroid is within the mission area, alternating with intervals in which it is out. Determining the transit times t tor and related positions r tor = r(t tor ) on the torus surface involves the solution of a non-linear problem: where the first equation expresses the condition that the asteroid is on the torus surface and the inequalities provide the necessary conditions on the time rate of change of distance from the torus surface. The most direct strategy for deciding if an asteroid meets the condition expressed by Eq. 7 is to find the time of the first transit t (I) tor . For this purpose, one needs to determine all the solutions of Problem 8, r (1) tor , r (2) tor , . . . , and then to solve Direct Kepler Problems (DKPs) for each of them with initial conditions at t = t 0 , so to find the related times t (1) tor , t (2) tor , . . . ; the smallest solution is the next asteroid transit time t (I) tor . Summing up, we have tor,k , . . . for k = 1 : N tot This solution is related to an incoming transit if the asteroid is initially out, or to an outgoing transit if the asteroid is inside. So the first asteroid selection can be performed: After this first selection, the number of asteroids taken into consideration in the following phases will be: • Potentially observable asteroids number: N ≤ N tot 4 Search phase

Initial assumptions
Once the asteroids within the reach of the mission have been identified, attention is focused onto the target search phase. The goal is to derive a strategy to be used at any time during the mission to determine the next destination to visit and the relevant transfer orbit. The strategy must depend only on the current state of the spacecraft-i.e. mission time t, spacecraft position r S (t) and velocity v S (t)-and so a local selection of asteroids is made. Consequently, the mission profile will be shaped step-by-step with each new choice. The spacecraft final trajectory will be a segmented curve composed by transfer orbit arcs from an asteroid to the next one. For a n-asteroids sequence with indexes {k 1 , k 2 , . . . , k n } encountered at times {t 1 , t 2 , . . . , t n }, we have: The search begins with the spacecraft positioned at the Earth: Propellant consumption for the escape maneuver will be determined later either by considering escape from the parking orbit (see Sect. 5) or by leaving it to the launcher to provide the necessary impulse: • First maneuver consumption v 1 = 0 with launcher v esc with escape phase.
As discussed in Sect. 2, flybys must take place within the mission area and the selected targets must be reached with a propellant consumption lower than the pre-defined maximum ( v < v max ). There may be situations in which multiple asteroids fulfill this requirement; in such case, selecting the target accessible with the lowest v is not necessarily a wise choice, as this does not necessarily lead to maximize the number of asteroids encountered during the mission. Instead, the preferred trajectory will be determined using a tree-structured strategy ( [18]). During each search every selected asteroid is first taken into account and then new target searches are performed starting from each one of them. In this way we obtain a complex branching of possible paths that the spacecraft can travel. The branching will be stopped at those nodes where either the propellant is exhausted, or the end of mission time has been exceeded, or simply no further reachable targets are found. For a n-asteroids sequence with consumption values { v 1 , v 2 , . . . , v n }, we have the following criterion: • Tree-structure branching stop criterion Finally, a global selection of the best mission profile is made by selecting the sequence of branches that, starting from the Earth, passes by the largest number of asteroids with the lowest propellant consumption.

Asteroid transit observation area
To reduce the computational cost, the search is initially focused only on one part of the mission area. Since the mission has a finite duration, it is preferable to intercept the asteroids as close as possible to the spacecraft, to minimize transfer times between two targets and thus to maximize the number of asteroids encountered. To this end, we consider a spherical observation area centered at the spacecraft's position within which the flyby must take place. The sphere must be large enough to hold a number of possible encounters and, at the same time, small enough to limit the transfer distance.
To determine the sphere's radius δ, we take as a reference a characteristic length of the system-namely, the distance between the spacecraft and the currently closest asteroidand multiply it by a magnification factor c 1 chosen at will: • Observation area radius where t * is the current mission time. In our case c 1 = 10 is assumed.
It is necessary to select only those asteroids that transit within the observation area and to calculate when this happens. Let r A be a generic asteroid position; the asteroid is on the surface of the sphere when the distance from the spacecraft d A is Therefore, it is inside the sphere when: An asteroid must be intercepted when it is simultaneously inside the observation area (sphere) and the mission area (torus): In the above formula, an interval of time equal to the asteroid period T has been set to include all the positions along a single orbit. Let S and T be respectively the sphere and torus three-dimensional domains; their intersection is defined as D = S ∩ T . An asteroid will be selected only if it is already inside D at the time t * , or if it is initially external but its entry is expected at time t < t * + T ; it will instead be discarded if is initially outside D but it will enter at time t > t * + T , or if it will never enter. Determining the transit times t D and positions r D (t D ) on D involves the solution of a non-linear problem: To determine if an asteroid respects the condition in Eq. 10 it is necessary to find the time of the first transit t The obtained solution is related to an incoming asteroid if it is initially out of D, or, on the contrary, to an outgoing asteroid if it is inside. Finally, the transiting asteroid can be selected: • Transiting asteroids selection After this step, the number of asteroids considered could decrease again: • Transiting asteroids number: N 1 ≤ N .

Flyby position
In this phase the previously calculated transit times t D,k , . . . are used to set the temporal domain I in which the asteroids are reachable. Among them, an optimal time for flyby is selected, if possible. Finally, all the data concerning transfer orbits and maneuvers are calculated.
The transit time intervals are found by re-ordering the transit times in the form t D,k , . . . , and by considering the current asteroid position r * . In particular, two different cases must be distinguished: (1) when the asteroid is initially inside the domain D, and (2) when it is initially outside: • Transit time intervals statement To minimize the transfer times, the search for optimal flyby conditions starts with the first of the aforementioned time intervals and moves on to the following ones if needed. To determine a transfer trajectory between two positions in a given time interval, a Lambert Problem (LP) must be solved. Let us consider the spacecraft position at current time r * S = r S (t * ) and the asteroid position at a generic flyby time r f = r(t f ). The goal is to determine t f ∈ I to find transfer trajectories from r * S to r f in a time interval equal to ( t) tr = t f − t * . Therefore, following [18], the following problem is solved: • Asteroid flyby position search The problem above and then a LP are solved for each asteroid, but only those that admit a solution are considered for treestructure branching, so that the nodes in the tree-structure include only the reachable asteroids: • Reachable asteroids number: N fin New search phases will be carried out for each of them considering the spacecraft state as shown in Sect. 4.1. On the other hand, if no asteroids can be reached according to the defined criteria, the size of the observation area is increased and the search phase restarts: • Observation area increase: In our case, c 2 = 2 is chosen; by considering a larger domain, the chances of finding favorable flyby positions are increased. The upper limit 2d max represents the maximum distance between two points inside the torus mission area. If increasing the observation area radius does not allows to find reachable asteroids, the tree-structure branching is interrupted. h Line of nodes The escape trajectory lies on the plane identified by the position over the parking orbit r S where the maneuver occurs, the position of the Earth and the direction of the vector v ⊕ ∞ (Fig. 5). In particular, it is possible to have two escape maneuvers on this plane which are distinguished by the orbit travel direction-i.e. clockwise or anti-clockwise-and the true anomaly swept for escaping. To indicate the orbit travel direction, the angular momentum unit vector is calculated for both cases: From now on the discussion will be valid for both cases, therefore no distinction will be made. First, the orbital elements of the escape hyperbola are determined. The escape orbit semi-major axis can be determined through conservation of mechanical energy: The escape orbit eccentricity can be calculated by solving a non-linear problem. Let's consider the true anomaly ν swept during escape: This must be equal to the difference between the hyperbolic excess and the initial true anomaly, respectively ν ∞ and ν esc : In Eq. 13, the escape eccentricity e esc is the only unknown term: The remaining orbital elements are determined using their definition. Note that the unit vectors refer to a geocentricequatorial reference system: With the escape orbit parameters known, it is now possible to calculate the associated propellant consumption. Let v par and v esc be respectively the parking and escape orbit velocities in r S ; the escape maneuver consumption is Therefore, the consumption for escape depends on two variables: the escape position r S and the escape orbit travel direction (± in Eq. 12). To optimize v esc , let's consider the following two functions: The two functions are optimized independently of each other; at the end, the best value between the two optimal ones is chosen.

Numerical solvers
To solve the non-linear problems introduced in the previous sections two solver algorithms have been implemented, one for root-finding problems and the other for optimization problems, both derived from [20]. These algorithms are designed for single-variable functions and are derivative-free. Contrary to most iterative numerical methods, they are designed so to be both very convergent and very robust. Brent's algorithms are a combination of a robust and a convergent method (see Table 3) that exploit the good properties of both to guarantee superlinear convergence independently from the starting values. Let f (x) be the single-variable function inherent to the problem. There are four parameters that need to be defined as inputs to Brent's algorithm. The tolerance tol is a combination of a relative tolerance ε and an absolute tolerance t: where x * is the known best approximation of solution. The parameter ε is the machine precision-equal to 2 −23 or 2 −52 for a 32-bit or a 64-bit processor, respectively-while t must be arbitrarily assumed to have the desired tolerance. In our case, we set t = 10 −8 . The remaining input parameters are the search interval edges, a and b, such that the solution is x * ∈ [a; b]. To identify them, (N + 1) equidistant control points x i (for i = 0 : N ) are taken within the f (x) domain (edges included) and the function is evaluated for each of them. The search interval edges are chosen using the following criteria: • Optimization problem As N increases, both accuracy and numerical effort grow. In our test cases (Sect. 7) it was assumed N = 100.

Test results
The software has been tested in two different ways to highlight how the trajectory depends on the mission start date. The first case includes a series of tests on 2-year missions spread over a period ranging from 2020 to 2040, while the second one includes 2-year missions in the period 2020 to 2025. In both cases only the mission time interval is changed, while the remaining input parameters are constant through the various executions; the escape maneuver is not considered. We focus here on the results obtained in a simple, constrained scenario. Investigation of the algorithm dependence on the input parameters is postponed to future works. Because of the limited computing power available and to prevent the generation of very large tree-structures with unacceptably long processing duration, the input values to the test case were arbitrarily selected to correspond to a 2-year mission with a modest quantity of propellant onboard ( Table 4). The maximum maneuver consumption was chosen to obtain a significant, albeit small, number of encounters, i.e., 16 in our test case.
To get an idea of the amount of data produced, the following parameters have been recorded and compared for each test: Testing was carried out on a computer equipped with an Intel ® Core ™ i7-3630QM processor. Software was compiled with the GNU FORTRAN compiler (GFortran), a free open source FORTRAN 95, FORTRAN 2003 and FORTRAN 2008 compiler. In all cases, each test run took not more than a few tens of minutes to complete. The results are reported and discussed in the following subsections. More complex mission cases can be addressed with acceptable computational time (tens of minutes or a few hours at most) by running the algorithm on present commonly available multi-core personal computers.

Test case 1
In this series of tests, twenty software runs were performed for mission epochs ranging from January 1, 2020, to January 1, 2040. The main results are shown in Table 5.
As can be seen from Fig. 6, the number of possible flybys has no clear dependence on the mission period. However, it can be noted that over a period of 20 years there are numerous opportunities for prospecting missions to nine or more asteroids, many of them even with the same launch. In particular, for a mission that starts on January 1, 2040, it is possible to have five different trajectories each of which intercepts twelve asteroids, for a total of 60 asteroids encounters in two years. The "best" of these, i.e. the one that exhibits the smaller total propellant consumption, is shown in Fig. 7 and the relative dates and consumptions are listed in Table 6. Figure 7 shows the projection of the 3-D trajectory on the plane

Test case 2
Also, in this case twenty tests were carried out, for a period ranging from January 1, 2020, to January 1, 2025. The main results are shown in Table 7.
The maximum number of asteroids found in a trajectory and the number of sequences with this feature are shown in Fig. 8 as a function of the mission period. Unlike case 1, the frequency of negative results is higher, probably due to the shorter time span between missions.
However, there are some cases of significant relevance for practical application to prospecting: e.g., it appears to be possible to perform a mission with nine spacecraft, each  Time is 12:00:00 UTC for all epochs directed towards twelve asteroids, for a total of 108 asteroids prospected in two years. Considering a more basic mission architecture with the use of a single spacecraft, it is possible to fly a trajectory passing by 15 asteroids (Fig. 9, Table 8). Like Fig. 7, also Fig. 9 shows the projection of the trajectory on the ecliptic plane, the out-of-plane displacement being too small for the scale of the drawing.

Conclusions
As expected, the 3-D algorithm yields better results than its 2-D predecessor [18], albeit at an increased computational cost due to the need to solve a larger number of non-linear problems. For the same v tot and the same mission area size, the 3-D version returned a flyby frequency of about 6 asteroids/year against the 4 asteroids/year obtained with the 2-D version. The number of possible trajectories with the same launch is also higher. In many cases, launch windows were found in which it was possible to travel more than 10 different trajectories with 10 asteroids each. We conclude that the procedure here presented provides a relatively simple The first consumption is zero due to the use of the launcher to provide propulsion for the escape. All times are UTC method to determine trajectory sequences that allow for a substantial number of flyby's. Identifying asteroids that are potentially rich in raw materials requires a considerable effort. It has been estimated ( [21]) that, if harvesting costs are taken into account, the cost of prospecting must be a small fraction (10%) of the material value in order for the mining venture to generate a profit. Therefore, considering an ideal asteroid with a size of 100 m rich in platinum group metals for an estimated value of US$ 5B, the prospecting mission cost should not exceed US$ 500M; this strongly suggests that small, relatively cheap satellites should be seriously considered for the task.
The results of this study show that a multi-asteroid prospecting mission with a fleet of small satellites is indeed an attractive possibility. According to [21], finding at least one rich asteroid with a probability of 99% requires the investigation of about 44 asteroids. In light of the results of this work we can confidently conclude that such target can be achieved with a reasonable number of small spacecraft, thus keeping launch cost and duration of the prospecting operations at a minimum.
Author Contributions All authors contributed to the study conception and design. The first draft of the manuscript was written by Giuseppe Cataldi and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Funding The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.