Accessibility of near-Earth asteroids and main-belt asteroids in a gravity-assisted multi-target mission

Accessibility of near-Earth asteroids and main-belt asteroids in a gravity-assisted multi-target mission

Planetary and Space Science 182 (2020) 104851 Contents lists available at ScienceDirect Planetary and Space Science journal homepage: www.elsevier.c...

4MB Sizes 0 Downloads 60 Views

Planetary and Space Science 182 (2020) 104851

Contents lists available at ScienceDirect

Planetary and Space Science journal homepage: www.elsevier.com/locate/pss

Accessibility of near-Earth asteroids and main-belt asteroids in a gravity-assisted multi-target mission Pan Sun 1, Hongwei Yang 2, Shuang Li *, 3 College of Astronautics, Nanjing University of Aeronautics and Astronautics, Nanjing, 210016, China

A R T I C L E I N F O

A B S T R A C T

Keywords: Near-earth asteroid Main-belt asteroid Asteroid accessibility Gravity assist Tisserand graph Trajectory optimization

To maximize the science return of asteroid exploration, multi-target asteroid missions are preferred. In this paper, the accessibility of asteroids in a mission, including both near-Earth asteroid (NEA) sample return and main-belt asteroid (MBA) rendezvous using gravity assists (GAs) is investigated. 1081 NEAs and 37 MBAs are chosen as candidates with given constraints. A Tisserand graph with modified tick marks is applied to give a straightforward and quick reference of the possible gravity-assist sequences. The patched conic method is adopted to compute the total velocity increments for six types of gravity-assist trajectories. The particle swarm optimization (PSO) algorithm is employed to minimize the total delta-V. To reduce the computational complexity of the optimization problem, only 50 optimal asteroid groups attained under the mission path of Earth-NEA-EGA1(the first Earth gravity assist)-MBA (ENEM) are reserved to be further optimized for the remaining other five trajectory paths. Feasible trajectory design results with launch dates between 2022 and 2028 are presented.

1. Introduction Asteroids are thought to be the primitive bodies in our solar system. They hold the significant clues for explaining both the formation of planets and the origin of life. While some asteroids could pose potential risks for the Earth, rich resources can be exploited from them in the future space missions (see (Green et al., 2013)). The asteroids with perihelia within 1.3 AU of the sun are known as near-Earth asteroids (NEAs) (see (Cheng et al., 1997)). Of all the asteroids, most are found in the main asteroid belt between the orbits of Mars and Jupiter, known as main-belt asteroids (MBAs) (see (Cheng et al., 1997)). According to the planetary science division of the NASA science mission directorate, missions are generally ranked in terms of science return as a flyby, rendezvous, land, rove, and ultimately performing a sample return (see (Dankanich et al., 2010)). Apart from the huge science return expected, it is not that difficult to accomplish the sample return from NEAs using current state-of-art chemical propulsion or electric propulsion (see (Dankanich et al., 2010)). Hayabusa accomplished the world’s first sample return mission (see (Kawaguchi et al., 2008)), while Hayabusa2 (see (Tsuda et al., 2013)) and OSIRIS-REx are in the process of a round-trip NEA exploration mission. Although there have been plenty

of documented studies focusing on NEAs, a bounty of science exists in the main asteroid belt (see (Dankanich et al., 2010)). The Galileo makes the first encountering of two MBAs (see (Chapman, 1994)), while Dawn accomplishes the first mission to reach two MBAs (see (Garner et al., 2013)). Furthermore, the ZhengHe mission (see (Zhang et al., 2019)) of China aiming at returning a sample from the Earth’s quasi-satellite 2016HO3 and rendezvousing with the main-belt comet 133P/Elst-Pizarro, is scheduled to be launched in 2022. Low-cost and feasible asteroid exploration missions with both a scientific and a technological character are favored in the deep space exploration. Extensive research findings considering asteroid exploration accessibility can be found in literatures. In (Lau and Hulkower, 1989) the authors presented a ranking by accessibility of NEAs for ballistic rendezvous and sample return trajectories. In (Qiao et al., 2006) the authors introduced the technique of EGAs to further reduce the launch energy and total delta-V. In (Perozzi et al., 2001) the authors provided basic targeting strategies by considering both scientific and technical information. In (Cui et al., 2010) the authors proposed a complete set of approaches to selecting mission targets and designing the transfer trajectory. Considering the dense asteroid population in the main asteroid belt and the potential more expanded exploration (see (Yen, 1982)), MBA

* Corresponding author. E-mail address: [email protected] (S. Li). 1 Ph.D. candidate, Advanced Space Technology Laboratory, and College of Astronautics. 2 Associate professor, Advanced Space Technology Laboratory, and College of Astronautics. Member AIAA. 3 Professor, Advanced Space Technology Laboratory, and College of Astronautics. https://doi.org/10.1016/j.pss.2020.104851 Received 26 September 2019; Received in revised form 2 December 2019; Accepted 21 January 2020 Available online 23 January 2020 0032-0633/© 2020 Elsevier Ltd. All rights reserved.

P. Sun et al.

Planetary and Space Science 182 (2020) 104851

authors apply a machine learning algorithm to the MBA accessibility assessment. To reduce the cost of a single launch and increase the science benefit for asteroid exploration, multi-target asteroid exploration missions have significant roles. In (Sears et al., 2004) the authors investigated the

exploration missions also greatly aroused researchers’ interest. In (Yen, 1982) the author gave candidate targets and detailed feasible MBA exploration mission options for the early 1990s using chemical propulsion. In (Chen et al., 2014) the authors further investigated the accessibility of MBAs via gravity assists. Recently, in (Shang and Liu, 2017) the 0.7

18

0.9<=a<0.95

0.95<=a<1.0

1.0<=a<1.05

0.9<=a<0.95

1.05<=a<1.1

0.95<=a<1.0

1.0<=a<1.05

1.05<=a<1.1

16

0.6

14

0.5

inclination(deg)

eccentricity

12

0.4

0.3

10 8 6

0.2 4

0.1 2

0 0.9

0.95

a)

1 semi-major axis of 1081 NEA(AU)

1.05

0 0.9

1.1

0.95

Semi-major axis vs eccentricity

b)

1 semi-major axis of 1081 NEA(AU)

1.05

1.1

Semi-major axis vs inclination

Fig. 1. a) Semi-major axis vs eccentricity, and b) Semi-major axis vs inclination for target 1081 NEAs. Table 1 Orbital elements of candidate MBAs [MJD: 58600]. MBA_num

Name

a/AU

e

i/deg

Ω/deg

ω/deg

M/deg

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37

Flora Harmonia Melpomene Victoria Euterpe Vesta Urania Nemausa Thyra Athamantis Iris Metis Ausonia Nausikaa Massalia Nysa Hebe Lutetia Fortuna Isis Parthenope Julia Amphitrite Astraea Aspasia Egeria Irene Pomona Thalia Fides Eunomia Io Proserpina Klotho Juno Bamberga Angelina

2.2017642 2.2672466 2.2966536 2.3343151 2.3466645 2.3614179 2.3655722 2.3658354 2.3795825 2.3823032 2.3853338 2.3856366 2.3954309 2.4031006 2.4097818 2.4234511 2.42516 2.4349582 2.4427107 2.4417036 2.4531094 2.5500122 2.5541136 2.5742489 2.5755453 2.5759811 2.5855674 2.5885078 2.6256835 2.6422834 2.6441003 2.6517654 2.6543309 2.6684905 2.6691496 2.6814471 2.6811249

0.1564993 0.0470213 0.2176744 0.2201716 0.1732259 0.0887214 0.127581 0.0675594 0.1931027 0.0615837 0.2312057 0.1231143 0.1269055 0.2450422 0.1420667 0.1479229 0.203007 0.1633561 0.1580469 0.2228231 0.1004723 0.184748 0.0726956 0.1910946 0.0725909 0.0851214 0.1665823 0.0809015 0.2349839 0.1756051 0.1860843 0.1945458 0.0901448 0.2571373 0.2569423 0.3406546 0.1255294

5.88696 4.2574 10.12873 8.37307 1.58371 7.14177 2.09575 9.97718 11.59296 9.45096 5.52365 5.57682 5.77558 6.79799 0.70875 3.70655 14.73791 3.06406 1.57378 8.51462 4.62989 16.1296 6.08252 5.36699 11.26816 16.53613 9.12165 5.52394 10.11426 3.07081 11.75243 11.96206 3.56342 11.7776 12.98892 11.1017 1.3097

110.88929 94.18878 150.38387 235.41018 94.78784 103.81084 307.46872 175.9785 308.79564 239.85332 259.56321 68.90857 337.73327 343.10857 206.10893 131.54757 138.6402 80.86565 211.14406 84.19731 125.5466 311.54939 356.34176 141.57661 242.15281 43.2219 86.12268 220.44169 66.84679 7.26904 292.93432 203.11353 45.77783 159.62429 169.85274 327.85103 309.11789

285.28745 269.73164 227.95081 69.64178 356.44995 150.72843 87.42605 2.58053 97.1419 140.12731 145.26509 6.41736 295.94408 30.67492 256.77317 343.33264 239.80747 249.91704 182.06496 237.30543 195.55035 45.38133 63.36319 358.68757 353.22363 80.54479 97.85894 339.53569 60.63681 62.84353 98.49868 123.03247 193.44993 268.70452 248.13861 44.28052 178.80564

194.88292 107.525 267.25439 133.33592 335.31609 95.862 29.09799 96.04493 16.74347 23.21224 140.41971 276.86164 139.55636 358.19327 117.69516 82.70872 86.19795 344.29353 197.33866 203.41429 278.93071 137.45193 284.23536 282.3663 139.30645 187.48857 164.93587 31.92099 335.82021 341.54257 283.38773 220.68087 105.72362 310.99343 34.92504 89.85985 38.70045

2

P. Sun et al.

Planetary and Space Science 182 (2020) 104851

maneuvers, the modified Tisserand graph for preliminary sequence analysis and the brief problem statement. Section 3 presents the optimal delta-V searching method for the six types of gravity-assist trajectories. In Section 4, feasible example mission lists obtained by PSO with launch dates between 2022 and 2028 are presented and discussed. Global optimal results for the mission of 2014 JR24 sample return and Urania rendezvous are analyzed in detail. Section 5 gives the final conclusion.

sample return of three NEAs for the Hera mission. In (Morimoto et al., 2004) the authors focused on two mission types for the “post-MUSES-C” asteroid exploration mission, one being the multiple rendezvous and sample return mission, and the other being the multiple flybys and sample return mission. In (Yamakawa et al., 2004) the authors studied the sample return missions of multiple NEAs under the solar electric propulsion for the orbital sequences both with and without EGA. In (Olympio, 2011) the authors obtained the formulation of the optimal control problem for the multiple asteroid encounter problem considering both flyby and rendezvous conditions. In (Song and Gong, 2019) the authors solve the multiple near-Earth asteroids rendezvous problem using the deep neural networks for solar sails. Also note that of the 10 Global Trajectory Optimization Competitions (GTOC) organized to date (see (Ceriotti and Vasile, 2010)), five are concerned about multiple asteroid exploration missions, such as multiple asteroid rendezvous, multiple asteroid sample return and multiple asteroid flyby, etc. In particular, in (Yang et al., 2015) the authors studied the low-cost transfer between a near-Earth asteroid and a main-belt asteroid using low thrust propulsion and multiple gravity assists. To solve the multi-target asteroid exploration mission, multiple gravity assists (see (Bayliss, 1971), (Sims, 1997), (Vasile and Pascale, 2006) and (Ceriotti and Vasile, 2010)), sequence optimization methods such as branch and bound (see (Yang et al., 2018)) and beam search (see (Wilt et al., 2010)), Tisserand graph (see (Strange and Longuski, 2002)), evolutionary algorithms including particle swarm optimization (PSO) (see (Kennedy and Eberhart, 1995), (Pontani and Conway, 2010) and (Bao et al., 2015)) algorithm and modified Ant Colony Algorithm (MACO) (see (Ceriotti and Vasile, 2010)), etc., are common techniques utilized in preliminary trajectory design. Moreover, some machine learning algorithms such as Gaussian process regression (see (Shang and Liu, 2017)) and deep neural networks (see (Song and Gong, 2019) and (Li et al., 2019)) provide fast optimal estimation solutions. Among them, multiple gravity assists help reduce total launch energy or transfer time in missions to high delta-V targets. The Tisserand graph is a graphic technique able to identify all feasible ballistic gravity-assist sequences to a target destination (see (Strange and Longuski, 2002)). The PSO algorithm is a very efficient, reliable and intuitive methodology (see (Pontani and Conway, 2010)) for global optimization. The accessibilities of the MBAs and the NEAs via gravity assists have been studied by (Chen et al., 2014) and (Qiao et al., 2006), respectively. The software that performs the NHATS (Near-Earth Object Human Space Flight Accessible Targets Study) round-trip trajectory calculations for NEAs also has been implemented at JPL/CNEOS.4 However, to the authors’ knowledge, few studies have been published on the accessibility analysis of asteroids in a mission including both NEA sample return and MBA rendezvous. Although (Yang et al., 2015) focus on the specific case of designing low thrust trajectories between a NEA and a MBA, detailed accessibility analysis of NEAs and MBAs is not presented. Different from their research works, this paper evaluates the accessibility of asteroids in a mission including one complete NEA sample return (by EGA1) and one complete MBA rendezvous using chemical propulsion, in which the scenario at Earth departure for the rendezvous of the NEA is considered, and then successively come the operations of the NEA rendezvous and sampling, the NEA sample return by EGA1, and the MBA rendezvous. The simplified two-body dynamical model capable of providing adequate preliminary design results is used. Furthermore, a Tisserand graph with modified tick marks is introduced to give a straightforward and quick reference of the feasible gravity-assist sequences. Considering the complexity of the optimization problem, the efficient PSO algorithm is employed to find the global optima for different mission paths. The paper is organized as follows. Section 2 gives a basic formulation of the problem, including the selection criteria for candidate NEAs and candidate MBAs, the model of gravity assist and impulsive delta-V

4

2. Problem formulation and analysis 2.1. Candidate NEAs and MBAs When it comes to the selection of target asteroids, both science value and technology feasibility need to be considered. In this paper, the semimajor axis, eccentricity, and orbital inclination of candidate NEAs and MBAs are chosen within the following constraints in the heliocentric ecliptic inertial frame 8 < 0:9 AU < ¼ aN <¼ 1:1 AU e < 0:6 : N iN < 15∘

(1)

8 < 2:0 AU < ¼ aM < ¼ 2:7 AU D > 100 km : M iM < 18∘

(2)

where 1 AU is 1 astronomical unit (1 AU ¼ 1.49597871  108 km). aN, eN and iN are the semi-major axis, eccentricity and orbital inclination of the NEA, respectively. aM, DM and iM are the semi-major axis, diameter and orbital inclination of the MBA, respectively. When the NEA is in the vicinity of Earth, the sample return mission can be completed at a lower cost. Restrictions on the orbital inclination can reduce the energy required to perform out-of-plane maneuvers, thus further reducing total fuel consumption. The constraint of eN < 0.6 can avoid the target NEA orbit being too oblate and hence reduce the NEA rendezvous difficulty. MBAs with diameters greater than 100 km are selected in this paper since they are considered as candidate targets for future asteroid missions. They were chosen as candidate targets for accessibility analysis of MBAs in (Chen et al., 2014) as well. Considering the mission cost, it is much easier to accomplish the rendezvous of the inner belt and the middle belt asteroid than the outer belt asteroid. Therefore, only the rendezvous missions of the inner and the middle belt asteroid are considered in this paper. Under the constraint conditions (1) and (2), 1081 candidate NEAs and 37 candidate MBAs are identified. The orbital characteristics of all

v v

vp

v

P

v Fig. 2. The reference frame of the gravity assist impulsive model.

https://cneos.jpl.nasa.gov/nhats/. 3

P. Sun et al.

Planetary and Space Science 182 (2020) 104851

3.5

3.5

Outer belt

Outer belt

3

3

Middle belt

Middle belt

2.5

2.5

Inner belt

Ra /AU

Inner belt

Ra /AU

Urania: 2.5 km/s

Urania: 1-7 km/s

2

2 Mars: 3 km/s

Mars: 1-9 km/s

1.5

1.5 Earth: 3.5 km/s

2014 JR24: 2 km/s

2014 JR24: 1-7 km/s 1

0.5

0

0.5

2014 JR24: 1 km/s

1

Earth: 1-7 km/s

Earth: 1 km/s

1

1.5

2

0.5

2.5

0

0.5

1

Rp /AU

1.5

2

2.5

Rp /AU

a) The v contours for 2014 JR24, Earth, Mars, and Urania

b) Illustration of the sequence analysis result

Fig. 3. a) The v∞ contours for 2014 JR24, Earth, Mars, and Urania, and b) Illustration of the sequence analysis result.

3

3.5 v P

2.5

+

M=

M0 - 2 Mmax

Outer belt 3

vp

Middle belt 2.5 2

Urania: 2.5 km/s

E0 ,

Inner belt

Ra /AU

Ra /AU

M= M0 - Mmax

M0

Mars: 3 km/s

2 Mars: 3 km/s

1.5

1.5 Earth: 3.5 km/s Earth: 3.5 km/s

1

1

Earth: 1 km/s

E= E0 + Emax

0.5

0.5 0.6

0.7

0.8

0.9

1

1.1

1.2

1.3

1.4

1.5

1.6

0

0.5

1

1.5

2

2.5

Rp /AU

Rp /AU

a) The v

2014 JR24: 2 km/s 2014 JR24: 1 km/s

contours with modified tick marks

b) Illustration of the modified sequence analysis result

Fig. 4. a) The v∞ contours with modified tick marks, and b) Illustration of the modified sequence analysis result.

v0

vn1 , vn 2

vega1T

vf

Earth

NEA

EGA1

MBA

selected NEAs (all retrieved from the Near-Earth Objects Dynamic Site5) are shown in Fig. 1. The orbital elements of all selected MBAs (all retrieved from the International Astronomical Union6) are given in Table 1. 2.2. Gravity assist and impulsive delta-V maneuvers

t (1)

0 Lam

1 Lam

2 Lam

t (2), t (3)

t (4)

3

4

The simplified impulsive model of gravity assist is considered in this paper. It is assumed that the spacecraft is only subject to the central gravitational force of the sun when it is outside the influence sphere of

t (5)

5 6

Fig. 5. Patched conic method for the ENEM trajectory. 4

https://newton.spacedys.com/neodys/[retrieved on 12 April 2019]. https://minorplanetcenter.net/iau/MPEph/[retrieved on 16 April 2019].

P. Sun et al.

Planetary and Space Science 182 (2020) 104851

Table 2 Global minimum J for the ENEM trajectory [J < 6.77 km/s]. N_50

NEA_name

NEA_num

MBA_num

C3(km2/s2)

J (km/s)

J00(km/s)

Tof(day)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

2004 2014 2012 2010 2010 2007 2012 2012 2008 2014 2008 2010 2013 2012 2011 2000 2013 2012 2016 2008 2004 2016 2014 2010 2019 2005 2016 2000 2010 2018 2010 2016 2010 2018 2015 2016 2012 2012 2007 2000 2014 2001 2017 2011 2014 2010 2017 2014 2017 2018

119 552 440 343 343 179 469 469 222 536 266 343 501 469 380 39 512 453 726 251 119 694 552 343 1059 139 694 39 325 1009 348 694 337 976 672 694 445 422 179 39 568 41 833 384 552 337 812 582 802 963

24 5 19 19 11 16 11 5 6 21 5 15 18 19 24 30 15 13 21 18 7 16 4 1 15 11 24 12 2 11 5 5 24 5 5 15 10 21 15 18 18 5 24 15 7 37 18 7 16 15

18 11.356 15.994 14.058 14.058 17.962 18 18 17.998 8.319 18 14.058 18 18 18 17.763 18 16.178 15.462 11.848 12.488 18 10.740 7.910 13.600 18 18 17.996 18 9.832 3.349 18 18 18 18 18 12.025 10.751 17.695 17.781 18 18 18 7.656 10.783 18 18 18 18 9.019

5.923 6.041 6.076 6.206 6.219 6.224 6.257 6.268 6.274 6.299 6.301 6.301 6.305 6.316 6.352 6.402 6.458 6.466 6.477 6.485 6.495 6.527 6.534 6.535 6.536 6.547 6.572 6.580 6.585 6.602 6.605 6.614 6.630 6.632 6.646 6.653 6.654 6.668 6.684 6.691 6.691 6.693 6.698 6.705 6.713 6.714 6.717 6.723 6.744 6.751

10.166 9.411 10.076 9.956 9.969 10.462 10.500 10.511 10.517 9.183 10.543 10.050 10.548 10.558 10.595 10.617 10.701 10.488 10.409 9.928 10.029 10.770 9.812 9.348 10.224 10.790 10.815 10.822 10.828 9.738 8.435 10.857 10.873 10.874 10.889 10.895 10.122 9.946 10.890 10.908 10.934 10.935 10.940 9.472 9.997 10.957 10.960 10.965 10.986 9.754

1309 891 882 1307 1393 1261 1127 1258 818 922 921 1255 928 1034 1109 1384 1018 1034 945 1267 911 895 1046 1327 1208 1306 914 924 1143 1365 1279 866 1226 1198 1047 887 945 962 1241 1235 948 1245 1235 1046 901 1035 799 1096 1058 1168

VJ1 JR24 FS35 TE55 TE55 BB VC26 VC26 CM74 FW32 YN2 TE55 PG10 VC26 DS SZ162 TT5 PB20 HF19 ST VJ1 CO246 JR24 TE55 DJ1 VN5 CO246 SZ162 JR34 SD2 UJ CO246 RF12 KP1 YA1 CO246 HK31 BB14 BB SZ162 QH33 BA16 FR ED12 JR24 RF12 BU UV34 BB7 GE2

v0

vn1 , vn 2

vega1T

vdsm1

vega 2T

vf

Earth

NEA

EGA1

DSM1

EGA2

MBA

0

t (2), t (3)

t (1)

Lam 5

1 t (4), Lam

E

,

E

2 , vega1T OP

3 t (5)

4 t (6)

Lam

6

t (7)

Lam 7

Fig. 6. Patched conic method for the ENEEM trajectory.

instantaneous velocity increment during the gravity assist. The advantages of this model lie in the simplicity of use and the efficiency of computation. It has been used in many research works (see (Sims, 1997), (Vasile and Pascale, 2006), (Ceriotti and Vasile, 2010), (Chen et al., 2014) and (Yang et al., 2015)). Based on the assumptions above, the following formulas can be obtained

the gravity-assist planet. When the spacecraft is in the planet’s sphere of influence, only the central gravitational force of the planet is considered. Considering that both the flight time and space when the spacecraft is in the heliocentric elliptical orbit are much larger than those when it is in the sphere of influence of the planet, the instantaneous positions of the spacecraft before and after the gravity assist are assumed to be equal to the position of the gravity-assist planet. The spacecraft obtains an 5

P. Sun et al.

Planetary and Space Science 182 (2020) 104851

v0

vn1 , vn 2

vega1T

vmga1T

vf

Earth

NEA

EGA1

MGA1

MBA

0

1

t (1)

2 t (4)

t (2), t (3) Lam

3 t (6)

t (5)

Lam

Lam 6

Fig. 7. Patched conic method for the ENEMM trajectory.

v ∞ ¼ v  vp

¼

(9)

v  ∞ v ∞

 vp ;j ¼ k  i  vp 



μp

μp þ rp v∞

(4)

(10)

2

where μp is the gravitational constant of the gravity-assist planet. rp ¼ RP þ hga is the hyperbolic periapsis radius about the gravity-assist planet, in which RP is the radius of the gravity-assist planet, hga is the minimum attitude above the planet surface during gravity assist, and hga  hmin ¼ 300 km (the admissible minimum altitude in this study). When hga equals hmin, the maximum turn angle δPmax is obtained. In summary, when rp and ψ are specified, the impulse Δvga gained during the gravity assist can be computed by

where v is the heliocentric velocity of the spacecraft immediately before gravity assist, vþ is the heliocentric velocity of the spacecraft immediately after gravity assist, and vp is the velocity of the gravity-assist planet. A reference frame P-ξηζ of gravity assist impulsive model is defined in Fig. 2. The origin of the frame is the center of mass of the gravity-assist planet. The ξ axis is aligned to the hyperbolic excess velocity of the spacecraft before gravity assist. The i, j and k are used to represent the unit vectors along the ξ, η and ζ axis, respectively. There exists the following relationship v ∞ i¼ v  ; k ∞

v ∞  vþ ∞ ¼ v∞ 2 cos δ

 δ ¼ 2 arcsin

(3)

vþ ∞ ¼ vþ  vp

(8)

where v∞ is the magnitude of the hyperbolic excess velocity of the spacecraft. ψ is the direction angle between the projection of v-∞ in the Pηξ plane and the ζ axis, and ψ 2[0,2π). δ is the turn angle between the hyperbolic excess velocity before and after gravity assist, and is determined by (Sims, 1997)

Lam

5

4

kv ∞ k ¼ kvþ ∞ k ¼ v∞

△vga ¼ vþ ∞  v ∞ ¼ vþ  v ¼ v∞ ððcos δ  1Þi þ sin δ sin ψ j þ sin δ cos ψ kÞ

(5)

Also, the heliocentric velocity v of the spacecraft immediately after gravity assist can be obtained

Case 1. It is assumed that no maneuver impulses are performed by the spacecraft during the gravity assist, and no impulses are executed after gravity assist. One has (see (Jiang et al., 2012) and (Yang et al., 2018)) v ∞ ¼ v∞ i

(6)

vþ ∞ ¼ v∞ ðcos δi þ sin δ sin ψ j þ sin δ cos ψ kÞ

(7)

vþ ¼ v þ △vga

(12)

Case 2. It is assumed that no impulses are performed by spacecraft during the gravity assist. However, there can be an extra maneuver impulse executed instantaneously after the gravity assist. Assume that the heliocentric velocity of spacecraft immediately

v0

vn1 , vn 2

vega1T

vdsm1

vega 2T

vmga1T

vf

Earth

NEA

EGA1

DSM1

EGA2

MGA1

MBA

0 t (1)

1 t (4), Lam

t (2), t (3) Lam

2 , v E ega1T OP

E,

6

3

5

4 t (6)

t (5)

t (7)

Lam

Lam

7

8

t (8) Lam

9

Fig. 8. Patched conic method for the ENEEMM trajectory.

v0

vn1 , vn 2

vega1T

vmga1T

vdsm1

vmga 2T

vf

Earth

NEA

EGA1

MGA1

DSM1

MGA2

MBA

0 t (1)

1

t (2), t (3)

t (4)

Lam

Lam 6

2 t (5), Lam

7

M

,

M

3 , vmga1T OP

5

4 t (6)

t (8)

t (7)

Lam

Lam

8

9

Fig. 9. Patched conic method for the ENEMMM trajectory.

v0

vn1 , vn 2

vega1T

vdsm1

vega 2T

vmga1T

vdsm1

vmga 2T

vf

Earth

NEA

EGA1

DSM1

EGA2

MGA1

DSM2

MGA2

MBA

0 t (1)

t (2), t (3)

Lam 8

(11)

þ

1 t (4), Lam

E

,

9

E

2 , vega1T OP

3 t (6)

t (5)

Lam 10

4 t (7), Lam

M

,

M

5 , vmga1T OP

11

Fig. 10. Patched conic method for ENEEMMM. 6

7

6 t (9)

t (8)

Lam

t (10)

Lam 12

P. Sun et al.

Planetary and Space Science 182 (2020) 104851

Table 3 Optimized solutions of the 50 target asteroid groups. ENEEM

ENEMM 2

N_50

J/km/s

C3/km /s

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

7.506 6.868 7.404 5.184 4.944 9.649 8.093 5.972 9.779 8.782 7.751 7.252 8.39 5.587 8.181 8.076 7.878 6.541 8.422 7.352 8.959 7.885 9.919 8.104 6.314 5.7 6.748 8.337 8.409 8.24 9.041 6.447 5.885 7.974 7.747 6.48 9.843 7.989 8.551 6.466 8.71 9.103 8.563 7.841 6.972 8.096 8.891 8.1 10.004 9.173

12.459 10.763 16.105 8.641 8.549 14.548 18 8.431 11.066 8.345 17.991 7.458 18 17.998 17.255 17.839 18 17.997 15.175 9.981 17.997 18 18 8.668 12.952 18 18 17.998 13.521 3.84 2.513 18 18 14.506 16.88 18 7.857 10.746 8.04 17.766 18 18 10.81 14.653 10.75 18 18 18 17.998 18

2

2

2

Tof/day

J/km/s

C3/km /s

Tof/day

1353 1499 1510 1738 1815 1293 1546 1416 1476 1405 1408 1584 1476 2350 1552 1824 1876 1813 1430 1775 1391 1582 1487 1451 1433 1872 1632 1394 1580 1665 1955 1912 1447 1810 1466 1605 1377 1471 1611 1780 1568 1595 1576 1986 1422 1479 1346 1623 1607 1325

7.822 5.521 7.571 8.304 8.101 9.541 8.991 5.289 7.246 7.535 6.08 9.284 10.568 7.673 8.707 9.431 6.763 7.094 7.791 7.32 8.733 6.011 6.236 4.826 5.829 10.307 6.743 7.877 8.395 8.106 8.231 7.544 8.224 9.288 7.431 6.804 7.582 7.993 9.074 7.102 8.914 9.676 10.724 9.376 6.808 7.729 8.613 8.304 11.585 8.115

18 11.476 13.588 4.94 3.996 7.725 18 17.342 16.718 15.094 18 7.877 18 17.485 18 18 18 12.606 6.97 8.606 18 18 11.218 7.748 14.543 18 18 17.793 7.139 9.594 6.775 18 18 18 18 18 12.883 7.407 16.815 18 14.232 15.707 12.319 18 11.194 18 18 18 18 3.626

1412 1255 1775 1658 1691 1634 1362 1399 1859 2046 1276 1772 1085 1433 1244 1528 1560 1853 2072 1612 1664 1584 2085 1377 1933 1666 1643 1410 1469 1613 1510 1419 1702 1833 1141 1585 1604 2076 1454 1634 1113 1449 1742 1529 2049 1716 1039 2024 1393 1278

(17)

(18)

To sum up, with the specified heliocentric velocity of spacecraft immediately before and after gravity assist, the optimal turn angle δ*, and the optimal direction angle ψ * can be computed by Eq. (16) and Eq. (17). Thus the hyperbolic periapsis radius rp and the maneuver impulse ΔvT performed by the spacecraft are given. The mapping relationship is as follows 

   △ vga ; △vT ; δ* ; ψ * ; rp ¼ GA T1 v ; vþ ; μp

(19)

2.3. Preliminary gravity-assist sequence analysis Tisserand graph is an efficient and generic technique used in MBA exploration missions (see (Strange and Longuski, 2002) (Chen et al., 2014), and (Yang et al., 2015)), which permits the quick identification of a ballistic path of gravity-assist planets. In this paper, a plot of Ra vs Rp (Ra-Rp plot) is employed. The fundamental principle is that both the apoapsis Ra and the periapsis Rp are calculated by the state vectors of spacecraft. Each pair of (Ra, Rp) determines an orbit in the ecliptic. The farthest point to the upper right on the v∞ contour corresponds to positive alignment of the spacecraft’s hyperbolic excess velocity vector vþ∞ after gravity assist with the planet’s velocity vp (αP ¼ 0 deg), which represents the highest heliocentric energy possible for encounter with the planet. Rotation of the v∞ vector away from this alignment corresponds to moving from the upper right to the bottom left on the v∞ contour toward negative alignment with the planet’s velocity vector (αP ¼ 180 deg), which corresponds to the lowest heliocentric energy. Usually, the point of αP ¼ 0 deg is taken as the first tick mark node, and the bottom left point of αP ¼ 180 deg is not marked. There are other n points with αP ¼ k⋅δPmax (k ¼ 1,2, …,n; k 2 Z*, Z* is the set of positive integers, and n⋅δPmax < 180 deg) marked on the v∞ contour. From a given tick mark on a contour, the spacecraft’s orbit may move a maximum distance indicated by the next tick mark (either up or down that contour) before violating the 300-km altitude constraint. A Ra-Rp plot combining several v∞ contours for NEA 2014 JR24, Earth, Mars and MBA Urania is shown in Fig. 3a. It is assumed that all of the planets are in circular and coplanar orbits. It is also assumed that the spacecraft is in the same plane as the planets (see (Strange and Longuski, 2002)). The contours in Fig. 3a start next to each planet’s name at 1 km/s and increase in steps of 2 km/s toward the upper left of the plot. All tick marks are spaced for the 300-km constraint. Fig. 3a shows that if there is an appropriate phasing between the Earth and 2014 JR24, it is not that difficult to transfer from Earth to 2014 JR24 and accomplish the NEA sample return mission. This is because an intersection point between the 1 km/s contour of Earth and 1 km/s contour of 2014 JR24 exists, indicating a potential low-energy transfer orbit. In the research work of MBA accessibility analysis, dual Mars gravity assists (MMGA) (see (Chen et al., 2014) and (Yang et al., 2015)) are effective to lower the launch energy from Earth after the subtask of NEA sample return and to reduce the rendezvous delta-V to match the orbit of target MBA. The Earth gravity assist is necessary for reducing the total impulse and to return the sample collected from NEA according to the research work of (Yang et al., 2015). The results in Fig. 3a and b are consistent with those aforementioned. Assuming that no Mars gravity assist (MGA) is considered when the spacecraft departs from Earth for Urania, at least a 6 km/s Earth launch velocity increment is required, and at least a 6 km/s arrival velocity

(13)

(14)

where v1, v2 and v3 denote the components of the velocity vector vþ-valong the directions of the orthonormal triad {i, j, k}. Then the maneuver impulse ΔvT after gravity assist is determined △vT ¼ ðv1  cos δ þ 1Þi þ ðv2  sin δ sin ψ Þj þ ðv3  sin δ cos ψ Þk

δ* ¼ minðδ; γÞ ! 1 þ v1 γ ¼ arccos qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1 þ v1 Þ2 þ v2 2 þ v3 2

where Δvga can be obtained from Case 1. Project the velocity vector vþ-vto the frame of P-ξηζ, and the following equation can be obtained vþ  v ¼ v1 i þ v2 j þ v3 k

(16)

The expression of γ is given by (see (Jiang et al., 2012))

before and after the gravity assist is v and vþ, respectively. The instantaneous maneuver impulse performed by spacecraft after gravity assist is ΔvT. It can be defined that vþ  v ¼ △vga þ △vT

ψ * ¼ arctan 2ðv2 ; v3 Þ

(15)

To reduce the total fuel consumption, δ and ψ are optimized for the minimum magnitude of the maneuver impulse ΔvT. The optimal results are as follows (see (Jiang et al., 2012)) 7

P. Sun et al.

Planetary and Space Science 182 (2020) 104851

Table 4 Optimized solutions of the 50 target asteroid groups. N_50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

ENEEMM

ENEMMM 2

J/km/s

C3/km /s

7.408 6.275 7.872 6.22 6.939 7.576 7.601 6.972 7.123 8.717 7.746 6.725 9.531 7.152 8.386 8.309 6.363 6.06 8.382 6.319 7.187 7.303 7.958 6.583 5.505 8.089 7.86 7.763 6.046 6.61 6.841 5.985 6.405 8.122 7.786 6.708 8.352 7.937 7.108 7.566 8.284 6.566 9.594 6.374 7.338 7.85 7.98 8.36 7.564 6.147

0.849 10.72 15.879 7.398 8.113 16.825 18 16.005 17.834 8.334 17.938 8.385 18 18 18 16.188 18 17.825 15.375 15.042 18 18 10.716 7.654 11.73 17.998 18 17.794 16.259 14.166 7.617 18 18 13.924 18 18 12.845 10.782 16.324 17.069 18 18 10.341 8.312 10.706 17.886 17.998 18 18 18

2

ENEEMMM 2

Tof/day

J/km/s

C3/km /s

2080 2143 1902 2148 2085 2714 2305 2332 2304 2241 1813 2188 2744 2287 2024 2047 2453 2533 2271 2048 2172 1991 2326 2607 2292 2110 2547 2741 2159 2160 1785 1878 2781 2147 2077 1809 2744 2297 2708 2416 2630 1824 2283 2451 2222 2748 2314 2420 2521 2417

7.22 4.339 6.74 4.691 7.989 5.898 5.857 4.499 7.19 6.39 5.707 4.561 6.667 5.212 7.149 6.346 5.13 6.711 5.992 4.896 5.448 5.398 4.889 4.742 4.539 6.522 7.036 5.969 4.647 8.144 5.489 5.49 5.98 6.133 7.434 5.518 8.135 6.881 5.533 5.371 8.295 7.63 6.188 5.313 3.495 6.484 6.409 5.661 6.091 6.024

17.979 11.096 13.68 7.28 4.184 16.685 18 17.998 16.765 11.269 18 7.299 18 17.994 17.999 17.997 18 16.248 13.686 11.922 17.835 18 10.774 7.358 14.153 18 18 18 14.522 9.571 9.541 18 17.744 18 15.823 18 13.341 12.733 16.866 17.851 10.332 13.204 18 7.683 11.594 17.794 18 18 18 18

2

Tof/day

J/km/s

C3/km2/s2

Tof/day

3069 2425 3284 2178 3249 2894 2414 2299 3534 2454 2307 2292 2885 2455 2355 2719 2589 2330 2295 2574 2301 3205 2375 2805 2126 3451 2879 3323 2771 3132 3288 2236 2843 3237 2958 3190 2904 2438 2838 3172 2464 3235 2435 2601 2289 2913 2789 2487 2698 2551

7.162 5.992 5.951 6.136 6.203 7.167 6.968 4.412 6.371 6.967 6.686 5.968 6.488 6.271 8.489 6.581 7.353 4.412 6.951 5.817 5.682 6.135 6.892 5.448 4.789 7.42 6.425 6.197 5.34 7.462 6.519 5.163 5.954 5.688 7.57 5.409 9.281 6.264 6.822 6.785 6.256 6.31 8.192 7.08 5.806 6.76 7.003 6.938 7.83 7.335

12.9 17.97 15.89 7.566 7.357 6.847 18 11.514 14.824 8.368 9.618 7.419 18 18 17.94 17.999 18 16.2 15.399 9.448 6.363 17.97 16.584 8.378 12.097 16.884 18 12.759 13.002 13.529 2.443 18 10.494 13.555 15.993 18 11.95 12.578 6.085 18 18 18 10.989 7.657 10.716 18 18 17.864 18 18

2618 2875 3338 3561 3506 3248 2896 2961 3579 3720 3255 3656 2903 2947 2741 2760 2937 3726 3733 3072 3662 3206 2843 2910 3134 3388 2845 2863 2867 3673 3345 2885 3207 3329 2880 3461 3038 3798 3186 3192 2937 3584 2844 3978 3919 3073 2766 2841 3832 3926

Table 5 Comparison results of the postlaunch delta-V with the work of (Chen et al., 2014). MBA_name

Flora Iris Vesta Victoria Parthenope Astraea

N_50

24 26 9 23 38 33

ENEM

DT

△v* p /km/s

△vp /km/s

4.652 3.568 4.148 5.505 4.084 4.019

4.223 3.565 4.095 4.538 4.16 4.127

Evp

0.102 0.001 0.013 0.213 0.018 0.026

N_50

24 5 9 23 10 27

ENEMM

MGA

△v* p /km/s

△vp /km/s

2.765 4.711 3.848 3.934 4.94 4.747

2.795 4.679 4.752 4.862 4.071 3.854

Evp

0.011 0.007 0.190 0.191 0.213 0.232

N_50

24 7 9 23 19 27

ENEMMM

MMGA

△v* p /km/s

△vp /km/s

3.22 4.376 3.768 3.562 3.909 4.166

1.58 3.872 3.372 3.152 3.337 4.169

Evp

1.038 0.130 0.117 0.130 0.171 0.001

km/s. Fig. 3b shows that about one Earth flyby and two Mars flybys are needed in designing the multi-target asteroid exploration trajectory starting from Earth, which includes subtasks of the 2014 JR24 sample return and Urania rendezvous. It can be known that before calculating the minimum number of flyby bodies needed to reach Urania, a rough

relative to Urania is required, both of which cause unacceptable fuel consumption. When both EGA and MGA are considered, the Urania rendezvous delta-V can be reduced to 2.5 km/s. Also, the Earth v∞ contour of 3.5 km/s has an intersection point with the Mars v∞ contour of 3 km/s. After the subtask of NEA sample return, which is accomplished by EGA1, the Earth departure hyperbolic excess velocity can be reduced to 2

8

P. Sun et al.

Planetary and Space Science 182 (2020) 104851

Table 6 Comparison results of the postlaunch delta-V with the work of (Chen et al., 2014). MBA_name

N_50

ENEEM

EGA

△v p /km/s

△vp /km/s

6.895 5.871 6.359 7.485 6.942 5.219

5.685 5.852 6.107 6.26 5.752 5.865

*

Flora Iris Vesta Victoria Parthenope Astraea

24 30 9 23 19 15

Evp

0.2128 0.0033 0.0413 0.1957 0.2069 0.1102

ENEEMM

EMGA

N_50

△v p /km/s

△vp /km/s

24 5 9 23 19 33

5.135 5.73 4.768 6.915 6.914 5.056

3.985 5.806 5.118 5.822 5.063 5.456

*

Evp

0.2886 0.0131 0.0684 0.1877 0.3655 0.0733

Table 7 Comparison results of the mission duration from the Earth (EGA1, in this paper) to the MBA. comparison value

TEM/day

MBA_name

Astraea Harmonia Victoria Athamantis Massalia

Chen et al.

This paper

Etem

type

variable

value/day

type

N_50

variable

value/day

DT MGA MMGA EGA EMGA

TOF

470.68 1006.49 1822.67 890.811 1264.24

ENEM ENEMM ENEMMM ENEEM ENEEMM

33 29 23 37 36

tEM

473 1015 1820 881 1261

a) The minimum, median, and maximum C3 of the 50

b) The minimum, median, and maximum J of the 50 asteroid

asteroids groups for six types of trajectories

groups for six types of trajectories

0.0049 0.0084 0.0015 0.011 0.0026

Fig. 11. a) The minimum, median, and maximum C3, and b) J of the 50 asteroids groups for six types of trajectories.

to Mars along Earth’s v∞ ¼ 3.5 km/s contour is less than the tick mark spacing for one flyby, and the distance to Urania along Mar’s v∞ ¼ 3 km/s contour is greater than the tick mark spacing for one flyby, but less than the spacing of two tick marks. In particular, it is much more straightforward and intuitive to give a quick reference of the number of flyby bodies needed to reach Urania, i.e., one Earth flyby and two Mars flybys, which is consistent with the result of the work of (Yang et al., 2015). The result can be directly obtained through the position distribution of the modified tick marks in Fig. 4b, instead of making a rough estimation according to the tick mark spacing in the Ra-Rp plot with conventional tick marks. The Ra-Rp plot with modified tick marks provides a quick and intuitive way to determine the minimum number of gravity-assist bodies. Furthermore, the analysis results found here can serve as a reference for the gravity-assist sequence analysis of other multi-target asteroid exploration missions.

estimate must be made according to the tick mark spacing along Earth’s v∞ ¼ 3.5 km/s contour and Mar’s v∞ ¼ 3 km/s contour, which is not intuitive and exact from Fig. 3b. In order to indicate exactly and accurately the number of flyby bodies in sequence analysis, a Ra-Rp plot with modified tick marks are introduced in Fig. 4a. Thereinto, the intersection point of Earth’s v∞ ¼ 3.5 km/s contour and Mar’s v∞ ¼ 3 km/s contour is labeled as the initial tick mark node. The angle value αP0 (for this case, αE0 for Earth and αM0 for Mars) of the current intersection point corresponding to the v∞ contours of two gravity-assist planets is computed. On the corresponding contour lines, taking the initial node αP0 as the starting point, the other n1 nodes with αP decreasing and increasing in steps of δPmax, i.e., αP ¼ αP0þk1⋅δPmax (k12Z, and Z is the set of integers; αP > 0 deg and αP < 180 deg; the total number of k1 eligible for these constraints is n1) are chosen as new tick marks. Based on the Ra-Rp plot with modified tick marks, a renewed gravityassist sequence analysis result is shown in Fig. 4b. Thereinto, the distance

9

P. Sun et al.

Planetary and Space Science 182 (2020) 104851

a) The minimum, median, and maximum Tof of the 50

b) The minimum, median, and maximum vf of the 50 asteroid

asteroids groups for six types of trajectories

groups for six types of trajectories

Fig. 12. a) The minimum, median, and maximum Tof, and b) Δvf of the 50 asteroids groups for six types of trajectories.

Table 8 The accessibility list for the ENEMMM trajectory.

Table 9 Three optimal example missions for the ENEMMM trajectory.

No.

N_50

NEA_num

MBA_num

J/km/s

Tof/day

N_50

45

2

8

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

45 2 8 25 12 29 4 24 23 20 17 14 44 40 22 21 31 32 36 39 48 11 7 6 28 33 19 50 49 34 43

552 552 469 1059 343 325 343 343 552 251 512 469 384 39 694 119 348 694 694 179 582 266 469 179 39 337 726 963 802 976 833

7 5 5 15 15 2 19 1 4 18 15 19 15 18 16 7 5 5 15 15 7 5 11 16 12 24 21 15 16 5 24

3.49517 4.338651 4.499152 4.538501 4.560669 4.646848 4.690841 4.741521 4.888554 4.896414 5.130327 5.212002 5.313218 5.371291 5.398322 5.448237 5.48937 5.490384 5.517572 5.532662 5.661066 5.707384 5.857119 5.898057 5.969188 5.980209 5.9918 6.023641 6.090947 6.133154 6.187509

2289.375 2425.183 2299.074 2126.284 2291.887 2770.734 2178.107 2804.557 2374.34 2574.225 2588.865 2455.005 2600.902 3172.502 3204.328 2300.911 3287.134 2235.804 3190.097 2838.531 2487.196 2307.28 2414.096 2893.246 3323.811 2842.683 2295.112 2550.299 2698.489 3237.135 2434.286

{NEA_num,MBA_num}

{552,7}

{552,5}

{469,5}

t(1) C3 t(2) △vn1 t(3) △vn2 t(4) △vega1T rp1 t(5) △vmga1T rp2 t(6) △vdsm1 t(7) △vmga2T rp3 t(8) △vf

{2025-05-13} 11.0135 {2026-01-15} 0.1733 {2026-5-10} 1.2251 {2026-11-23} 9.9E-11 12383.93 {2027-9-9} 1.66E-11 3734.55 {2029-2-26} 0.2265 {2030-5-26} 9.02E-10 3837.47 {2031-8-20} 1.8703

{2025-05-13} 11.09551 {2026-02-01} 0.147383 {2026-06-05} 1.504497 {2026-12-16} 1.69E-06 20776.8 {2027-12-03} 8.89E-05 4867.141 {2029-07-17} 0.434335 {2030-08-23} 0.984093 3695 {2032-01-02} 1.268252

{2025-11-01} 17.99823 {2026-01-12} 0.563228 {2026-03-03} 0.916151 {2026-12-02} 8.47E-09 23232.01 {2028-03-22} 0.000645 12802.82 {2029-06-23} 0.987317 {2030-10-02} 0.928766 3695 {2032-02-17} 1.103045

ENEEMMM: 18 feasible examples.

3. Optimal delta-V searching method To design the multi-target asteroid exploration trajectory, a set of following assumptions are made: 1) Launch is assumed to take place between 2022/01/01 and 2028/01/ 01. t0s ¼ 2022/01/01, and t0f ¼ 2028/01/01. 2) The maximum launch energy is C3max ¼ 18 km2/s2. 3) The initial Earth launch velocity increment provided by spacecraft is

2.4. Problem statement To evaluate the accessibility of NEAs and MBAs in a multi-target mission including both NEA sample return (by EGA1) and MBA rendezvous, a patched conic or simple two-body orbit is used, and thus multiimpulse trajectories are calculated with launch dates between 2022 and 2028. Considering the huge number of candidate asteroids, the long period of a multi-target asteroid exploration mission and the difficulty in a multi-segment orbit optimization problem, the technique of gravity assist is introduced, and the PSO algorithm is employed to search for the global minimum velocity increments of the multi-impulse trajectories.

pffiffiffiffiffiffi 0; ð△v00pffiffiffiffiffiffi 18 km=sÞ pffiffiffiffiffiffi . Δv00 ¼ ||Vd0△v00  18 ð△v00 > 18 km=sÞ VE0||, where Vd0 is the spacecraft’s initial velocity vector on the Earth to the NEA transfer leg, and VE0 is the Earth velocity vector at departure time. 

△v0 ¼

4) The launch energy provided by a launch vehicle is assumed to be C3 ¼

10

P. Sun et al.

Planetary and Space Science 182 (2020) 104851

Table 10 The accessibility list for the ENEEMMM trajectory. No.

N_50

NEA_num

MBA_num

J/km/s

Tof/day

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

18 8 25 32 29 36 24 21 34 45 20 3 33 12 2 22 4 28

453 469 1059 694 325 694 343 119 976 552 251 440 337 343 552 694 343 39

13 5 15 5 2 15 1 7 5 7 18 19 24 15 5 16 19 12

4.411766 4.412489 4.789292 5.163008 5.339875 5.408616 5.447624 5.682222 5.688408 5.806042 5.817286 5.950987 5.953508 5.967523 5.992351 6.135148 6.136274 6.196678

3725.426 2960.359 3133.77 2885.168 2866.744 3461.081 2909.633 3662.165 3328.41 3919.122 3071.862 3338.173 3206.456 3655.672 2874.942 3206.867 3561.032 2863.175

Fig. 13. Ballistic transfer trajectory for the ENEMMM trajectory.

quickly as possible. Different from NEAs, the orbits of MBAs are distant from Earth (see (Yang et al., 2015)). Thus, for MBA rendezvous, single and multiple gravity assists, such as EGA, MGA, MMGA and Earth-Mars gravity assist (EMGA), are studied to lower the mission cost (see (Chen et al., 2014)). To further investigate the combined effect of EGA and MMGA, an extra type of gravity-assist trajectory, namely Earth-Mars-Mars gravity assist (EMMGA), is considered. For all above considerations, the following six types of gravity-assist trajectories are studied in this paper.

Table 11 Three optimal example missions for the ENEEMMM trajectory. N_50

18

8

25

{NEA_num,MBA_num}

{453,13}

{469,5}

{1059,15}

t(1) C3 t(2) △vn1 t(3) △vn2 t(4) △vega1T rp1 t(5) △vdsm1 t(6) △vega2T rp2 t(7) △vmga1T rp3 t(8) △vdsm2 t(9) △vmga2T rp4 t(10) △vf

{2025-02-11} 16.20009 {2025-11-01} 0.077334 {2026-01-09} 0.818981 {2026-08-14} 0.002172 261854.2 {2028-01-28} 1.053868 {2028-12-11} 4.09E-07 22631.58 {2029-10-30} 1.38E-05 5388.763 {2032-01-01} 0.000107 {2033-08-04} 0.00022 3769.202 {2035-04-26} 2.45907

{2023-11-24} 11.5139 {2024-10-24} 1.553151 {2025-03-02} 0.128646 {2025-11-12} 0.136053 31134.68 {2026-06-19} 1.74E-06 {2026-11-12} 4.77E-08 6972.314 {2027-12-02} 1.21E-06 4641.295 {2029-06-17} 0.41135 {2030-08-22} 0.905365 3695 {2032-01-02} 1.277922

{2024-03-25} 12.09707 {2024-10-30} 0.433864 {2024-12-20} 0.654873 {2025-08-30} 7.09E-05 30474.77 {2026-06-24} 0.401239 {2027-02-15} 0.000129 14789.54 {2028-04-24} 0.000197 4202.02 {2029-12-21} 0.029782 {2031-03-20} 0.528936 3695 {2032-10-23} 2.740201

Type 1: Earth-NEA(sampling)-EGA1(sample return)-MBA (ENEM); Type 2: Earth-NEA-EGA1-DSM1(the first deep space maneuver(DSM))-EGA2(the second EGA)-MBA (ENEEM); Type 3: Earth-NEA-EGA1-MGA1(the first MGA)-MBA (ENEMM); Type 4: Earth-NEA-EGA1-DSM1-EGA2-MGA1-MBA (ENEEMM); Type 5: Earth-NEA-EGA1-MGA1-DSM1-MGA2(the second MGA)MBA (ENEMMM); Type 6: Earth-NEA-EGA1-DSM1-EGA2-MGA1-DSM2(the second DSM)-MGA2-MBA (ENEEMMM). The PSO algorithm is applied in the optimization of the six types of trajectories. All the default parameter settings for the PSO algorithm are the same as those in (Chen et al., 2014) except that both the values of cognitive and social weights are updated within the range of [0.5, 2.5]. To make a more concise and efficient introduction to the design and optimization process for the six types of trajectories, some extra common definitions and general attention are given first:

(

5) 6) 7) 8) 9)

pffiffiffiffiffiffi △v200 ; ð△v00  18 km=sÞ . pffiffiffiffiffiffi 18; ð△v00 > 18 km=sÞ The NEA sampling duration tsam is assumed to be in the range of [50,150] days. The attitude above the surface of all gravity-assist planets is assumed to be greater than 300 km, i.e., hga300 km. The total mission duration is assumed to be less than 4000 days, i.e., Tof4000 days. The initial mass of the spacecraft is assumed to be m0 ¼ 2000 kg, and the specific impulse of the chemical propulsion is set to be Isp ¼ 350 s. The orbits of Earth and Mars are computed with the DE421 ephemeris from JPL. The transfer trajectory of the spacecraft is propagated with the two-body dynamical model.

1) t(1) is the initial Earth departure time. t(2) is the NEA arrival time. t(3) is the NEA departure time after completing the sampling operation. t(4) is the time of the EGA1. The specific definition of the variable t(num) (num>4, num 2 Z*) should refer to that given for each type of gravity-assist trajectory in the following. 2) Assuming that the solving process for one type of gravity-assist trajectory can be decomposed into ct_seg (ct_seg>0, ct_seg 2 Z*) segments, ranging from Segment 0 to Segment ct_seg-1. tofi (i2[0,ct_seg-1], i 2 Z*) denotes the flight time of Segment i. tofi equal 0 if Segment i denotes the solving process of Earth or Mars flyby. This is because the spacecraft is assumed to be perturbed instantaneously by the gravityassist body during the flyby. 3) Δvn1 is the magnitude of NEA arrival hyperbolic excess velocity vector. Δvn2 is the magnitude of NEA departure hyperbolic excess velocity vector after finishing the sampling operation. Δvf is the magnitude of terminal brake impulse required to reach the orbit of the target MBA. Δvega1T andΔvega2T are the maneuver impulses

To accomplish the combined mission at a lower cost and with a shorter duration, only the direct transfer trajectories from Earth to NEAs are considered. EGA1 is used to return the sample collected from NEAs as

11

P. Sun et al.

Planetary and Space Science 182 (2020) 104851

a)

X-Y plot

b)

X-Z plot

Fig. 14. a) X-Y, and b) X-Z plot for the ENEMMM trajectory.

a)

The semi-major axis a

b)

The eccentricity e

Fig. 15. a) The semi-major axis, and b) eccentricity for the optimized ENEMMM trajectory over time.

process can be decomposed into ct_seg ¼ 5 segments. Segment 0, 1 and 2 are computed with a Lambert solver (Lam in Fig. 5, for short). Segment 3 denotes the NEA sampling process. Segment 4 denotes the solving process of Earth flyby. The state vector of Earth at t(1), [RE0,VE0]T, and the state vector of NEA at t(2), [Rn1,Vn1]T, are computed by ephemerides. With a lambert solver, the spacecraft state vectors corresponding to the initial time and the final time of Segment i (i2{0,1,2}) are determined by the following mapping relationships

executed instantaneously after the EGA1 and EGA2, respectively. Δvdsm1 andΔvdsm2 are the magnitudes of the impulse for DSM1 and DSM2, respectively. Δvmga1T andΔvmga2T are the magnitudes of the impulse executed instantaneously after MGA1 and MGA2, respectively.

3.1. ENEM Considering the huge number of candidate NEAs and MBAs to be optimized, Fig. 5 firstly shows the most straightforward type of trajectory, where t(5) is the final MBA rendezvous time. The mission solving 12

ðVd0 ; Va0 Þ ¼ Lam solðRE0 ; Rn1 ; tof0 Þ

(20)

ðVd1 ; Va1 Þ ¼ Lam solðRn2 ; RE1 ; tof1 Þ

(21)

P. Sun et al.

Planetary and Space Science 182 (2020) 104851

a)

The inclination i

b)

The RAAN

Fig. 16. a) The inclination, and b) right ascension of ascending node for the optimized ENEMMM trajectory over time.

a)

The argument of periapsis

b)

The true anomaly f

Fig. 17. a) The argument of periapsis, and b) true anomaly for the optimized ENEMMM trajectory over time.

  ðVd2 ; Va2 Þ ¼ Lam sol RE1 ; Rf ; tof2

where μE is the gravitational constant of Earth, Δvega1 is the impulse provided by EGA1, δ*E and ψ *E are the optimal gravity-assist angles, and rp is the flyby radius of the EGA1. Thus, the following equations are satisfied

(22)

where Vdi, i2{0,1,2}, is the spacecraft’s velocity vector corresponding to the initial time of Segment i, while Vai, i2{0,1,2}, is the spacecraft’s velocity vector corresponding to the final time of Segment i. [Rn2,Vn2]T is the NEA state vector at t(3), [RE1,VE1]T is the Earth state vector at t(4), and [Rf,Vf]T is the MBA state vector at t(5). Based on Eq. (19), the following mapping relationship of Segment 4 can be deduced 



△ vega1 ; △vega1T ; δE ; ψ E ; rp ¼ GA T1ðVa1 ; Vd2 ; μE Þ *

*

(23) 13

△vn1 ¼ kVn1  Va0 k

(24)

△vn2 ¼ kVd1  Vn2 k

(25)

  △vf ¼ Vf  Va2 

(26)

P. Sun et al.

Planetary and Space Science 182 (2020) 104851

a)

The distance

b) The heliocentric velocity

Fig. 18. a) The distance, and b) heliocentric velocity over time. Case 2: ENEEMMM

a)

3D transfer trajectory for the ENEEMMM trajectory

b) X-Y plot

Fig. 19. a) 3D, and b) X-Y plot of the optimized ENEEMMM trajectory.

This paper aims to complete the multi-target asteroid exploration mission with the minimum fuel consumption imparted by spacecraft. Thus, the optimal delta-V index can be chosen as J ¼ △ v0 þ △vn1 þ △vn2 þ △vega1 T þ △vf

tðkÞ ¼ tðk  1Þ þ 50 þ 650  xðkÞ ; k 2 f2; 4; 5g; ; In this case, the 50 target groups corresponding to the 50 optimal solutions are listed in Table 2, which are reserved to be further optimized for other five types of gravity-assist trajectories. In Table 2, NEA_num is the original number of all the 1081 candidate NEAs, C3 is the launch energy, J is the mission performance index, J00 is the total velocity increment required throughout the mission, i.e., J00 ¼ Jþ△v00-△v0, Tof is the total time of flight, and N_50 represents the ranking by J of all the 50 target asteroid groups. N_50 will also be used in the rest of the paper, with which the corresponding NEA_name, NEA_num, MBA_num and MBA_name can be obtained by referring to Tables 1–2. Additionally, there are 35 different NEAs involved in the 50 optimal solutions. Refer to

(27)

There are totally 5 optimal variables denoted as x(i) (i2[1,5], and i 2 Z*)2[0,1]. Meanwhile, the mission durations of Segment 0, 1 and 2 are all constrained in the range of [50,700] days. The variables required for patched conic method are expressed by these optimal variables   tð1Þ ¼ t0s þ t0f  t0s  xð1Þ ; tð3Þ ¼ tð2Þ þ 50 þ 100  xð3Þ ; 14

P. Sun et al.

Planetary and Space Science 182 (2020) 104851

The optimal delta-V index can be chosen as

Table 12 Optimized results for the ENEEMMM trajectory. variable

value

variable

value

variable

value

t(1)

{2025-0517} 10.716

rp1

231069.527

rp3

21800.202

t(5)

{2028-0319} 1.284

t(8)

{2032-0704} 1.106

C3 t(2) △vn1 t(3) △vn2 t(4) △vega1T

{2025-1208} 0.262

△vdsm1

{2026-0130} 0.804 {2026-1018} 2.074E-06

△vega2 T

t(6)

rp2 t(7) △vmga1T

{2029-0102} 7.72E-10 14288.658 {2029-1213} 0.000677

△vdsm2 t(9) △vmga2T rp4 t(10) △vf

J ¼ △ v0 þ △vn1 þ △vn2 þ △vega1T þ △vdsm1 þ △vega2T þ △vf

There are totally 12 optimal variables denoted as x(i) (i2[1,12], and i 2 Z*)2[0,1]. Meanwhile, the mission durations of Segment 0, 1, 2, 3 and 4 are all constrained in the range of [50,700] days. The 3 variables for patched conic method including t(1), t(2) and t(3) are the same as those for the ENEM trajectory. The remaining other 9 variables are expressed by these optimal variables

{2034-0321} 0.180

tðkÞ ¼ tðk  1Þ þ 50 þ 650  xðkÞ ; k 2 ½4; 7; k 2 Z ; The magnitude of Δvega1T: △vega1T ¼ 1:5  xð8Þ; Two angle variables for Δvega1T: beta1 ¼ 2π  xð9Þ, beta2 ¼ π ðxð10Þ  0:5Þ; Two gravity-assist angles of the EGA1: ψ E ¼ 2π  xð11Þ, δE ¼ δEmax  xð12Þ;

3695 {2036-0208} 2.170

Appendix A for the specific orbital elements of the 35 NEAs. Also note that the total delta-V provided by spacecraft of all the ENEM trajectories for the 50 optimal target groups is less than 6.77 km/s, which means that the final spacecraft mass mf after reaching the target MBA is greater than 278.2 kg. That is, the optimized results may be feasible and significant for the multi-target asteroid exploration mission using traditional chemical propulsion.

3.3. ENEMM In this case, one Mars gravity assist is considered. The patched conic method for the ENEEM trajectory is illustrated in Fig. 7, where t(5) is the time of the MGA1, and t(6) is the final MBA rendezvous time. The mission solving process can be decomposed into ct_seg ¼ 7 segments. Segment 0, 1, 2 and 3 are solved by the Lambert method. Segment 4 denotes the NEA sampling process. Segment 5 and 6 denote the processes of the EGA1 and the MGA1, respectively. The rp1 and rp2 are defined as the flyby radii of the EGA1 and the MGA1, respectively. The optimal delta-V index can be chosen as

3.2. ENEEM Considering that the spacecraft and the target MBA may not be the right places to enable a low arrival delta-V after the EGA1, DSM1 is assumed between the two Earth gravity assists. Hence, the total delta-V may be decreased with the help of the DSM1. The patched conic method for the ENEEM trajectory is illustrated in Fig. 6, where t(5) is the time of the DSM1, t(6) is the time of the EGA2, and t(7) is the final MBA rendezvous time. The mission solving process can be decomposed into ct_seg ¼ 8 segments. Segment 0, 1 and 5 are the same as the Segment 0, 1 and 3 for the ENEM trajectory, respectively. Segment 3 and Segment 4 are solved by the Lambert method. Segment 7 denotes the solving process of the EGA2. The transfer leg of Segment 2 is solved by the orbit propagation (OP in Fig. 6, for short). The rp1 and rp2 are defined as the flyby radii of the EGA1 and the EGA2, respectively.

a)

(28)

J ¼ △ v0 þ △vn1 þ △vn2 þ △vega1T þ △vmga1T þ △vf

(29)

There are 6 optimal variables denoted as x(i) (i2[1,6], and i 2 Z*)2 [0,1]. Meanwhile, the mission durations of Segment 0, 1, 2, 3 and 4 are all constrained in the range of [50,700] days. The 3 variables for patched conic method including t(1), t(2) and t(3) are the same as those for the ENEM trajectory. The remaining other 3 variables are expressed by these optimal variables tðkÞ ¼ tðk  1Þ þ 50 þ 650 ⋅ xðkÞ ; k 2 ½4; 6; k 2 Z ;

The semi-major axis a

b)

The inclination i

Fig. 20. a) The semi-major axis, and b) inclination for the optimized ENEEMMM trajectory over time. 15

P. Sun et al.

Planetary and Space Science 182 (2020) 104851

3.4. ENEEMM

3.6. ENEEMMM

In this case, both one EGA and one MGA are considered. The patched conic method for the ENEEMM trajectory is illustrated in Fig. 8, where t(5) is the time of the DSM1, t(6) is the time of the EGA2, t(7) is the time of the MGA1, and t(8) is the final MBA rendezvous time. The mission solving process can be decomposed into ct_seg ¼ 10 segments. Segment 0, 1, 3, 4 and 5 are solved by the Lambert method. The transfer leg of Segment 2 is solved by the orbit propagation. Segment 6 denotes the NEA sampling process. Segment 7, 8 and 9 denote the processes of the EGA1, the EGA2 and the MGA1, respectively. The rp1, rp2 and rp3 are defined as the flyby radii of the EGA1, the EGA2 and the MGA1, respectively. The optimal delta-V index can be chosen as

In this case, both one EGA and dual Mars gravity assists are introduced. The patched conic method for the ENEEMMM trajectory is illustrated in Fig. 10, where t(5) is the time of the DSM1, t(6) is the time of the EGA2, t(7) is the time of the MGA1, t(8) is the time of the DSM2, t(9) is the time of the MGA2, and t(10) is the final MBA rendezvous time. The mission solving process can be decomposed into ct_seg ¼ 13 segments. Segment 0, 1, 3, 4, 6 and 7 are solved by the Lambert method. The transfer legs of Segment 2 and Segment 5 are solved by the orbit propagation. Segment 8 denotes the NEA sampling process. Segment 9, 10, 11 and 12 denote the processes of the EGA1, the EGA2, the MGA1 and the MGA2, respectively. The rp1, rp2, rp3 and rp4 are defined as the flyby radii of the EGA1, the EGA2, the MGA1 and the MGA2, respectively. The optimal delta-V index can be chosen as

J ¼ △ v0 þ △vn1 þ △vn2 þ △vega1T þ △vdsm1 þ △vega2T þ △vmga1T þ △vf

J ¼ △ v0 þ △vn1 þ △vn2 þ △vega1T þ △vdsm1 þ △vega2T þ △vmga1T

(30)

þ △vdsm2 þ △vmga2T þ △vf

There are totally 13 optimal variables denoted as x(i) (i2[1,13], and i 2 Z*)2[0,1]. Meanwhile, the mission durations of Segment 0, 1, 2, 3 and 4 are all constrained in the range of [50,700] days. The 3 variables for patched conic method including t(1), t(2) and t(3) are the same as those for the ENEM trajectory. The remaining other 10 variables are expressed by these optimal variables

(32) There are totally 20 optimal variables denoted as x(i) (i2[1,20], and i 2 Z*)2[0,1]. The mission durations of Segment 0, 1, 2, 3, 4, 6 and 7 are all constrained in the range of [50,700] days, while the mission duration on the MGA1 to DSM1 transfer leg of Segment 5 is constrained in the range of [350,1350] days. The 3 variables for patched conic method including t(1), t(2) and t(3) are the same as those for the ENEM trajectory. The remaining other 17 variables are expressed by these optimal variables

tðkÞ ¼ tðk  1Þ þ 50 þ 650  xðkÞ ; k 2 ½4; 8; k 2 Z ; Two gravity-assist angles of the EGA1: ψ E ¼ 2π  xð9Þ, δE ¼ δEmax  xð10Þ; The magnitude of Δvega1T: △vega1T ¼ 1:5  xð11Þ; Two angle variables for Δvega1T: beta1 ¼ 2π  xð12Þ, ...;

tðkÞ ¼ tðk  1Þ þ 50 þ 650  xðkÞ ; k 2 f4; 5; 6; 7; 9; 10g; ; tð8Þ ¼ tð7Þ þ 350 þ 1000  xð8Þ ;

3.5. ENEMMM

The magnitude of Δvega1T: △vega1T ¼ 1:5  xð11Þ; Two angle variables for Δvega1T: beta1 ¼ 2π  xð12Þ, beta2 π ðxð13Þ  0:5Þ; Two gravity-assist angles of the EGA1: ψ E ¼ 2π  xð14Þ, δE δEmax  xð15Þ; The magnitude of Δvmga1T: △vmga1T ¼ 1:5  xð16Þ; Two angle variables for Δvmga1T: beta3 ¼ 2π  xð17Þ, beta4 π ðxð18Þ  0:5Þ; Two gravity-assist angles of the MGA1: ψ M ¼ 2π  xð19Þ, δM δMmax  xð20Þ;

In this case, dual Mars gravity assists are considered. The patched conic method for the ENEEMM trajectory is illustrated in Fig. 9, where t(5) is the time of the MGA1, t(6) is the time of the DSM1, t(7) is the time of the MGA2, and t(8) is the final MBA rendezvous time. The mission solving process can be decomposed into ct_seg ¼ 10 segments. Segment 0, 1, 2, 4 and 5 are solved by the Lambert method. The transfer leg of Segment 3 is solved by the orbit propagation. Segment 6 denotes the NEA sampling process. Segment 7, 8 and 9 denote the processes of the EGA1, the MGA1 and the MGA2, respectively. The rp1, rp2 and rp3 are defined as the flyby radii of the EGA1, the MGA1 and the MGA2, respectively. The optimal delta-V index can be chosen as

4. Results and discussion

J ¼ △ v0 þ △vn1 þ △vn2 þ △vega1T þ △vmga1T þ △vdsm1 þ △vmga2T

4.1. Data quality validation and preliminary analysis

¼ ¼

¼ ¼

þ △vf Table 3 and Table 4 list the global minimum delta-V results of all the other five types of gravity-assist trajectories for the 50 optimal asteroid groups in Table 2. Table 5 and Table 6 show the comparison results of the postlaunch delta-V Δvp for some same candidate MBAs displayed in the work of (Chen et al., 2014). DT, MGA, MMGA, EGA and EMGA correspond to the direct transfer trajectory and four types of gravity-assist trajectories to MBAs in the work of Chen et al. respectively. In this paper, Δv*p is defined as Δv*p ¼ J-(Δv0þΔvn1þΔvn2þΔvega1T), and Evp is defined as Evp ¼ |Δv*p -Δvp |/Δvp, denoting the relative deviation between Δvp and Δv*p. It can be computed that the average relative deviation values for the five cases compared are 0.0621, 0.140, 0.2646, 0.1284 and 0.0867, respectively. Note that the value of Evp ¼ 1.038 for Flora is much larger than that for any other MBA within the same comparison group of ENEMMM and MMGA. The reason for the significant difference between the results is due to the intrinsic population-based stochastic characteristic of the PSO algorithm, the different performance index to be optimized and the different mission assumptions given in the two papers.

(31) There are totally 13 optimal variables denoted as x(i) (i2[1,13], and i 2 Z*)2[0,1]. The mission durations of Segment 0, 1, 2, 4 and 5 are all constrained in the range of [50,700] days, while the mission duration on the MGA1 to the DSM1 transfer leg of Segment 3 is constrained in the range of [350,1350] days, considering that the orbital period of Mars is about 686.98 days. The 3 variables for patched conic method including t(1), t(2) and t(3) are the same as those for the ENEM trajectory. The remaining other 10 variables are expressed by these optimal variables tðkÞ ¼ tðk  1Þ þ 50 þ 650 ⋅ xðkÞ ; k 2 f4; 5; 7; 8g; ; tð6Þ ¼ tð5Þ þ 350 þ 1000 ⋅ xð6Þ ; Two gravity-assist angles of the MGA1: ψ M ¼ 2π ⋅ xð9Þ, δM ¼ δMmax ⋅ xð10Þ; The magnitude of Δvmga1T: △vmga1T ¼ 1:5 ⋅ xð11Þ; Two angle variables for Δvega1T: beta1 ¼ 2π ⋅ xð12Þ, beta2 ¼ π ðxð13Þ  0:5Þ; 16

P. Sun et al.

Planetary and Space Science 182 (2020) 104851

there are totally 31 feasible examples found, which is the most among the 6 cases. 18 feasible examples are obtained for the ENEEMMM trajectory, while the total number of feasible solutions gained for the remaining other four types of trajectories is 19 (see Appendix B). In particular, the global optimal solution is acquired for the ENEMMM trajectory. ENEMMM: 31 feasible examples.

Considering that the comparison results for other cases are basically the same, this outlier of Evp ¼ 1.038 will not affect the overall comparison of optimization results. Table 7 shows the comparison results of the mission duration TEM on the Earth to MBA transfer leg (the EGA1 to MBA transfer leg in this paper) for the 5 same candidate MBAs both in the work of Chen et al. and this paper. The comparison value TEM corresponds to the TOF in Tables 4 and 5 in the work of Chen et al. In this paper, TEM corresponds to tEM ¼ Tof-(tof0þtof1þtsam). The definitions of Tof, tof0, tof1 and tsam are the same as those defined in Section 3. Etem ¼ |TOF-tEM|/TOF characterizes the relative deviation between the mission duration from Earth to the target MBA in the two papers. As shown in Table 7, for the cases analyzed, the relative deviation in the total flight time of the MBA rendezvous mission is less than 1%, validating the optimality and feasibility of the results in this paper. Figs. 11 and 12 illustrate the optimized minimum, median, and maximum value of the launch energy C3, performance index J, time of flight Tof and terminal rendezvous delta-V Δvf of the 50 target asteroid groups for each trajectory type. As shown in Fig. 11a, the overall launch energy for the ENEEMMM trajectory, median value being 15.11 km2/s2, is obviously smaller than that for the other five types of trajectories. However, the overall launch energy for the ENEM trajectory is high. The result is consistent with the work of (Chen et al., 2014), i.e., the gravity of Mars is helpful to deliver a spacecraft from the Earth to the MBA with low launch energy whether the subtask of NEA sample return exists in the whole mission. Fig. 11b shows that the single technique of one Earth gravity assist or one Mars gravity assist does not work well in reducing the total mission cost, while the overall fuel consumption imparted by the spacecraft for the ENEMMM and ENEEMMM trajectory is smaller than that of other cases. Also, the minimum performance index for the ENEMMM trajectory reaches the global minimum value of 3.4952 km/s. However, the time of flight for the ENEMMM trajectory or the ENEEMMM trajectory is almost twice or three times of that for the ENEM trajectory. The existence of the time-of-flight penalty may cause an increase in the total mission cost. Fig. 12b shows that both the techniques of one MGA and dual MGAs enable a low arrival delta-V to reach the orbit of the target MBA, which is consistent with the work of (Chen et al., 2014). To sum up, the trajectory type of ENEMMM may offer better solutions for a multi-target asteroid exploration mission including the subtask of NEA sample return and the subtask of MBA rendezvous. In this case, the total fuel consumption can be much lower. However, it cannot be ignored that a time-of-flight penalty exists and the mean time of flight obtained in this paper for the ENEMMM trajectory is greater than 2500 days.

4.3. A sample return mission to 2014 JR24 and rendezvous mission to Urania For the feasible lists provided in Section 4.2, it is shown that there are feasible examples under the asteroid group of 2014 JR24 (NEA_num: 552) and Urania (MBA_num: 7) for the ENEMMM and ENEEMMM trajectory. Also note that the example mission targeting the asteroid group of 2014 JR24 and Urania for the ENEMMM trajectory is the global optimal solution among all the cases. Here the trajectory design results for the exploration mission including sample return of the near-Earth asteroid 2014 JR24 and rendezvous of the main-belt asteroid Urania are presented and analyzed. Case 1: ENEMMM J ¼ 3.495 km/s, Tof ¼ 2289 days. The ballistic transfer trajectory is depicted in Fig. 13, and the trajectory details can be found in Table 9. Fig. 14 shows both the X-Y and XZ projections of the same transfer trajectory. The mission duration (tof0þtsam þ tof1) for the sample return of 2014 JR24 is 559 days, which accounts for 24.41% of the total time of flight Tof. There is a 115 day sampling period before departing back to Earth. After the gravity assist of Earth (EGA1), it takes 1371 another days to reach the orbit of Urania. Figs. 15–17 show the semi-major axis a, the eccentricity e, the inclination i, the right ascension of ascending node (RAAN) Ω, the argument of periapsis ω, and the true anomaly f over time, for the optimized ENEMMM transfer trajectory, respectively. Fig. 15a presents that both the two MGAs and the EGA1 contribute to the significant increase in the semi-major axis, helping to increase the orbit energy to reach the target MBA. As shown in Fig. 16a, the orbital inclination is substantially increased by the EGA1, while both the MGA1 and MGA2 help to decrease the inclination to that of the Urania orbit. That is, both the two MGAs and the EGA1 serve to adjust the vertical orientation of the orbital plane. Fig. 15b shows that the transfer trajectory gets more oblate after the MGA2. Figs. 16b and 17a illustrate that the MGA1 has the greatest influence on the right ascension of ascending node and the argument of periapsis, playing a significant role in adjusting the horizontal orientation of the orbital plane. Fig. 17b shows that the DSM1 is executed in the vicinity of the aphelion. The multi-target asteroid exploration mission is assumed to launch just at 1 AU from the sun, with a minimum distance of nearly 0.83 AU on the sample return leg (from 2014 JR24 to Earth) and a maximum distance of nearly 2.66 AU near the end of the mission. The spacecraft has a maximum distance from Earth of 3.299 AU, but the spacecraft is on the same side of the sun during the sampling phase of the mission with an Earth-Spacecraft distance of 2.78 AU at arrival and 2.64 AU at departure. Similar to the work of (Dankanich et al., 2010), the spacecraft distances from the sun, Earth, Mars, 2014 JR24 and Urania are shown in Fig. 18a. Fig. 18b illustrates that the spacecraft gets the most velocity increment of 3.0758 km/s by EGA1. J ¼ 5.806 km/s, Tof ¼ 3919 days. Fig. 19 shows the 3-D and X-Y projections of the optimized ENEEMMM transfer trajectory. The trajectory details can be found in Table 12. Also note that the units of time, velocity are the same as those defined before. The semi-major axis a and the inclination i over time for the same transfer trajectory are shown in Fig. 20. Fig. 20a shows that both the EGA2 and the two MGAs contribute to increasing the semi-major axis to reach that of the Urania orbit. Fig. 20b presents that EGA1 helps reduce the orbital inclination to that of 2014 JR24, and the difference in

4.2. Feasible trajectory design results between 2022 and 2028 Here, the initial mass of the spacecraft is assumed to be m0 ¼ 2000 kg, the specific impulse of the chemical propulsion is set to be Isp ¼ 350 s, and the final spacecraft mass mf is required to be greater than mf0 ¼ 328.5 kg. Thus the performance index J to be minimized has to satisfy the following constraint  J < Isp g0 log10

m0 mf0

 ¼ 6:2 ðkm=sÞ

(33)

in which g0 is the standard sea level acceleration of gravity. Under the assumptions above, all feasible mission lists and corresponding feasible asteroid groups for the six types of trajectories are presented. Tables 8–9 and Tables 10–11 show the feasible mission lists and the 3 optimal mission examples for the ENEMMM trajectory and the ENEEMMM trajectory, respectively. The accessibility lists and corresponding optimal mission examples for the remaining other four types of trajectories are shown in Appendix B. In these tables, No. represents the ranking by the performance index J; N_50, NEA_num, and MBA_num have the same definitions as those given above; the units of time and velocity are UTCG and km/s, respectively. For the ENEMMM trajectory, 17

P. Sun et al.

Planetary and Space Science 182 (2020) 104851

rendezvous, apart from the direct transfer, EGA, MGA, EMGA and MMGA trajectories, an extra type of gravity-assist trajectory, EMMGA, is further investigated to explore the combined effect of EGA and MMGA. The PSO algorithm is employed to optimize the six types of gravity-assist trajectories. Through numerical simulation, the accessible asteroid group lists and feasible multi-target missions within the assumed 2022–2028 launch window are presented. The conclusions about the gravity-assisted multiasteroid exploration mission are as follows. For the subtask of NEA sample return, EGA1 not only helps to return the sample collected from the NEA as quickly as possible but also plays a critical role in adjusting the vertical orientation of the orbital plane. For the subtask of MBA rendezvous, dual Mars gravity assists are able to contribute to remarkable decrease in postlaunch delta-V as both the MGA1 and MGA2 make a tremendous contribution to the increase of the orbital energy, and hence reduce the terminal MBA rendezvous difficulty, while the single technique of one EGA or MGA does not work well, which is consistent with the conclusion in (Chen et al., 2014). The total mission cost imparted by spacecraft can be significantly reduced for the ENEMMM and ENEEMMM trajectory despite a time-of-flight penalty. The global minimum velocity increment of J ¼ 3.495 km/s is attained in the example mission of 2014 JR24 sample return and Urania rendezvous for the ENEMMM trajectory, in which only one EGA is needed.

the inclination between the Earth orbit and the Urania orbit is mostly reduced after EGA2 and MGA1. That is, both EGA2 and MGA1 play a fundamental role in adjusting the vertical orientation of the orbital plane. 4.4. Discussion Within the 2022–2028 launch window, the accessibility lists containing 68 feasible example missions using chemical propulsion are presented. This study has demonstrated that it may be possible to perform the mission of NEA sample return (via EGA) and MBA rendezvous using chemical rockets. There are different feasible solutions corresponding to the same target NEA and MBA, which are optimized for different types of gravity-assist trajectories. Also there are different feasible solutions corresponding to different target NEAs and MBAs for the same type of gravity-assist trajectory. The results found here can be taken as a guide in the selection of the feasible candidate asteroids and mission lists for future multi-target asteroid exploration missions. The global optimal solution is obtained for the ENEMMM trajectory in the mission of 2014 JR24 sample return and Urania rendezvous. Dual Mars gravity assists are required in the preliminary trajectory design, which reduces the total mission cost with low launch energy and terminal rendezvous delta-V. The result is consistent with the work of (Chen et al., 2014). As a result of the two numerical test cases above on the mission of 2014 JR24 sample return and Urania rendezvous, the second MGA (MGA2) shows better performances in the increase of the semi-major axis a and eccentricity e than the EGA1 for the ENEMMM trajectory (also better than the EGA2 for the ENEEMMM trajectory). However, it is important to mention that the effect of the EGA (corresponding to the EGA1 for the ENEMMM trajectory, and EGA1 and EGA2 for the ENEEMMM trajectory, respectively) on the orbital inclination is particularly evident. In fact, EGA not only shows significant influence on the orbital inclination, but also helps to return the sample collected from the near-Earth asteroid.

Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. CRediT authorship contribution statement Pan Sun: Conceptualization, Methodology, Software, Writing - original draft. Hongwei Yang: Methodology, Software, Writing - review & editing. Shuang Li: Conceptualization, Methodology, Writing - review & editing, Supervision.

5. Conclusion Acknowledgments A multi-target asteroid exploration mission of NEA sample return (via EGA) and MBA rendezvous using chemical propulsion and multiple gravity assists has been studied. The accessibility of 1081 candidate NEAs and 37 candidate MBAs with diameters larger than 100 km is evaluated by minimizing the velocity increments imparted by the spacecraft. A modified Tisserand graph with redefined tick marks is presented to give an intuitive and quick reference of possible gravity-assist sequences and the minimum number of gravity-assist bodies. For the subtask of MBA

This work was partially supported by the National Natural Science Foundation of China (Grant No. 11972182), sponsored by Qing Lan Project, Natural Science Foundation of Jiangsu Province (BK20180410), the Young Elite Scientists Sponsorship Program by CAST (2018QNRC001) and Funded by Science and Technology on Space Intelligent Control Laboratory (Grant No. KGJZDSYS-2018-11). The authors fully appreciate their financial supports.

Appendix A. The orbital elements of the 35 NEAs In this appendix, the specific orbital elements7 of the 35 NEAs chosen from the 50 optimal solutions for the ENEM trajectory are given in Table A. 1. Table A 1 NEAs chosen from the optimized ENEM trajectory [MJD: 58400] No.

NEA_num

a/AU

e

i/deg

Ω/deg

ω/deg

M/deg

1 2 3 4 5 6 7 8

119 552 440 343 179 469 222 536

0.943691 1.06645 1.09895 0.949421 0.933767 0.951627 1.08907 0.930751

0.164097 0.118351 0.118409 0.127384 0.140786 0.194972 0.146844 0.160615

1.294 0.93 2.337 1.155 3.532 1.488 0.855 1.31

233.325 48.89 186.55 198.101 297.666 43.112 321.5 138.966

332.279 246.427 42.257 304.996 301.044 119.44 242.82 260.325

271.221 309.387 208.608 118.413 257.452 38.633 87.841 154.059

(continued on next column)

7

https://newton.spacedys.com/neodys/[retrieved on 12 April 2019]. 18

P. Sun et al.

Planetary and Space Science 182 (2020) 104851

Table A 1 (continued ) No.

NEA_num

a/AU

e

i/deg

Ω/deg

ω/deg

M/deg

9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

266 501 380 39 512 453 726 251 694 1059 139 325 1009 348 337 976 672 445 422 568 41 833 384 812 582 802 963

1.03833 1.07819 1.03501 0.929333 1.05819 1.05411 0.917135 0.964063 0.999366 1.00307 0.945352 0.959323 0.930411 0.945551 1.06045 1.01896 1.07065 1.07376 1.0644 1.08524 0.940205 1.0632 1.09813 1.01969 1.06343 0.939931 1.08631

0.220677 0.226646 0.229307 0.168071 0.092285 0.094843 0.146505 0.125853 0.125196 0.116902 0.232746 0.144695 0.112119 0.094056 0.188281 0.138575 0.187149 0.120828 0.099355 0.184881 0.137591 0.221413 0.134086 0.251142 0.125382 0.211183 0.226315

3.49 0.107 0.277 0.892 7.341 5.839 0.777 1.908 6.357 1.672 2.063 0.72 3.707 0.371 0.883 1.385 6.58 2.206 2.645 2.832 5.773 4.372 2.817 1.504 8.16 5.756 0.475

88.519 3.858 346.426 13.404 11.912 142.911 205.039 189.48 136.802 353.002 48.221 37.002 2.046 132.479 163.844 236.246 266.382 222.673 316.917 335.989 115.446 354.412 349.974 75.546 206.549 127.816 42.359

263.982 210.719 71.341 132.165 77.278 50.014 205.91 291.204 119.066 270.645 116.174 316.278 137.684 73.12 267.558 147.44 86.329 96.704 255.039 264.452 243.301 84.159 117.596 159.264 101.388 119.575 74.906

163.966 296.054 156.096 273.06 150.504 7.818 86.413 111.775 137.86 128.935 290.426 234.248 246.383 62.556 80.465 359.842 263.47 198.254 320.406 303.098 261.434 225.473 260.937 153.519 289.136 213.256 216.142

Appendix B. The accessibility lists and optimal mission examples for four types of trajectories In this appendix, the 3 optimal mission examples for the ENEM trajectory are shown in Table B. 1. Table B. 2-3, Table B. 4-5, and Table B. 6-7 present the feasible mission lists and the 3 optimal mission examples for the ENEEM trajectory, the ENEMM trajectory and the ENEEMM trajectory, respectively. ENEM: 3 feasible examples. Table B 1 Example missions for the ENEM trajectory N_50

1

2

3

{NEA_num,MBA_num}

{119,24}

{552,5}

{440,19}

t(1) C3 t(2) △vn1 t(3) △vn2 t(4) △vega1T

{2026-12-28} 18 {2027-10-18} 0.676613307 {2028-03-16} 1.555100041 {2029-01-25} 0.001858871 6678 {2030-07-29} 3.68977593

{2025-05-12} 11.35618 {2026-02-24} 0.132528 {2026-04-15} 1.104272 {2026-11-18} 0.980596 6678 {2027-10-20} 3.823502

{2027-03-15} 15.99386 {2027-12-01} 0.380041 {2028-01-20} 1.420284 {2028-09-07} 0.138789 6678 {2029-08-13} 4.137281

rp1

t(5) △vf

ENEEM: 6 feasible examples.

Table B 2 The accessibility list for the ENEEM trajectory No.

N_50

NEA_num

MBA_num

J/km/s

Tof/day

1 2 3 4 5 6

5 4 14 26 33 8

343 343 469 139 337 469

11 19 19 11 24 5

4.943742 5.183609 5.587287 5.700038 5.88452 5.972297

1815.078 1737.98 2350.022 1872.185 1447.613 1415.365

19

P. Sun et al.

Planetary and Space Science 182 (2020) 104851

Table B 3 The 3 optimal example missions for the ENEEM trajectory N_50

5

4

14

{NEA_num,MBA_num}

{343,11}

{343,19}

{469,19}

t(1) C3 t(2) △vn1 t(3) △vn2 t(4) △vega1T rp1 t(5) △vdsm1 t(6) △vega2T rp2 t(7) △vf

{2022-07-16} 8.5492 {2022-11-14} 0.7112 {2023-01-31} 0.5384 {2023-11-02} 0.03395 7236.86 {2025-02-20} 0.0753 {2025-11-07} 8.5438E-8 6678.040 {2027-07-05} 3.585

{2022-07-16} 8.6412 {2022-11-01} 0.6992 {2023-01-19} 0.6279 {2023-11-04} 2.515E-6 9152.52 {2024-11-17} 0.0947 {2025-11-11} 0 6678.000 {2027-04-19} 3.762

{2024-11-19} 17.9982 {2025-09-26} 0.8497 {2026-01-24} 0.6821 {2026-11-21} 2.581E-6 6921.15 {2028-02-05} 0.0339 {2029-11-16} 5.9008E-6 9110.292 {2031-04-27} 4.022

ENEMM: 5 feasible examples.

Table B 4 The accessibility list for the ENEMM trajectory No.

N_50

NEA_num

MBA_num

J/km/s

Tof/day

1 2 3 4 5

24 8 2 25 22

343 469 552 1059 694

1 5 5 15 16

4.825886 5.288581 5.521142 5.829097 6.01089

1376.605 1399.405 1254.6 1932.77 1583.704

Table B 5 Three optimal example missions for the ENEMM trajectory N_50

24

8

{NEA_num,MBA_num}

{343,1}

{469,5}

2 {552,5}

t(1) C3 t(2) △vn1 t(3) △vn2 t(4) △vega1T rp1 t(5) △vmga1T rp2 t(6) △vf

{2022-10-19} 7.748242 {2023-08-20} 0.326063 {2023-10-29} 1.215239 {2024-10-19} 0.519469 6678 {2025-04-22} 0.552916 3695 {2026-07-27} 2.212199

{2024-11-21} 17.34178 {2025-09-13} 0.856296 {2025-12-08} 0.805096 {2026-11-12} 0 6829.723 {2027-04-12} 0 3772.022 {2028-09-20} 3.62719

{2025-05-11} 11.47579 {2026-03-04} 0.129075 {2026-05-04} 1.220301 {2026-11-26} 0.50494 6678 {2027-04-15} 3.13E-07 3695 {2028-10-17} 3.666825

ENEEMM: 5 feasible examples.

Table B 6 The accessibility list for the ENEEMM trajectory No.

N_50

NEA_num

MBA_num

J/km/s

Tof/day

1 2 3 4 5

25 32 29 18 50

1059 694 325 453 963

15 5 2 13 15

5.505304 5.984626 6.045876 6.05994 6.147016

2292.118 1878.044 2158.494 2533.665 2416.557

20

P. Sun et al.

Planetary and Space Science 182 (2020) 104851

Table B 7 Three optimal example missions for the ENEMM trajectory N_50

25

32

29

{NEA_num,MBA_num}

{1059,15}

{694,5}

{325,2}

t(1) C3 t(2) △vn1 t(3) △vn2 t(4) △vega1T rp1 t(5) △vdsm1 t(6) △vega2T rp2 t(7) △vmga1T rp3 t(8) △vf

{2026-03-19} 11.73001 {2027-01-17} 0.361365 {2027-06-05} 0.449278 {2028-03-01} 1.35E-07 144288.8 {2028-08-06} 1.090237 {2029-05-22} 0 6715.588 {2030-08-24} 0.352145 3695 {2032-06-27} 3.252279

{2026-02-07} 18 {2026-06-15} 0.520576 {2026-10-30} 1.007939 {2027-08-12} 3.59E-07 39773.93 {2028-04-30} 0.00846 {2029-01-20} 3.33E-06 17010.24 {2030-01-27} 1.572712 3695 {2031-03-31} 2.615843

{2026-04-19} 16.25926 {2026-10-10} 0.503624 {2026-12-16} 2.094875 {2027-12-24} 4.4E-05 9981.113 {2028-09-03} 0.10916 {2029-06-07} 0 10990.39 {2030-09-17} 0 4433.122 {2032-03-17} 3.338172

Olympio, J.T., 2011. Optimal control problem for low-thrust multiple asteroid tour missions. J. Guid. Contr. Dynam. 34, 1709–1720. https://doi.org/10.2514/1.53339. Perozzi, E., Rossi, A., Valsecchi, G.B., 2001. Basic targeting strategies for rendezvous and flyby missions to the near-Earth asteroids. Planet. Space Sci. 49, 3–22. https:// doi.org/10.1016/S0032-0633(00)00124-0. Pontani, M., Conway, B.A., 2010. Particle swarm optimization applied to space trajectories. J. Guid. Contr. Dynam. 33, 1429–1441. https://doi.org/10.2514/ 1.48475. Qiao, D., Cui, H., Cui, P., 2006. Evaluating accessibility of near-Earth asteroids via Earth gravity assists. J. Guid. Contr. Dynam. 29, 502–505. https://doi.org/10.2514/ 1.16757. Sears, D., Allen, C., Britt, D., Brownlee, D., Franzen, M., Gefert, L., Gorovan, S., Pieters, C., Preble, J., Scheeres, D., Scott, E., 2004. The Hera mission: multiple near-Earth asteroid sample return. Adv. Space Res. 34, 2270–2275. https://doi.org/10.1016/ j.asr.2003.05.059. Shang, H., Liu, Y., 2017. Assessing accessibility of main-belt asteroids based on Gaussian process regression. J. Guid. Contr. Dynam. 40, 1144–1154. https://doi.org/10.2514/ 1.G000576. Sims, J.A., 1997. Delta-V Gravity-Assist Trajectory Design: Theory and Practice. Ph.D. Dissertation. School of Aeronautics and Astronautics. Purdue Univ, West Lafayette, IN, pp. 106–108. Song, Y., Gong, S., 2019. Solar-sail trajectory design for multiple near-Earth asteroid exploration based on deep neural networks. Aero. Sci. Technol. 91, 28–40. https:// doi.org/10.1016/j.ast.2019.04.056. Strange, N.J., Longuski, J.M., 2002. Graphical method for gravity-assist trajectory design. J. Spacecraft Rockets 39, 9–16. https://doi.org/10.2514/2.3800. Tsuda, Y., Yoshikawa, M., Abe, M., Minamino, H., Nakazawa, S., 2013. System design of the Hayabusa 2—asteroid sample return mission to 1999 JU3. Acta Astronaut. 91, 356–362. https://doi.org/10.1016/j.actaastro.2013.06.028. Vasile, M., Pascale, P. De, 2006. Preliminary design of multiple gravity-assist trajectories. J. Spacecraft Rockets 43, 794–805. https://doi.org/10.2514/1.17413. Wilt, C.M., Thayer, J.T., Ruml, W., 2010. A comparison of greedy search algorithms. In: The Third Annual Symposium on Combinatorial Search. Association for the Advancement of Artificial Intelligence. Yamakawa, H., Kawakatsu, Y., Morimoto, M., Abe, M., Yano, H., 2004. Low-thrust trajectory design for a low-cost multiple asteroid rendezvous and sample return mission. In: The 55th International Astronautical Congress. IAC-04-Q, p. 20. https:// doi.org/10.2514/6.IAC-04-Q.P.20. Yang, H., Li, J., Baoyin, H., 2015. Low-cost transfer between asteroids with distant orbits using multiple gravity assists. Adv. Space Res. 56, 837–847. https://doi.org/ 10.1016/j.asr.2015.05.013. Yang, H., Tang, G., Jiang, F., 2018. Optimization of observing sequence based on nominal trajectories of symmetric observing configuration. Astrodynamics 2, 25–37. https:// doi.org/10.1007/s42064-017-0009-2. Yen, C.L., 1982. Main-belt asteroid exploration: mission options for the 1990s. In: The AIAA/AAS Astrodynamics Conference. AIAA Paper, 1982-1463. Zhang, X., Huang, J., Wang, T., Huo, Z., 2019. ZhengHe-A mission to a near-Earth asteroid and a main belt comet. In: The 50th Lunar and Planetary Science Conference.

References Bao, C., Yang, H., Barsbold, B., Baoyin, H., 2015. Capturing near-Earth asteroids into bounded Earth orbits using gravity assist. Astrophys. Space Sci. 360, 224–234. https://doi.org/10.1007/s10509-015-2581-3. Bayliss, S., 1971. Precision targeting for multiple swingby interplanetary trajectories. J. Spacecraft Rockets 8, 927–931. https://doi.org/10.2514/3.59749. Ceriotti, M., Vasile, M., 2010. Automated multigravity assist trajectory planning with a modified Ant Colony algorithm. J. Aero. Comput. Inf. Commun. 7, 261–293. https:// doi.org/10.2514/1.48448. Chapman, C.R., 1994. The Galileo encounters with Gaspra and Ida. In: SymposiumInternational Astronomical Union, pp. 357–365. https://doi.org/10.1017/ s0074180900046647. Chen, Y., Baoyin, H., Li, J., 2014. Accessibility of main-belt asteroids via gravity assists. J. Guid. Contr. Dynam. 37, 623–632. https://doi.org/10.2514/1.58935. Cheng, A.F., Santo, A.G., Heeres, K.J., Landshof, J.A., Farquhar, R.W., Gold, R.E., Lee, S.C., 1997. Near-earth asteroid rendezvous: mission overview. J. Geophys. Res.: Planets 102 (10), 23695–23708. https://doi.org/10.1029/96JE03364. Cui, P., Qiao, D., Cui, H., Luan, E., 2010. Target selection and transfer trajectories design for exploring asteroid mission. Sci. China Technol. Sci. 53, 1150–1158. https:// doi.org/10.1007/s11431-010-0007-6. Dankanich, J.W., Landau, D., Martini, M.C., Oleson, S.R., Rivkin, A., 2010. Main belt asteroid sample return mission design. In: The 46th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit. American Institute of Aeronautics and Astronautics. Garner, C.E., Rayman, M.M., Whiffen, G.J., Brophy, J.R., Mikes, S.C., 2013. Ion propulsion: an enabling technology for the Dawn mission. In: The 23rd AAS/AIAA Spaceflight Mechanics Meeting. American Astronomical Society. AAS 13-342. Green, S.F., Barucci, M.A., Michel, P., B€ ohnhardt, H., Brucato, J.R., Dotto, E., Ehrenfreund, P., Franchi, I.A., Lara, L.M., Marty, B., Koschny, D., Agnolon, D., Chalex, R., Martin, P., Romstedt, J., Andert, T.P., Cremonese, G., Groussin, O., Josset, J.L., 2013. MarcoPolo-R: sample return from NEA (341843) 2008 EV5. In: European Planetary Science Congress 2013, pp. EPSC2013–E2776. Jiang, F., Baoyin, H., Li, J., 2012. Practical techniques for low-thrust trajectory optimization with homotopic approach. J. Guid. Contr. Dynam. 35, 245–258. https:// doi.org/10.2514/1.52476. Kawaguchi, J.I., Fujiwara, A., Uesugi, T., 2008. Hayabusa—its technology and science accomplishment summary and Hayabusa-2. Acta Astronaut. 62, 639–647. https:// doi.org/10.1016/j.actaastro.2008.01.028. Kennedy, J., Eberhart, R., 1995. Particle swarm optimization. In: IEEE International Conference on Neural Networks. Institute of Electrical and Electronics Engineers Inc, Piscataway, NJ, USA. https://doi.org/10.1109/ICNN.1995.488968. Lau, C.O., Hulkower, N.D., 1989. Accessibility of near-Earth asteroids. J. Guid. Contr. Dynam. 10, 225–232. https://doi.org/10.2514/3.20207. Li, H., Chen, S., Izzo, D., Baoyin, H., 2019. Deep Networks as Approximators of Optimal Transfers Solutions in Multitarget Missions arXiv preprint arXiv:1902.00250. Morimoto, M., Yamakawa, M., Yoshikawa, M., Abe, M., Yano, H., 2004. Trajectory design of multiple asteroid sample return missions. Adv. Space Res. 34, 2281–2285. https:// doi.org/10.1016/j.asr.2003.10.055.

21