Abstract

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-deficient GEO targets and FDs for the purpose of delivering the fuel stored in FDs to GEO targets. The objective is to find 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 find the optimal decision variable for a given outer-loop solution. Finally, numerical simulations are presented to illustrate the effectiveness and validity of the proposed approach.

1. Introduction

In the past few years, on-orbit servicing compromising refueling, repairs, and upgrades has acquired significant attention because it would enhance the operational life and capability of space systems [18]. Long et al. [9] summarized five broad categories of benefits of servicing: (1) reduce risk of mission failure, (2) reduce mission cost, (3) increase mission performance, (4) improve mission flexibility, and (5) enable new missions.

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 selection 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 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 [2022].

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 [2325]. 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 [2730] 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 many-to-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.

2. 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. 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, Ssc and FDs are used to refuel these GEO targets. These capacitated Ssc have to move from the FDs and visit the GEO targets out of candidates for the purpose of delivering the fuel stored in the FDs to fuel-deficient GEO targets. Each time the Ssc departs from an FD, the initial fuel mass should not exceed the fuel capacity . 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 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 # Ssc can be divided into several subtasks. In the th GEO target refueling of # Ssc, the target serial number and the duration of this subtask will be denoted by and , 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 be the decision variable. () represents that # Ssc will return to # FD and get refueled. 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(2)The fuel cost is time-dependant and related to the transfer strategy(3)The initial locations of FDs are not fixed, and they are instead design variables(4)The GEO targets are moving, while, of course, the customers are stationary(5)There are infinite possible trajectories to transfer from one GEO target to another(6)The objective is to minimize the fuel cost and mission duration simultaneously

To sum up, the MGSRP studied in this paper can be stated as follows: Ssc and FDs running on the GEO belt are employed to visit and refuel GEO targets out of , and the following key problems should be resolved: (1)Location Selection. Selecting the initial locations of FDs(2)Task Assignment. Identifying the GEO targets that should be visited by each Ssc(3)Mission Sequence. Deciding the visitation order for each Ssc (including visiting FDs).(4)Time Distribution. Deciding the travel time shared for each orbit transfer

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.

3. Optimization Model

3.1. Design Variables

There are two types of design variables of the MGSRP. The first type consists of discrete variables , , and , given by where , , and are the refueling order, given mission duration and decision variable of # Ssc, respectively. The second type consists of continuous variables , , and , given by where , , and represent the orbit inclinations, RAANs, and longitudes of FDs, respectively.

Obviously, variable identifies the targets that should be refueled by each Ssc, as well as the refueling order. It can also be represented by where represents a permutation of GEO targets, and gives the number of GEO targets assigned to each Ssc. In this paper, and are optimized directly rather than .

3.2. Objectives

The two objective functions are where represents the duration of refueling GEO targets by # Ssc. As Ssc could execute the refueling process simultaneously and independently, taking the maximum value of 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 , given by where is the fuel received when the Ssc departs from the FD for the th time, is the total fuel cost of # Ssc, and is the number of fuel transfers between # Ssc and FDs. The details of evaluating objective function are given in Section 4.

3.3. Constraints

We introduce a lower bound and an upper bound 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 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 is optimized in the inner-loop optimization problem. The objective of the inner-loop optimization is to minimize the total fuel cost. The variables , , , , , and 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

4. Transfer Problem

In order to evaluate the fuel cost with a given solution, , , , , , and , 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.

4.1. Time Relationship

Let , be the initial and ending time of the th subtask of # Ssc, respectively. It follows that where 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 th subtask is the ending time of th subtask; thus,

4.2. Transfer Strategy

It is assumed that each Ssc is equipped with a high thrust engine and the maneuvers are impulsive. The total velocity increment () 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.

Case 1. .

Two impulses are exerted for Ssc to transfer to a GEO target, given by (in the following, superscript and subscript are omitted for simplicity). where is exerted for plane change, , are the first and second impulse of phasing maneuvers, and . The magnitude of the two impulses can be computed by where is the angel between and . The velocity increment is .

Let , , be the orbit inclination, RAAN, and longitude of the GEO target, respectively. Let , , , 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 . For circular GEO, the velocity increment required for a plane change of is [14] where 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 is the refueling time of the GEO target and and are the orbital transfer time and awaiting time of Ssc, respectively. It can be seen from Figure 2 that the length of arc () is

The first impulse is exerted when Ssc fly through the crossing point of Ssc orbit and target orbit. The awaiting time on Ssc’s initial orbit could be computed by arc length . If , then , else , where is the orbit period of GEO. When the GEO target refueling is completed, can be computed by

The longitude difference is measured from the target to the Ssc, given by , and . If , then and Ssc’s initial velocity are in the same direction, and , else they are in opposite direction, and .

Let be the number of phasing orbit revolutions for the Ssc, and let be the number of whole revolutions for the GEO target. and should satisfy where is the orbit period of the phasing orbit. Therefore, the orbit transfer time can be calculated by .

Instead of optimizing , directly, the optimal and are determined such that the is minimized. Since is minimized, the closer is to , the minimization problem is equivalent to minimizing . According to Equation (23), we have

Given that varies from to , it follows that and should be equal and . Otherwise, , and is an integer. It is obvious that . Obviously, the greater and , the smaller the will be. However, they are also constrained by Equation (13), given by where floor () means rounding to the nearest integers less than or equal to .

Once and are obtained, then Equation (24) could be used to determine , and then the semimajor axis of phasing orbit can be determined easily. The magnitude of phasing maneuvers can be obtained based on the vis-via equation, given in [32] where is the Earth’s gravitational constant and is the orbit radius of GEO.

Case 2. .

In this case, a total of four impulses are exerted of Ssc to transfer to a GEO target and an FD, given by where , are exerted for plane change, , , and are phasing maneuvers, and , . Denote , as the angel between and and and , respectively. The velocity increment is .

Let , , and (, , and ) be the orbit inclination, RAAN, and longitude of the GEO target (FD), respectively. Let , , , and be the orbit inclination, RAAN, longitude, and argument of latitude of Ssc, respectively. Similar to case 1, we have where , 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 and are the refueling times of the GEO target and Ssc, respectively. and , , , and 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 and are

If , then , else . When the GEO target refueling is completed, can be computed by

If , then , else . When the Ssc is refueled by the FD, can be computed by

The longitude difference , are defined by , , and . If , then,, else ().

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 , be the number of revolutions of phasing orbits 1 and 2. Let and be the number of whole revolutions for the GEO target and the FD. , , , and should satisfy where , are the periods of phasing orbits 1 and 2, respectively. Therefore, the orbit transfer time and can be calculated by , .

Similarly, instead of optimizing the impulses of the Ssc directly, the optimal , , , and are determined such that the is minimized. The objective is to minimize the , which is equivalent to minimizing , given by

Similar to case 1, given that and vary from to , it follows that and . Then, we have

Let , given by

Inspecting Equation (35) shows that the minimization problem is merely a function of integers and ; therefore, it can be formulated as

If we relax the integer constraint of and , the solution to Equation (37) can be obtained by setting the partial derivative of with respect to , , and to zero.

Then, we have

Besides, and are also guaranteed such that this solution corresponds to a minimum rather than a maximum. The possible optimal integer could be or , and is given by . Therefore, the optimal and could be found easily by solving the minimization problem expressed by Equation (37). Once , , , and are obtained, then Equation (33) could be used to determine and .

To facilitate the readers to understand our proposed transfer strategy for case 2 (), the block diagram of the scheme to obtain Ssc’s optimal impulses is given in Figure 3.

4.3. Fuel Cost

Once the th () fuel transfer between # Ssc and FDs is finished, the Ssc will visit and refuel GEO targets and then return to an FD. At the th time of # Ssc departing from the FD, the velocity increment required for the Ssc to visit the th GEO target and return to an FD are denoted by () and , respectively. For a given solution, and could be calculated by the method given in Subsection 4.2.

The th time of # Ssc departing from the FD, let be the remaining weight of the Ssc after refueling the th GEO target, given by where is the initial weight of # Ssc the th time departing from the FD, is the fuel demand of the GEO target, is the specific impulse, is the Earth gravity constant. Let 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 ; then, .Then, the fuel weight required to be stored in # Ssc each time it departs from an FD is given by , and the total fuel cost could be obtained by Equation (6). The pseudocode for calculating the total fuel weight of # Ssc received from the FD is given in Algorithm 1, in which find () returns linear indices corresponding to the entries of that are greater than 0, and returns the number of elements. If exceeds the capacity of the Ssc , then the corresponding solution is infeasible, and will set to be a relatively high penalty cost.

1: ;
2: , , ;
3: for h =1 to
4: , ;
5: for j = Q to 1
6:  ;
7: end for
8: ;
9: if
10:  ;
11: else
12:  , break; {Si is infeasible}
13: end if
14: end for

5. 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. For the inner-loop optimization, we propose a branch and bound (B&B) method which could quickly locate the optimal in terms of with fewer computation costs than exhaustive search (ES).

5.1. 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. , , , , , and are concatenated to form a particle. The main steps of MOHPSO are given in Algorithm 2.

1) Initialization.
 a. Initialize the population POP.
 b. Initialize the memory of each particle, .
 c. Evaluate each of the particles in POP.
 d. Store the non-dominated solutions in the repository Rep.
2) WHILE maximum number of cycles has not been reached, do
 a. Update particles using crossover and mutation:
  ;
  ;
  ;
  where 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: .
 d. Increment the loop counter.
3) END WHILE.
4) RETURN Rep.

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.

5.2. 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 of # Ssc should be solved. A straightforward approach to optimize 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. , are the number of rows and columns of .

1: ();
2: while true do
3: , ;
4: if then
5:  break;
6: else
7:  , ;
8: end if
9: end while
10: discard the solutions in Setthat did not satisfy the constraint ;
11: calculate Mfueli by Algorithm 1, find the optimal Si in Set.

The number of feasible decision variable amounts to , and the total number of feasible solutions of the inner-loop problem is . 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 , 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 , calculate their corresponding fuel costs and find the optimal ; otherwise, is optimized as follows. Initially, and (). If , is expanded, . The partial solutions in are not comparable and the best partial solutions in could be found. is updated by . This process is repeated until , and then is expanded to , and (). Finally, the optimal full solution could be found in Set.

1: (), ;
2: while true do
3: ;
4: if then
5:  calculate Mfueli by Algorithm 1, find the optimal Sopt in Set;
6:  return Sopt ; break;
7: else
8:  if then
9:   ;{k denotes a column vector ()}
10:   calculate Mfueli by Algorithm 1, find the best Soptk in Setk;
11:   update Set by ;
12:  else
13:   expand Set by , ();
14:   calculate Mfueli by Algorithm 1, find the best Sopt in Set;
15:   return Sopt ; break;
16:  end if
17: end if
18: end while

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 days. The orbit parameters and fuel demands of ten GEO targets are given in Table 1. The mass parameters for the Ssc are kg and kg. The initial of each Ssc is 0 deg. The thruster parameter is 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.

6.2. Comparision of ES and B&B

The inner-loop problem of # Ssc is solved by ES and B&B and both of them are sure to find the optimal . 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 inner-loop problem at a relatively low computational cost.

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. Figures 710 show the angular momentum projections and longitudes of the given 10 GEO targets, as well as the selected Pareto-optimal solutions of four cases. Considering the length of this paper, the details of the solutions of cases 1-3 are not provided.

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 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 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 ΔT2 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 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.

7. 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.

Acknowledgments

This work was partially supported by the Construction and Scientific Research Project of Zhejiang Province under Grant (number 2021K132), the Professional Development Program for Visiting Scholars in Zhejiang Province Colleges and Universities under Grant (number FX2021158), and the Key Research Project of Zhejiang Lab under Grant (number G2021NB0AL03).