Multiobjective Mission Planning for Multiple Geosynchronous Spacecraft Refueling

. This paper studies the multiple geosynchronous spacecraft refueling problem (MGSRP) with multiple servicing spacecraft (Ssc) and fuel depots (FDs). In the mission scenario, multiple Ssc and FDs are parked in the geosynchronous Earth orbit (GEO) initially. Ssc start from FDs and maneuver to visit and refuel multiple GEO targets with known demands. These capacitated Ssc are expected to rendezvous with fuel-de ﬁ cient GEO targets and FDs for the purpose of delivering the fuel stored in FDs to GEO targets. The objective is to ﬁ nd a set of Pareto-optimal solutions with minimum fuel cost and mission duration. The MGSRP is a much more complex variant of multidepot vehicle routing problems mixing discrete and continuous variables. A two-nested optimization model is built. We propose a new multiobjective hybrid particle swarm optimization to solve the outer-loop problem, and the design variables are the refueling sequence, task assignment, time distribution, and locations of FDs. In the inner-loop problem, branch and bound method is used to ﬁ nd the optimal decision variable for a given outer-loop solution. Finally, numerical simulations are presented to illustrate the e ﬀ ectiveness and validity of the proposed approach.

The best economy will be achieved if multiple spacecraft can be serviced in a single mission [10,11]. Mission planning problem arising from servicing multiple spacecraft began to attract attention recently. Multispacecraft on-orbit service missions can be divided into three main types. The first type uses a single servicing spacecraft to service multiple targets. Cerf [12] investigated the space debris collecting mission aimed at removing 5 heavy low Earth orbit (LEO) debris per year to stabilize the debris population. A branch and bound algorithm was applied to optimize debris selec-tion and orbit maneuvers. In our previous research, we studied multiple geosynchronous Earth orbit (GEO) spacecraft inspection mission, in which the chaser spacecraft is parked in an equator-high elliptical orbit initially, and two maneuvers are exerted at perigee for each visual inspection. The objective is to optimize the visitation order and time regarding fuel cost [13]. Alfriend et al. [14] developed a method to optimize the order for visiting a set of GEO spacecraft with small inclinations, and finding the path of minimum fuel cost is similar to the travelling salesman problem (TSP). Zhang et al. [15] studied the optimization of a near circular LEO long-duration multiple spacecraft refueling mission considering the J 2 perturbation and time-window constraints with one active servicing spacecraft.
The second type uses multiple servicing spacecraft to service multiple targets. Bo and Feng [16] proposed a new scheme for refueling spacecraft constellation, and the new pattern is based on formation flying. Both single-supplier refueling and double-supplier refueling were studied. Yu et al. [17] studied the mission planning of GEO debris removal with multiple servicing spacecraft (Ssc). Considering it as a hybrid optimal control problem, a mathematical model was proposed. The third type uses a peer-to-peer (P2P) strategy in which all the spacecraft can be used as an active one or a passive one. Shen and Tsiotras [18] studied the P2P refueling problem, which seeks fuel equalization among the spacecraft in the same constellation. In the scenario, all the spacecraft were capable of visiting and refueling each other using two-impulse, multirevolution transfers. The objective function of the P2P refueling problem reflects a balance between achieving fuel equalization and minimizing the fuel cost. Later on, Dutta et al. [19] researched on the P2P refueling problem further by developing a nonlinear programming solver capable of determining low-thrust P2P maneuvers. In addition, several strategies and corresponding algorithms have been proposed, such as asynchronous P2P refueling, egalitarian P2P, and cooperative P2P refueling [20][21][22].
For many mission planning problems of on-orbit servicing, it is often necessary to optimize fuel cost and mission duration simultaneously that are generally conflicting with each other. Thus, there is not one unique optimal solution but a set of Pareto-optimal solutions [23][24][25]. Madakat et al. [26] proposed a biobjective time-dependent TSP model for the debris removal mission and used a branch and bound (B&B) method to deal with it. The GEO debris removal mission planning problem investigated by Yu et al. is also biobjective. Over the past 20 years, a variety of evolutionary algorithms have been developed to solve multiobjective problems, and most of them are able to find a set of Pareto-optimal solutions in a single run. The representative and widely used methods [27][28][29][30] are the Pareto archive evolutionary strategy (PAES), the nondominated sorting genetic algorithm (NSGA) and its improved version NSGA-II, multiobjective particle swarm optimization (MOPSO), and so on.
In this paper, we focus on multiple geosynchronous spacecraft refueling problems (MGSRP), and the reasons are as follows. Hundreds of high-value geosynchronous spacecraft are parked in the GEO, with an altitude of nearly 35,786 km, and have the same angular velocity with the Earth's [9]. Onboard fuel is a main factor that influences the lifetime of a GEO spacecraft because each GEO spacecraft requires approximately 52 m/s for station-keeping per year, and it may need to maneuver to cover a new location in which fuel cost is inevitable [6]. The lifespan can be extended if on-orbit refueling is conducted. In addition, on-orbit refueling is perceived to be of the lowest risk because it is generally carried out when the onboard fuel of a spacecraft is almost exhausted. Hence, on-orbit refueling of GEO spacecraft, as a means to extend their lifetimes rather than replacing them with new ones, is promising.
The refueling strategies can be divided into three main categories: one-to-many [15], many-to-many [31], and peer-to-peer [19]. For the one-to-many strategy, the number of fuel-deficient targets that could be refueled is quite small due to the limited fuel capacity of Ssc. Similarly, the manyto-many strategy will fail as well if the fuel demands of targets exceed the capacity of all the Ssc. The P2P strategy cannot be used for GEO spacecraft as not all of them have the ability of imparting and receiving fuel. In order to overcome the shortcomings of the above three strategies, multiple Ssc and fuel depots (FDs) are employed. The fuel stored in FDs is delivered to multiple GEO targets by Ssc. The purpose of this paper is to optimize a multiple GEO target refueling mission with the objective of minimizing fuel cost and mission duration at the same time. Four key problems which are location selection, task assignment, mission sequence, and time distribution are solved by a proposed approach using nested multiobjective hybrid particle swarm optimization (MOHPSO) and B&B.
The rest of this paper is organized as follows. In Section 2, we present a general description of the MGSRP with multiple Ssc and FDs. After that, in Section 3, a two-nested optimization model is built. Later on, we solved the transfer problem with a given solution. Then, MOHPSO and B&B are used to solve the MGSRP. Finally, four cases are employed to demonstrate the validity and effectiveness of the proposed two nested resolution method.

Mission Scenario
The multiple geosynchronous spacecraft refueling problem (MGSRP) involves three kinds of spacecraft, namely, servicing spacecraft (Ssc) that deliver fuel; fuel depots (FDs) that provide fuel; and fuel-deficient GEO spacecraft that accept fuel. N fuel-deficient GEO spacecraft with different inclinations, phase angles, and right ascension of ascending node (RAAN) angles are at the ends of lifetimes. In order to extend their lifetimes, m Ssc and p FDs are used to refuel these GEO targets. These capacitated Ssc have to move from the FDs and visit the n GEO targets out of N candidates for the purpose of delivering the fuel stored in the FDs to fueldeficient GEO targets. Each time the Ssc departs from an FD, the initial fuel mass should not exceed the fuel capacity C. It is assumed that idle Ssc is not permitted and the refueling mission is conducted in a noncooperative way, namely, only the Ssc perform orbit transfers. The Ssc will return to any FD whenever necessary to get replenished. Obviously, in the process of refueling mission, Ssc may be in one of five states: parked in GEO, in transit to a GEO target, refueling a GEO target, in transit to an FD, or receiving fuel from an FD. FDs may be in one of two states: parked in GEO or refueling an Ssc. In the mission scenario, each GEO target is refueled no more than once; however, the Ssc can be refueled by the FDs for several times. When the whole mission is completed, the Ssc should finally return to one of p FDs waiting for the next refueling mission. Obviously, the fuel cost is closely related to the locations of FDs. Hence, RAANs, orbit inclinations, and phase angles of FDs are continuous design variables.
The refueling mission of #i Ssc can be divided into several subtasks. In the jth GEO target refueling of #i Ssc, the target serial number and the duration of this subtask will be denoted by x i j ∈ f1, 2, ⋯, Ng and Δt i j , respectively. When a refueling process is completed, the Ssc can either go back to a FD and get refueled or transfer to refuel another GEO target. Let s i j ∈ f0, 1, ⋯, pg be the decision variable. s i j = k 2 International Journal of Aerospace Engineering (1 ≤ k ≤ p) represents that #i Ssc will return to #k FD and get refueled. s i j = 0 represents that it will maneuver to refuel the next target. Figure 1 illustrates an instance of the MGSRP problem with 10 candidate GEO targets (represented by solid circles). Two Ssc and two FDs are employed to refuel 8 GEO targets. It can be seen in Figure 1 that after #1 Ssc departs from #1 FD, it visits #3, #1, and #6 targets in order and then transfers to #2 FD and gets refueled. Next, #1 Ssc visit #4 and #7 targets in order. Finally, it transfers to #2 FD, waiting for another mission. #2 Ssc starts from #2 FD, visits #5, #9, and #2 targets in order, and finally comes back to #2 FD.
The MGSRP takes two objectives into account: minimizing fuel cost and duration of the refueling mission. The MGSRP does not have a single optimal solution because the design objectives are conflicting. Instead, there exists a set of solutions which are nondominated, known as Pareto-optimal solutions. Therefore, the mission planning problem is to find a set of Pareto-optimal solutions so as to grasp the tradeoff between the two objectives.
The archetype of the MGSRP is the multidepot vehicle routing problem (MDVRP), which is a variant of the basic vehicle routing problem (VRP). The VRP is concerned with delivering goods to a set of customers with known demands through vehicle routes that begin and finish at the depot with minimum cost. VRP and its variant MDVRP are classified as being NP-hard. The MGSRP is clearly NP-hard as a much more complex variant of the MDVRP. The differences between the MGSRP and the MDVRP are summarized here.
(1) Only n GEO targets among the N candidates must be visited To sum up, the MGSRP studied in this paper can be stated as follows: m Ssc and p FDs running on the GEO belt are employed to visit and refuel n GEO targets out of N, and the following key problems should be resolved: The details of rendezvous and dock, refueling operations are out of our considerations. The goal of this paper is to find a set of Pareto-optimal solutions with minimum fuel cost and mission duration.

Optimization Model
3.1. Design Variables. There are two types of design variables of the MGSRP. The first type consists of discrete variables X, ΔT, and S, given by where X i , ΔT i , and S i are the refueling order, given mission duration and decision variable of #i Ssc, respectively. The second type consists of continuous variables I, Ω, and L, given by where I, Ω, and L represent the orbit inclinations, RAANs, and longitudes of p FDs, respectively. Obviously, variable X identifies the targets that should be refueled by each Ssc, as well as the refueling order. It can also be represented by where X all represents a permutation of N GEO targets, and R gives the number of GEO targets assigned to each Ssc. In this paper, X all and R are optimized directly rather than X.

Objectives. The two objective functions are
where ∑ r i j=1 Δt i j represents the duration of refueling r i GEO targets by #i Ssc. As m Ssc could execute the refueling process simultaneously and independently, taking the maximum value of ∑ r i j=1 Δt i j , i = 1, 2, ⋯m indicates the total mission duration. Hence, the first objective is to minimize the total mission duration, and the second is to minimize the fuel cost Mfuel, given by where mfuel i h is the fuel received when the Ssc departs from the FD for the hth time, Mfuel i is the total fuel cost of #i Ssc, and H i is the number of fuel transfers between #i Ssc and FDs. The details of evaluating objective function f 2 are given in Section 4.

Constraints.
We introduce a lower bound ΔT min and an upper bound ΔT max on the given duration of each subtask: When the selected GEO targets of a Ssc have been refueled, the Ssc should finally return to one of the FDs; thus, Each Ssc should at least refuel one GEO target as idle Ssc is not permitted. In addition, a total of n GEO targets should be refueled. We have The fuel received from FD is limited by the capacity of Ssc.
3.4. Two-Nested Optimization Model. In this subsection, a two-nested optimization model is proposed for the MGSRP. The decision variable S is optimized in the inner-loop optimization problem. The objective of the inner-loop optimization is to minimize the total fuel cost. The variables X all , R, ΔT, I, Ω, and L are sought in the outer-loop optimization problem; the first objective is to minimize the duration of the refueling mission, and the second is to minimize the total fuel cost. The two-nested optimization model for the MGSRP is formulated as find X all , R, ΔT, I, Ω, L min f 1 , f 2 s:t:

Transfer Problem
In order to evaluate the fuel cost with a given solution, X, ΔT , S, I, Ω, and L, we should solve several transfer problems by optimizing the trajectory from one GEO target to another. In this section, the maneuver strategy of Ssc is described in detail, and we propose a simple method to obtain the optimal impulses of the Ssc, as well as the fuel weight each time it leaves the FDs.

Time Relationship.
Let t i j,0 , t i j,f be the initial and ending time of the jth subtask of #i Ssc, respectively. It follows that where ΔT i j is the actual duration of the subtask, including the orbit transfers between GEO targets and FDs and the operations applied to each of them, constrained by where Clearly, the initial time of jth subtask is the ending time of ðj − 1Þth subtask; thus, 4 International Journal of Aerospace Engineering

Transfer Strategy.
It is assumed that each Ssc is equipped with a high thrust engine and the maneuvers are impulsive. The total velocity increment (Δv) consists of two components, the plane change and in-plane transfer. Plane change impulse should adjust both the orbit inclination and RAAN of the Ssc at the same time. Two-impulse phasing maneuver is applied to in-plane transfer, which is used to adjust the Ssc's phase to match with that of the GEO target. The plane change impulse can be coupled with the first phasing change impulse, so that we essentially have a two-impulse transfer.
The procedure of a subtask can be classified into the following two cases, as well as the method to obtain the optimal impulses of the Ssc with a given solution.
Two impulses are exerted for Ssc to transfer to a GEO target, given by (in the following, superscript i and subscript j are omitted for simplicity).
where Δv out is exerted for plane change, Δv in1 , Δv in2 are the first and second impulse of phasing maneuvers, and Δv in2 = −Δv in1 . The magnitude of the two impulses can be computed by where θ is the angel between Δv out and Δv in1 . The velocity increment is Δv = jΔv 1 j + jΔv 2 j. Let I t , Ω t , f t be the orbit inclination, RAAN, and longitude of the GEO target, respectively. Let I s , Ω s , f s , u s be the orbit inclination, RAAN, longitude, and argument of latitude of Ssc, respectively. Denote the angle between the two orbit planes as γ (see Figure 2); then, we have where ΔΩ is the difference in RAAN, and ΔΩ ∈ ½0, π. For circular GEO, the velocity increment required for a plane change of γ is [14] Δv out j j= 2v sin where v is the orbit velocity of the GEO spacecraft. In this case, the actual duration of a subtask is made up of three parts, given by where ΔT ref is the refueling time of the GEO target and Δ T trans1 and ΔT trans2 are the orbital transfer time and awaiting time of Ssc, respectively. It can be seen from Figure 2 that the length of arc c (0 < c < π) is c = cos −1 −cos I t + cos γ cos I s sin γ sin I s : ð21Þ The first impulse Δv 1 is exerted when Ssc fly through the crossing point of Ssc orbit and target orbit. The awaiting time on Ssc's initial orbit ΔT trans2 could be computed by arc length c. If u s > c, then ΔT trans2 = ð2π − u s + c/2πÞT g , else ΔT trans2 = ðc − u s /2πÞT g , where T g is the orbit period of GEO. When the GEO target refueling is completed, u s can be computed by The longitude difference Δf is measured from the target to the Ssc, given by Δf = f s − f t , and Δf ∈ ½−π, π. If Δf > 0, then Δv in1 and Ssc's initial velocity are in the same direction, and θ = ðπ/2Þ + ðγ/2Þ, else they are in opposite direction, and θ = ðπ/2Þ − ðγ/2Þ.
Let N s be the number of phasing orbit revolutions for the Ssc, and let N g be the number of whole revolutions for the GEO target. N s and N g should satisfy where T s is the orbit period of the phasing orbit. Therefore, the orbit transfer time ΔT trans1 can be calculated by Δ T trans1 = N s T s . Instead of optimizing Δv in1 , Δv in2 directly, the optimal N s and N g are determined such that the Δv is minimized. Since Δv is minimized, the closer T s is to T g , the minimization problem is equivalent to minimizing jT s − T g j. According to Equation (23), we have Given that Δf varies from −π to π, it follows that N s and N g should be equal and jT s − T g j = jΔf /2πN s jT g . Otherwise, jT s − T g j = jðΔf + 2παÞ/2πN s jT g , and α is an integer. It is obvious that jΔf /2πN s jT g ≤ jΔf + 2πα/2πN s jT g . Obviously, the greater N s and N g , the smaller the Δv will be. However, they are also constrained by Equation (13), given by where floor (A) means rounding A to the nearest integers less than or equal to A. Once N s and N g are obtained, then Equation (24) could be used to determine T s , and then the semimajor axis of phasing orbit a s can be determined easily. The magnitude of phasing maneuvers can be obtained based on the vis-via equation, given in [32] Δv where μ is the Earth's gravitational constant and a g is the orbit radius of GEO.
Let I t , Ω t , and f t (I p , Ω p , and f p ) be the orbit inclination, RAAN, and longitude of the GEO target (FD), respectively. Let I s , Ω s , f s , and u s be the orbit inclination, RAAN, longitude, and argument of latitude of Ssc, respectively. Similar to case 1, we have where γ 1 , γ 2 are the angle between the two orbit planes of Ssc and target, target, and FD, respectively. In this case, the actual duration of a subtask is made up of four parts, given by where ΔT ref and ΔT reff are the refueling times of the GEO target and Ssc, respectively. and ΔT trans1 , ΔT trans2 , ΔT trans3 , and ΔT trans4 are the orbital transfer time and awaiting time for visiting the GEO target and the FD, respectively. Similar to case 1, the length of arcs c 1 and c 2 are c 1 = cos −1 −cos I t + cos γ 1 cos I s sin γ 1 sin I s , c 2 = cos −1 −cos I p + cos γ 2 cos I t sin γ 2 sin I t : ð30Þ If u s > c 1 , then ΔT trans2 = ð2π − u s + c 1 /2πÞT g , else Δ T trans2 = ðc 1 − u s /2πÞT g . When the GEO target refueling is completed, u s can be computed by If u s > c 2 , then ΔT trans4 = ð2π − u s + c 2 /2πÞT g , else Δ T trans4 = ðc 2 − u s /2πÞT g . When the Ssc is refueled by the FD, u s can be computed by The longitude difference Δf 1 , Δf 2 are defined by Δf 1 = f s − f t , Δf 2 = f t − f p , and Δf 1 , Δf 2 ∈ ð−π, πÞ. If Δf k > 0, then , θ k = ðπ/2Þ + ðγ k /2Þ, else θ k = ðπ/2Þ − ðγ k /2Þ (k = 1, 2).
The two transfer orbits that the Ssc used to visit the GEO target and the FD are called the phasing orbits 1 and 2, respectively. Let N s1 , N s2 be the number of revolutions of phasing orbits 1 and 2. Let N g and N p be the number of whole revolutions for the GEO target and the FD. N s1 , N s2 , N g , and N p should satisfy where T s1 , T s2 are the periods of phasing orbits 1 and 2, respectively. Therefore, the orbit transfer time ΔT trans1 and ΔT trans3 can be calculated by ΔT trans1 = N s1 T s1 , ΔT trans3 = N s2 T s2 . Similarly, instead of optimizing the impulses of the Ssc directly, the optimal N s1 , N s2 , N g , and N p are determined such that the Δv is minimized. The objective is to minimize the Δv, which is equivalent to minimizing jT s1 − T g j + jT s2 − T g j, given by International Journal of Aerospace Engineering Similar to case 1, given that Δf 1 and Δf 2 vary from −π to π, it follows that N s1 = N g and N s2 = N p . Then, we have Let N s = N s1 + N s2 , given by Inspecting Equation (35) shows that the minimization problem is merely a function of integers N s1 and N s2 ; therefore, it can be formulated as If we relax the integer constraint of N s1 and N s2 , the solution to Equation (37) can be obtained by setting the partial derivative of F = f + λðN s1 ′ + N s2 ′ − N s Þ with respect to N s1 ′ , N s2 ′ , and λ to zero.
Then, we have Besides, ∂ 2 F/∂N ′ 2 s1 > 0 and ∂ 2 F/∂N ′ 2 s2 > 0 are also guaranteed such that this solution corresponds to a minimum rather than a maximum. The possible optimal integer N s1 could be floorðN s1 ′ Þ or floorðN s1 ′ Þ + 1, and N s2 is given by N s2 = N s − N s1 . Therefore, the optimal N s1 and N s2 could be found easily by solving the minimization problem expressed by Equation (37). Once N s1 , N s2 , N g , and N p are obtained, then Equation (33) could be used to determine T s1 and T s2 .
To facilitate the readers to understand our proposed transfer strategy for case 2 (s i j > 0), the block diagram of the scheme to obtain Ssc's optimal impulses is given in Figure 3. The hth time of #i Ssc departing from the FD, let m i h,j be the remaining weight of the Ssc after refueling the jth GEO target, given by

Fuel
where m i h,0 is the initial weight of #i Ssc the hth time departing from the FD, δm i h,j is the fuel demand of the GEO target, I sp is the specific impulse, g 0 is the Earth gravity constant. Let m i h,f be the weight of the Ssc returning to a FD, given by For simplicity, and without loss of generality, it is assumed that the fuel is exhausted when the Ssc goes back to an FD. Denote the weight of the Ssc permanent structure as m dry ; then, m i h,f = m dry .Then, the fuel weight required to be stored in #i Ssc each time it departs from an FD is given  Figure 3: Schematic of the method to obtain the optimal impulses of Ssc.  (6). The pseudocode for calculating the total fuel weight of #i Ssc received from the FD is given in Algorithm 1, in which find (S i > 0) returns linear indices corresponding to the entries of S i that are greater than 0, and numelðsÞ returns the number of elements. If mfuel i h exceeds the capacity of the Ssc C, then the corresponding solution is infeasible, and Mfuel i will set to be a relatively high penalty cost.

Resolution Approach
The presence of two objectives in the outer-loop optimization gives rise to a set of Pareto-optimal solutions. One of these solutions cannot be said to be better than the other in terms of the two objectives. For this reason, we are required to find as many Pareto-optimal solutions as possible. A multiobjective hybrid particle swarm optimization (MOHPSO) is adopted to solve the outer-loop optimization.  where Rep½h is a value that is taken from the repository Rep and the selection method could be found in [27]. In each case, if the new particle is better than the old one, then it will be accepted. b. Update the particles in Rep: Insert all the currently non-dominated particles into Rep and eliminate the dominated particles from Rep. When the repository gets full, we apply a secondary criterion for retention: those located in less populated areas of objective space are given priority over those lying in highly populated regions. c. Update the memory of each particle: if the current particle is better than the one contained in its memory, it is updated using: Pbest½i = POP½i.  International Journal of Aerospace Engineering For the inner-loop optimization, we propose a branch and bound (B&B) method which could quickly locate the optimal S i in terms of Mfuel i with fewer computation costs than exhaustive search (ES).

Outer-Loop
Optimization Algorithm. A multiobjective particle swarm optimization (MOPSO) has been proposed in literature [27]. A MOHPSO algorithm, which is a combination of genetic algorithm (GA) and MOPSO, is employed to deal with the outer-loop optimization. Figure 4 shows a pictorial overview of the MOHPSO. In MOHPSO, particles are updated by genetic operators, i.e., selection, crossover, and mutation. X all , R, ΔT, I, Ω, and L are concatenated to form a particle. The main steps of MOHPSO are given in Algorithm 2.
In each fuel cost evaluation, it is required to have an inner-loop problem solver that provides the optimal fuel cost for each solution generated during the optimization of the MOHPSO outer-loop problem solver and that returns the fuel cost to the outer-loop problem solver.

Inner-Loop Optimization Algorithm.
To evaluate the fuel cost of a given outer-loop solution, the inner-loop problem of finding the optimal decision variable S i of #i Ssc should be solved. A straightforward approach to optimize S i would be an exhaustive search (ES), that is, the enumeration of all the feasible solutions and the choice of the best one. The pseudocode of ES is presented in Algorithm 3. r, c are the number of rows and columns of set.
The number of feasible decision variable S i amounts to pðp + 1Þ r i −1 , and the total number of feasible solutions of the inner-loop problem is ∑ m i=1 pðp + 1Þ r i −1 . ES becomes rapidly infeasible because computation time grows exponentially with the instance size. Therefore, ES is practically restricted only to small-size problems.
By a deep insight into the inner-loop problem of the MGSRP, we present a B&B method for finding the optimal S i , which is much more efficient than ES. In the B&B algorithm, the solution space of each node is divided into a number of subspaces. If it can be established that a subspace cannot contain the optimal solution, the whole subspace is discarded, else it is stored in the pool of live nodes. The search terminates when there is no unexplored part of the solution space left, and the optimal solution is then the one recorded as the "current best".
The B&B method is presented in Algorithm 4. If there is only one target, then Set = ½1 ; 2;⋯;p, calculate their corresponding fuel costs and find the optimal S i ; otherwise, S i is optimized as follows. Initially, Set = ½Set 0 ; Set 1 ;⋯;Set p and Set k = k (k = 0, 1, ⋯, p). If c < r j − 1, Set k is expanded, Set k = ½Set, k r×1 . The partial solutions in Set 0 are not comparable and the best partial solutions S optk in Set k could be found. Set is updated by Set = ½Set 0 ; S opt1 ;⋯;S optp . This process is repeated until c = r j − 1 , and then Set is expanded to Set = ½ Set 1 ;⋯;Set p , and Set k = ½Set, k r×1 (k = 1, ⋯, p). Finally, the optimal full solution could be found in Set. 6. Numerical Simulations and Analysis 6.1. Problem Configuration. To demonstrate the effectiveness and validity of the proposed method, four different cases are studied in this section. In the four cases, one Ssc and one FD, one Ssc and two FDs, two Ssc and one FD, and two Ssc and two FDs are used to refuel eight out of ten GEO targets, respectively. In cases 1 and 2, #1 Ssc starts from #1 FD. In case 3, both #1 and #2 Ssc start from #1 FD. In case 4, #1 and #2 Ssc start from #1 FD and #2 FD, respectively. It is assumed that the given duration of each subtask is an integer and the bound on it is Δt i j ∈ ½7, 21 days. The orbit parameters and fuel demands of ten GEO targets are given in Table 1. The mass parameters for the Ssc are m dry = 500 kg and C = 1000 kg. The initial u s of each Ssc is 0 deg. The thruster parameter is I sp g 0 = 3200 m/s. In the numerical simulations, the MOHPSO used a population of 100 particles, a repository size of 100 particles, and an iteration number of 100.

Comparision of ES and B&B.
The inner-loop problem of #i Ssc is solved by ES and B&B and both of them are sure to find the optimal S i . Figure 5 compares the computation time of ES and B&B on a 3.2 GHZ Intel Core i5-4460 personal computer with a memory of 4 GB. It can be seen that the time consumed by B&B is much shorter than that of ES, especially for a large number of GEO targets. Therefore, it can be concluded that B&B is able to deal with the innerloop problem at a relatively low computational cost. 1: Set k = k (k = 0, 1, ⋯, p); 2: while true do 3: Set = ½Set 0 ; Set 1 ;⋯;Set p , ½r, c = sizeðSetÞ; 4: if r i == c then 5: break; 6: else 7: Set k = ½Set, k r×1 , k = 0, 1, ⋯, p; 8: end if 9: end while 10: discard the solutions in Setthat did not satisfy the constraint s i r i > 0; 11: calculate Mfuel i by Algorithm 1, find the optimal S i in Set.
Algorithm 3: ES for finding the optimal S i . 9 International Journal of Aerospace Engineering 6.3. Pareto Solution Sets. Considering the stochastic characteristic of the MOHPSO, 30 independent runs for each case are executed, and all of the Pareto-optimal solutions obtained in 30 runs are compared and revised by deleting the dominated and repeating solutions. Then, the revised Pareto-optimal solutions of 30 runs are regarded as the final Pareto-optimal solution set. For cases 1-4, each run takes about 7 min, 24 min, 3.5 min, and 10 min, respectively. The Pareto-optimal solutions of four cases are shown in Figure 6. It can be seen that the fuel cost could be greatly reduced by extending the mission time for all cases. It is clear that case 4 performs better than cases 1-3 in overall fuel cost and mission duration, without taking the extra cost of Ssc and FD themselves into account. From the results of four cases, we may also conclude that with the increment of the amount of FDs (Ssc), the fuel cost (mission duration) decreases.
For cases 1 and 2 (3 and 4), Pareto-optimal solutions with 80 (45) days of mission time are listed in Table 2. The optimal solution of the 45-day MGSRP with two Ssc and two FDs is described in detail in the following. #1 and #2 FDs are located at the positions with inclinations 5.57 deg, 1.54 deg, RAANs 105.12 deg, 242.19 deg, and longitudes 141.06 deg, 179.48 deg, respectively. As is shown in Figure 10, the first time #1 Ssc departs from #1 PD, it maneuvers to refuel #7 target, and then it returns to #1PD and gets replenished. When #1Ssc departs from #1 PD for the second time, it visits and refuels #1, #3, and #10 targets in order. The first time #2 Ssc departs from #2 PD, it   10 International Journal of Aerospace Engineering   11 International Journal of Aerospace Engineering maneuvers to refuel #2, #4, and #5 targets in order, and then it returns to #2PD and gets replenished. When #2Ssc departs from #2 PD for the second time, it visits and refuels #8 target. At last, #1 Ssc and #2 Ssc return to #1 FD and #2 FD, respectively. Each time #1 Ssc departs from #1 FD, the amount of fuel stored is 250.9 kg and 778.2 kg, respectively. Each time #2 Ssc departs from #2 FD, the amount of fuel stored is 848.9 kg and 284.8 kg, respectively. All this fuel is delivered to GEO targets and consumed by orbital transfers. The optimal time distribution ΔT 1 of #1 Ssc is 10, 12, 12, and 11, and the actual ending times of each subtask are 9.32, 21.71, 33.73, and 44.33 days, respectively. The optimal time distribution ΔT 2 of #2 Ssc is 11, 10, 13, and 11, and the actual ending times of each subtask are 10.08, 20.42, 33.36, and 44.34 days, respectively. The optimal maneuver solutions of #1 and #2 Ssc are shown in Figure 11. #1 Ssc exerts

12
International Journal of Aerospace Engineering 12 impulsive maneuvers to visit 4 GEO targets and #1 FD, in which 4 maneuvers are exerted to go back to #1 FD. #2 Ssc exerts 9 impulsive maneuvers to visit 4 GEO targets and #2 FD, in which 4 maneuvers are exerted to go back to #2 FD.

Conclusions
In this paper, we investigate the MGSRP with multiple Ssc and FDs. We propose a two nested-loops optimization model and MOHPSO and B&B are employed to solve the outer-loop and inner-loop problems, respectively. Simulation results show that B&B is capable of finding the optimal inner-loop variable with much less computation cost than ES. Four different cases are studied in the numerical simulations, namely the multiple GEO targets refueling task is accomplished with one Ssc and one FD, one Ssc and two FDs, two Ssc and one FD, and two Ssc and two FDs. It can be seen that cases with two FDs consume less fuel than those with one FD for the same mission duration, whereas cases with two Ssc consume less time than those with one Ssc for the same fuel cost. Case 4 performs better than cases 1-3 in overall fuel cost and mission duration. However, it is worth mentioning that if we take the extra cost of Ssc and FD themselves into account, two Ssc and two FDs may not be the best choice.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no competing interests. 13 International Journal of Aerospace Engineering