# Mission Planning of GEO Active Debris Removal Based on Revolver Mode

Authors: Yong Zhao; Yi Cao; Yang Chen; Zhijun Chen; Yuzhu Bai
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

The mission planning of active debris removal (ADR) of revolver mode on geosynchronous orbit (GEO) is studied in this paper. It is assumed that there are one service satellite, one space depot, and some pieces of space debris in the ADR mission. The service satellite firstly rendezvouses with the debris and then releases the thruster deorbit kits (TDKs), which are carried with the satellite, to push the debris to the graveyard orbit. Space depot will provide replenishment for the service satellite. The purpose of this mission planning is to optimize the ADR sequence of the service satellite, which represents the chronological order, in which the service satellite approaches different debris. In this paper, the mission cost will be stated firstly, and then a mathematical optimization model is proposed. ADR sequence and orbital transfer time are used as designed variables, whereas the fuel consumption in the whole mission is regarded as objective for optimizing, and a specific number of TDKs is also a new constraint. Then, two-level optimization is proposed to solve the mission planning problem, which is low-level for finding optimal transfer orbit using accelerated particle swarm optimization (APSO) algorithm and up-level for finding best mission sequence using immune genetic (IGA) algorithm. Numerical simulations are carried out to demonstrate the effectiveness of the model and the optimization method. Results show that TDK number influences the fuel consumption through impacting the replenishing frequency and TDK redundancy. To reduce fuel consumption, the TDK number should be optimized and designed with suitable replenishing frequency and minimum TDK redundancy.


---

## Body

## 1. Introduction

Due to the human space activities and active spacecraft increasing, the quantities of failed/decayed satellites and space debris grow at the same time and become the major source of the space pollutions [1]. According to the reports from Space Surveillance Net (SSN), more than 18,000 on-orbit space targets have been detected and catalogued, and 90% of them are space debris [2]. Moreover, the geosynchronous orbit (GEO) is an important orbital resource but crowded by debris with different shape and size, which will threaten the safety of on-orbit active spacecraft. By 2016, there are 1,484 objects on GEO, and only 31.7% of them are active satellites. According to European Space Agency (ESA), the collision probability between the debris and satellites will reach 3.7% till 2030 if the effective debris removal measure is not taken [3]. In recent decades, the active debris removal (ADR) has been come out and regarded as an effective measure for this terrible situation. There are some methods to realize ADR mission: space robot arm catching [4], space rope net catching [5], increasing the resistance [6], and the method of laser cleaning [7]. The ADR requires that one service satellite rendezvouses with different on-orbit debris and remove them [8, 9]. Obviously, according to the background of ADR mission, it is a multispacecraft rendezvous problem scientifically. The mission planning method can be raised as an important research topic for this scientific problem, which could reduce the fuel consumption and also find the rational sequence for multispacecraft rendezvous in the ADR mission.Research on mission planning in ADR task has been investigated for several years. Braun [10] analyzed different types of propulsive systems carried by service satellite and obtained the optimal mission sequence by exhaustive method. Considering the fuel and time constraints, Madakat [11] optimized the mission planning for debris removal on low earth orbit (LEO). Olympio [12] studied the optimal sequence for removing debris on sun-synchronous orbit for low-thrust orbit maneuver, and the results showed that the proposed method can reduce fuel consumption effectively under the influence of J2 perturbation. Yu [13, 14] studied the debris removal sequence on GEO and LEO with multiple service satellites. The previous literatures studied the mission planning in different debris removal scenarios with the orbit at different heights. Mostly, they focused on the phoenix mode [15] for removing debris; that is, the service satellite dragged debris to graveyard orbit firstly and then maneuvered to next debris and continued this action sequence. In this paper, the revolver mode [15] is utilized to realize ADR, the mission scenario of which can be shown in Figure 1. The scenario includes one service satellite, one space depot, and some space debris on GEO. Different from phoenix mode, the revolver mode avoids the service satellite transfer between graveyard orbit and GEO-debris’ orbit. The service satellite carries several thruster deorbiting kits (TDKs) [16] and attaches them on the surface of debris. Then, the TDK will open the propeller and push the debris to the graveyard orbit. In this way, the service satellite only needs to maneuver among the target debris, and the fuel consumption can be effectively lessened.Figure 1
Mission scenario of active debris removal.Mission planning for revolver mode brings new problem: the TDK is a new constraint. If fuel or TDK is used, the ADR mission will fail. This paper quotes the concept of space depot in ADR mission to provide replenishment for service satellite. The orbital depot, which is used to store TDK and fuel, is defined as a system that enables refueling of spacecraft elements in space [17]. When depleted of the TDKs, service satellite will return to the space depot for replenishment. One of the most prominent counterarguments to the propellant depot is centered on the thermal management of cryogenic fuels, and numerous studies and experiments have been done to solve it [18]. Recently, the propellant depot in space has proved practicable. Since the launch in December 2018, NASA’s Robotic Refueling Mission 3 (RRM3) demonstrated the first ever long-term storage of cryogenic fluid with zero boil off, having successfully stored cryogenic fluid for four months [19]. Under the circumstances of revolver mode, the number of TDKs and the sequence of the service’s maneuver should be optimized obviously. The purpose of this paper is to find the optimal sequence in the revolver mode for ADR mission and analyze the influence of the number of TDKs on the mission planning.So far, several reported works have focused on the optimization of mission planning. Some approaches are the nested loops, including the outer-loop to optimize the mission sequence and inner loop to find the optimal transfer orbit. Wall and Conway [20] used genetic algorithm (GA) to solve the double outer-loop and inner-loop problem for ADR mission, and it was proved to be efficient. Ross and D’Souza [21] introduced an HOC method to solve the space mission planning problem and proposed a formalism that can free mission planners to focus on high-level decision making by automating and optimizing the details of the inner loops. Du [22] used the MIGA-SQP algorithm to tackle the nested loops in optimal scheduling of multispacecraft refueling. Motivated by Du’s work of nested loops optimization, this paper integrates accelerated particle swarm optimization (APSO) algorithm with immune genetic algorithm (IGA) for mission planning problem. The IGA algorithm and APSO algorithm are used to solve outer loop and inner loop, respectively. In the previous study, PSO was a stochastic population-based optimization method proposed by Kennedy and Eberhart [23]. PSO is a promising evolutionary algorithm with a strong global search capability and simple effective operators [24]. The particle position represents the value of designed variables and particle velocity means the changing rate of variables during optimization process. The optimal solution is found through the moving of particle. Moreover, APSO algorithm aims to fully utilize the global exploration of PSO and the rapid local exploitation of the gradient-based method to speed up the global search process [25]. The IGA is the genetic algorithm inspired by the biological immune system. It can effectively deal with mission planning problem [26]. Compared with genetic algorithms, IGA has good global searching capability through inhibition of repetitive antibodies. Meanwhile, good antibodies are preserved by memory cells to obtain good convergence [27, 28].In this paper, we aim to investigate the mission planning problem considering the constrain of fuel consumption and number of TDKs. Compared to existing results, the major contributions and differences of this study are summarized as follows: (1) the mission planning for ADR mission in revolver mode is studied, and the mathematical model considering the number constraint of TDKs is modelled; (2) two-level optimization is organized to solve mission planning problem, low-level for finding optimal transfer orbit using APSO and up-level for finding best removing sequence using IGA.This paper is organized as follows. The mission scenario is described firstly, and mission costs are stated in Section2. Then, ADR optimization method is organized to seek the optimal solution in Section 3. Two-level model came out for optimization, including APSO algorithm for low-level optimization and IGA for the up-level. Following that, numerical simulations are illustrated in Section 4 to demonstrate the efficiency of the proposed solution method. Influence on fuel consumption of different quantities of TDKs carried by service satellite is discussed in this part as well. Finally, some appropriate conclusions are presented in Section 5.


## 2. Mission Scenario and Mathematical Model

### 2.1. Mission Scenario

There are three different types of spacecraft in revolver mode, including service satellite, space depot, and space debris. As shown in Figure1, ignoring the process of launching, the whole process for ADR mission can be divided into four steps. Firstly, the service satellite maneuvers to the position of the first debris. Then, the satellite extends the robotic arms to grab the debris and installs TDK on the surface of debris. Next, two robotic arms are retracted, one for grabbing space debris and the other for attaching TDK on debris. Finally, TDK pushes the debris to graveyard orbit and service satellite maneuvers to the next debris until all pieces of debris are removed. When the fuels of the service or the number of TDKs are insufficient, service satellite will move to space depot for replenishment. In the mission, space depot is located on fixed location in orbit and does not maneuver due to its huge mass and limited fuel.Summing up, the background of revolver mode studied in this paper is that substantial pieces of debris existed on GEO and one service satellite running among them. The service satellite needs to rendezvous with each piece of debris by sequences and attach TDKs to remove them. The removing sequence and transfer orbit need to be planned and optimized to reduce the fuel consumption.


### 2.2. Two-Impulse Rendezvous Model

The two-body model and Earth J2000 coordinate are utilized to depict the relative and absolute motion of service satellite around the Earth. It should be noted that the third-body perturbation and solar pressure have significant influence on the inclination and eccentricity of GEO target, respectively, the fact of which has been shown in relevant references [29, 30]. In the analysis of two references [31, 32], scholars pointed out that, in a short period (within one year), the impact of perturbation on GEO target is very small and can be ignored, and the optimal rendezvous order does not change. Therefore, the impact of atmospheric drag can be neglected, while since the ADR task is set within 5 days, it is assumed that the influence of other perturbations is ignored, including third-body perturbation and solar radiation pressure. Hence, it may be reasonable to neglect perturbations for simplifying the orbital model. Furthermore, the complexity of the problem will increase if the uncertainty of debris is taken into consideration. Therefore, the perturbations and position error of debris are ignored in this paper.Two-impulse rendezvous maneuvering model is adopted. Given the orbital elements of the two targetsE1,E2 and the orbital transfer time t, the transfer orbit between the two targets can be obtained by solving the loopy Lambert problem [33].Then, the impulse required for the orbital transfer is(1)ΔV=lambertE1,E2,t,


### 2.3. Mission Costs

In the revolver mode, fuel, time, and number of TDKs are the main costs. Suppose that the quantity of debris to be cleaned isn and the TDKs are labelled as numbers 1,2,3,…,n, respectively. DefineΔTi, ΔVi, andΔMi as the transfer time, impulse, and fuel consumption referring to i-th debris, respectively.During the process of ADR mission, two cases will occur:The first case is that fuels and quantities of TDKs are sufficient, and the service maneuvers to next debris directly. Fuel consumption for removingi-th debris can be calculated by the Tsiolkovsky formula:(2)ΔMi=Mi1−exp−ΔViIspgwhereMiis the initial mass of the service before removing i-th debris. Ispdenotes specific impulse, and g is the gravitational constant.The second case is that fuel or TDKs are insufficient. Service satellite should be transferred to space depot for replenishment and then maneuver to the next debris. Fuel consumption is obtained by(3)ΔMi=ΔMi−1s+ΔMsi,whereΔMi−1s and ΔMsi represent the fuel consumption of transferring to depot and to i-th debris, respectively. The reference [16] gives full analysis of how much fuel is required to place an object with the TDK into the graveyard orbit and weight; therefore, the fuel consumption constraint of TDK meets the result of the argument in this paper.Denote the time for debris removing and replenishing asTg and Tf, respectively. The time cost can be given as(4)ΔTi=ΔTi+Tg+q⋅Tf,where q∈0,1 is the binary decision parameter, while q=1 represents that service satellite transfers to space depot for replenishment, and q=0 represents satellite transfers to the next debris directly.Definek as the quantity of TDKs the service satellite carries. Suppose that each debris should be removed by using one TDK. Therefore, the TDK is reduced by n after removing all the debris.In summary, the fuel, time, and TDK consumption can be calculated as follows.(5)ΔM=∑i=1nΔMi,ΔT=∑i=1nΔTi,Δk=n.


## 2.1. Mission Scenario

There are three different types of spacecraft in revolver mode, including service satellite, space depot, and space debris. As shown in Figure1, ignoring the process of launching, the whole process for ADR mission can be divided into four steps. Firstly, the service satellite maneuvers to the position of the first debris. Then, the satellite extends the robotic arms to grab the debris and installs TDK on the surface of debris. Next, two robotic arms are retracted, one for grabbing space debris and the other for attaching TDK on debris. Finally, TDK pushes the debris to graveyard orbit and service satellite maneuvers to the next debris until all pieces of debris are removed. When the fuels of the service or the number of TDKs are insufficient, service satellite will move to space depot for replenishment. In the mission, space depot is located on fixed location in orbit and does not maneuver due to its huge mass and limited fuel.Summing up, the background of revolver mode studied in this paper is that substantial pieces of debris existed on GEO and one service satellite running among them. The service satellite needs to rendezvous with each piece of debris by sequences and attach TDKs to remove them. The removing sequence and transfer orbit need to be planned and optimized to reduce the fuel consumption.


## 2.2. Two-Impulse Rendezvous Model

The two-body model and Earth J2000 coordinate are utilized to depict the relative and absolute motion of service satellite around the Earth. It should be noted that the third-body perturbation and solar pressure have significant influence on the inclination and eccentricity of GEO target, respectively, the fact of which has been shown in relevant references [29, 30]. In the analysis of two references [31, 32], scholars pointed out that, in a short period (within one year), the impact of perturbation on GEO target is very small and can be ignored, and the optimal rendezvous order does not change. Therefore, the impact of atmospheric drag can be neglected, while since the ADR task is set within 5 days, it is assumed that the influence of other perturbations is ignored, including third-body perturbation and solar radiation pressure. Hence, it may be reasonable to neglect perturbations for simplifying the orbital model. Furthermore, the complexity of the problem will increase if the uncertainty of debris is taken into consideration. Therefore, the perturbations and position error of debris are ignored in this paper.Two-impulse rendezvous maneuvering model is adopted. Given the orbital elements of the two targetsE1,E2 and the orbital transfer time t, the transfer orbit between the two targets can be obtained by solving the loopy Lambert problem [33].Then, the impulse required for the orbital transfer is(1)ΔV=lambertE1,E2,t,


## 2.3. Mission Costs

In the revolver mode, fuel, time, and number of TDKs are the main costs. Suppose that the quantity of debris to be cleaned isn and the TDKs are labelled as numbers 1,2,3,…,n, respectively. DefineΔTi, ΔVi, andΔMi as the transfer time, impulse, and fuel consumption referring to i-th debris, respectively.During the process of ADR mission, two cases will occur:The first case is that fuels and quantities of TDKs are sufficient, and the service maneuvers to next debris directly. Fuel consumption for removingi-th debris can be calculated by the Tsiolkovsky formula:(2)ΔMi=Mi1−exp−ΔViIspgwhereMiis the initial mass of the service before removing i-th debris. Ispdenotes specific impulse, and g is the gravitational constant.The second case is that fuel or TDKs are insufficient. Service satellite should be transferred to space depot for replenishment and then maneuver to the next debris. Fuel consumption is obtained by(3)ΔMi=ΔMi−1s+ΔMsi,whereΔMi−1s and ΔMsi represent the fuel consumption of transferring to depot and to i-th debris, respectively. The reference [16] gives full analysis of how much fuel is required to place an object with the TDK into the graveyard orbit and weight; therefore, the fuel consumption constraint of TDK meets the result of the argument in this paper.Denote the time for debris removing and replenishing asTg and Tf, respectively. The time cost can be given as(4)ΔTi=ΔTi+Tg+q⋅Tf,where q∈0,1 is the binary decision parameter, while q=1 represents that service satellite transfers to space depot for replenishment, and q=0 represents satellite transfers to the next debris directly.Definek as the quantity of TDKs the service satellite carries. Suppose that each debris should be removed by using one TDK. Therefore, the TDK is reduced by n after removing all the debris.In summary, the fuel, time, and TDK consumption can be calculated as follows.(5)ΔM=∑i=1nΔMi,ΔT=∑i=1nΔTi,Δk=n.


## 3. ADR Optimization Method

On the basis of the maneuvering strategy developed in Section2, the mission planning problem for ADR can be treated as an optimization problem for the purpose of minimizing the total fuel consumption.


### 3.1. Optimization Model

#### 3.1.1. Design Variables

For mission planning problem, it is necessary to determine the ADR sequence [34]. The design variables consist of removal order X, binary decision variable S, and time distribution T. The X and S are integers, and T is a continuous number, respectively. The variables are of the form(6)X=x1,x2,…,xn,S=s1,s2,…,sn,T=Δt1,Δt2,…Δtm,tw1,tw2,…,twm,For example, the removing orderX=2,5,4,3,1 denotes that 5 pieces of debris need to be cleaned, and service satellite removes the debris as the order 2 ⟶ 5 ⟶ 4 ⟶ 3 ⟶ 1. S=0,0,1,0,1means that service satellite transfers to space depot for replenishment before removing the 3-th and 5-th debris. In T, Δt1,Δt2,…Δtmindicate the time period for each maneuver, and tw1,tw2,…,twmare the waiting time before orbit transfer.In the ADR mission, the objective is to minimize the total fuel consumption,MinΔM.


#### 3.1.2. Constraints

(1) Fuel Constraint. The service satellite cannot maneuver if the remaining fuel is insufficient, and the removing task will be interrupted. The constraint for fuel can be stated as(7)Mfi−ΔMi≥0,i=1,2,…,k−1,where Mfiis remaining fuel before removing the i-th debris.(2) TDK Number Constraint. TDK is being consumed during the mission. Therefore, the number of TDKs is required to be sufficient and can be stated as(8)kleft,i≥0,where kleft,i is the number of left TDKs the service satellite carries before removing the i-th debris.(3) Time Distribution Constraint. The ADR mission usually requires timeliness. The constraint is(9)Δti+twi≤tmax.The integrated penalty function is designed to deal with the constraints of fuel and number of TDKs, which is constructed as(10)GX,S=∑i=1n−1φMfi−ΔMii+1+∑i=1n−1φkleft,i−1,where the penalty function φ⋅ is defined as(11)φx=ax2−bx+c,x<0,0,x≥0,a,b,c>0The penalty function is adopted as objective for optimization to solve the problem of satellite layout optimization designed by Chen [25]. a,b,c are the given constants. In this paper, a,b,c are 1, 2, 0 respectively. Mfi−ΔMii+1<0 or kleft,i<0 means that the constraint is not satisfied, and the penalty function GX,S will add to the objective function.


#### 3.1.3. Two-Level Optimization Model

Method for the ADR problem can be divided to the one-level and two-level optimization. The one-level optimization means that all designed variables are optimized directly. The two-level optimization includes the up-level to optimize the mission sequence and the low-level to find the optimal transfer orbit. As shown in Figure2, low-level optimization is nested in the up-level one; in other words, the up-level objective function is calculated by optimizing the low-level one. The up-level planning optimizes the sequence, while the low-level planning optimizes the object-to-object transfer impulse. Zhou [35] has demonstrated that the results of two-level optimization are better than the ones of one-level optimization. The reason is that the solution space of the up-level optimization is much smaller than that of the one-level optimization. In low-level optimization, APSO is used, which integrates a local acceleration mechanism based on gradient search in iterations, to solve optimal transfer impulse, and it can get faster calculation and maintain global search capability. In up-level optimization, IGA is an improved algorithm that inherits the genetic algorithm and integrates immune algorithm into it. IGA has the ability to generate diverse antibodies, self-regulating mechanisms, and immune memory. These mechanisms can get global search capability without falling into the local optimum. Through antigen memory, the search speed can be accelerated, and search capability of the genetic algorithm can be improved. In this way, IGA is more suitable for sequence optimization search problem of ADR.Figure 2
Two-level optimization.The two-level optimization model can be given as(12)findX=x1,x2,…,xnS=s1,s2,…,snminΔMs.t.xi≤nxi∈1,2,…,nxi≠xji,j∈1,2,…,n,i≠jsi∈0,1i∈1,2,…,nfindT=Δt1,Δt2,…Δtm,tw1,tw2,…,twmminΔMs.t.Δti+twi≤tmaxMfi−ΔMii+1≥0,i=1,2,…,n−1kleft,i≥0.


### 3.2. Optimizer Algorithm Design

#### 3.2.1. Low-Level Optimization

The low-level optimization is to determine the optimal transfer time for service satellite, the variable of which isT=Δt1,Δt2,…Δtm,tw1,tw2,…,twm. The objective for optimization is to minimize fuel consumption, shown in equation (5).In this article, the novel accelerated particle swarm optimization (APSO) algorithm is adopted in the low-level optimization [25]. It aims to fully utilize the global exploration of PSO and the rapid local exploitation of the gradient-based method to speed up the global search process. The traditional PSO algorithm has problems such as slow convergence speed, low computational efficiency, and inability to maintain the global search ability in high-dimensional problems. This paper uses PSO as basic framework and integrates a gradient-based search in its iterations, and in this way, APSO algorithm is proposed. When the number of iterations reaches a certain number Tacc, the current position of each particle is used as an initial point for the sequential quadratic programming (SQP) algorithm to exploit its local best and update its position until all the particles in the swarm have been accelerated. The process of APSO is shown in the right part of Figure 3. In this way, APSO can get faster calculation speed when optimizing object-to-object impulse.Figure 3
Flowchart of two-level optimization.


#### 3.2.2. Up-Level Optimization

The immune genetic algorithm (IGA) is adopted for up-level optimization. Compared with traditional GA, IGA is an improved genetic algorithm based on biological immune mechanism. The objective function of the problem corresponds to the antigen of the invading organism, and the solution corresponds to antibody produced by the immune system. The algorithm has the following features: (1) diverse antibodies: in IGA algorithm, the ability of the immune system to generate a large number of antibodies to resist various antigens is simulated. This mechanism is used to improve the global search ability of the genetic algorithm rather than falling into the local optimum. (2) Self-regulation mechanism: the immune system has a mechanism to maintain immune balance. By inhibiting and promoting antibodies, it can self-regulate to produce an appropriate amount of antibodies. Using this function, it can improve local search ability of genetic algorithm. (3) Immune memory cell function: some of the cells that produce affinity will be preserved as memory cells. For similar antigens that will invade in the future, the corresponding memory cells will be quickly stimulated to produce a large amount of antibodies. Using this memory function, it can speed up search speed and improve the overall search ability of genetic algorithm.In IGA, antigen stands for objective function, and antibody represents designed variables. The affinity between the antigen and the antibody indicates the matching degree between objective function and optimal variables. Each generation selects excellent antibodies by their affinity as memory cells. The affinity can be calculated as follows:(13)exc=1fit∗p+con∗1−p,0<p<1,where fit and con represent the objective function and antibody concentration, respectively, and p is the evaluation parameter. From equation (13), the calculation for affinity is not only based on the objective function, but also based on the antibody concentration. The high concentration will cause a decrease in affinity. This operation will keep individual varieties and avoid premature convergence.Code, select, crossover, and mutation are used in IGA. Each antibody is coded asC=X,S. The roulette wheel method is used in selection [36]. The one-point method is utilized in crossover. Select one position k in antibody and two different antibodies randomly. The mutation means, select two genes in variable X and exchange them, and then select one gene in variable S and change it to another value. The process of IGA is shown in the left part of Figure 3.Table 1
Parameters for APSO algorithm.
ParametersValueDescriptionsN30Particle sizeTmax100Max number of iterationTacc70The iteration of performing particle accelerationwmax0.8Max inertia weightwmin0.4Min inertia weightc11.5Cognitive factor 1c21.5Cognitive factor 2Table 2
Parameters for IGA algorithm.
ItemvalueDescriptionsNpop50Pop sizeNmem20Memory cells numberTmax100Max generationsPcro0.7Rate of crossoverPmut0.1Rate of mutationp0.95Evaluation parameterIn summary, the procedure of two-level optimization is summarized as follows:Step 1.
:set parameters (shown in Tables1 and 2) to initialize APSO and IGAStep 2.
:calculate the fitness function using APSO in low-level optimizationStep 3.
:calculate the affinity based on the fitness value and antibody concentrationStep 4.
:select the memory cells and operate the process of selection, crossover, and mutation in IGA of up-level optimizationStep 5.
:repeat the Steps2–5 until reaching maximum iterationsThe whole process of two-level optimization is shown in Figure3.


## 3.1. Optimization Model

### 3.1.1. Design Variables

For mission planning problem, it is necessary to determine the ADR sequence [34]. The design variables consist of removal order X, binary decision variable S, and time distribution T. The X and S are integers, and T is a continuous number, respectively. The variables are of the form(6)X=x1,x2,…,xn,S=s1,s2,…,sn,T=Δt1,Δt2,…Δtm,tw1,tw2,…,twm,For example, the removing orderX=2,5,4,3,1 denotes that 5 pieces of debris need to be cleaned, and service satellite removes the debris as the order 2 ⟶ 5 ⟶ 4 ⟶ 3 ⟶ 1. S=0,0,1,0,1means that service satellite transfers to space depot for replenishment before removing the 3-th and 5-th debris. In T, Δt1,Δt2,…Δtmindicate the time period for each maneuver, and tw1,tw2,…,twmare the waiting time before orbit transfer.In the ADR mission, the objective is to minimize the total fuel consumption,MinΔM.


### 3.1.2. Constraints

(1) Fuel Constraint. The service satellite cannot maneuver if the remaining fuel is insufficient, and the removing task will be interrupted. The constraint for fuel can be stated as(7)Mfi−ΔMi≥0,i=1,2,…,k−1,where Mfiis remaining fuel before removing the i-th debris.(2) TDK Number Constraint. TDK is being consumed during the mission. Therefore, the number of TDKs is required to be sufficient and can be stated as(8)kleft,i≥0,where kleft,i is the number of left TDKs the service satellite carries before removing the i-th debris.(3) Time Distribution Constraint. The ADR mission usually requires timeliness. The constraint is(9)Δti+twi≤tmax.The integrated penalty function is designed to deal with the constraints of fuel and number of TDKs, which is constructed as(10)GX,S=∑i=1n−1φMfi−ΔMii+1+∑i=1n−1φkleft,i−1,where the penalty function φ⋅ is defined as(11)φx=ax2−bx+c,x<0,0,x≥0,a,b,c>0The penalty function is adopted as objective for optimization to solve the problem of satellite layout optimization designed by Chen [25]. a,b,c are the given constants. In this paper, a,b,c are 1, 2, 0 respectively. Mfi−ΔMii+1<0 or kleft,i<0 means that the constraint is not satisfied, and the penalty function GX,S will add to the objective function.


### 3.1.3. Two-Level Optimization Model

Method for the ADR problem can be divided to the one-level and two-level optimization. The one-level optimization means that all designed variables are optimized directly. The two-level optimization includes the up-level to optimize the mission sequence and the low-level to find the optimal transfer orbit. As shown in Figure2, low-level optimization is nested in the up-level one; in other words, the up-level objective function is calculated by optimizing the low-level one. The up-level planning optimizes the sequence, while the low-level planning optimizes the object-to-object transfer impulse. Zhou [35] has demonstrated that the results of two-level optimization are better than the ones of one-level optimization. The reason is that the solution space of the up-level optimization is much smaller than that of the one-level optimization. In low-level optimization, APSO is used, which integrates a local acceleration mechanism based on gradient search in iterations, to solve optimal transfer impulse, and it can get faster calculation and maintain global search capability. In up-level optimization, IGA is an improved algorithm that inherits the genetic algorithm and integrates immune algorithm into it. IGA has the ability to generate diverse antibodies, self-regulating mechanisms, and immune memory. These mechanisms can get global search capability without falling into the local optimum. Through antigen memory, the search speed can be accelerated, and search capability of the genetic algorithm can be improved. In this way, IGA is more suitable for sequence optimization search problem of ADR.Figure 2
Two-level optimization.The two-level optimization model can be given as(12)findX=x1,x2,…,xnS=s1,s2,…,snminΔMs.t.xi≤nxi∈1,2,…,nxi≠xji,j∈1,2,…,n,i≠jsi∈0,1i∈1,2,…,nfindT=Δt1,Δt2,…Δtm,tw1,tw2,…,twmminΔMs.t.Δti+twi≤tmaxMfi−ΔMii+1≥0,i=1,2,…,n−1kleft,i≥0.


## 3.1.1. Design Variables

For mission planning problem, it is necessary to determine the ADR sequence [34]. The design variables consist of removal order X, binary decision variable S, and time distribution T. The X and S are integers, and T is a continuous number, respectively. The variables are of the form(6)X=x1,x2,…,xn,S=s1,s2,…,sn,T=Δt1,Δt2,…Δtm,tw1,tw2,…,twm,For example, the removing orderX=2,5,4,3,1 denotes that 5 pieces of debris need to be cleaned, and service satellite removes the debris as the order 2 ⟶ 5 ⟶ 4 ⟶ 3 ⟶ 1. S=0,0,1,0,1means that service satellite transfers to space depot for replenishment before removing the 3-th and 5-th debris. In T, Δt1,Δt2,…Δtmindicate the time period for each maneuver, and tw1,tw2,…,twmare the waiting time before orbit transfer.In the ADR mission, the objective is to minimize the total fuel consumption,MinΔM.


## 3.1.2. Constraints

(1) Fuel Constraint. The service satellite cannot maneuver if the remaining fuel is insufficient, and the removing task will be interrupted. The constraint for fuel can be stated as(7)Mfi−ΔMi≥0,i=1,2,…,k−1,where Mfiis remaining fuel before removing the i-th debris.(2) TDK Number Constraint. TDK is being consumed during the mission. Therefore, the number of TDKs is required to be sufficient and can be stated as(8)kleft,i≥0,where kleft,i is the number of left TDKs the service satellite carries before removing the i-th debris.(3) Time Distribution Constraint. The ADR mission usually requires timeliness. The constraint is(9)Δti+twi≤tmax.The integrated penalty function is designed to deal with the constraints of fuel and number of TDKs, which is constructed as(10)GX,S=∑i=1n−1φMfi−ΔMii+1+∑i=1n−1φkleft,i−1,where the penalty function φ⋅ is defined as(11)φx=ax2−bx+c,x<0,0,x≥0,a,b,c>0The penalty function is adopted as objective for optimization to solve the problem of satellite layout optimization designed by Chen [25]. a,b,c are the given constants. In this paper, a,b,c are 1, 2, 0 respectively. Mfi−ΔMii+1<0 or kleft,i<0 means that the constraint is not satisfied, and the penalty function GX,S will add to the objective function.


## 3.1.3. Two-Level Optimization Model

Method for the ADR problem can be divided to the one-level and two-level optimization. The one-level optimization means that all designed variables are optimized directly. The two-level optimization includes the up-level to optimize the mission sequence and the low-level to find the optimal transfer orbit. As shown in Figure2, low-level optimization is nested in the up-level one; in other words, the up-level objective function is calculated by optimizing the low-level one. The up-level planning optimizes the sequence, while the low-level planning optimizes the object-to-object transfer impulse. Zhou [35] has demonstrated that the results of two-level optimization are better than the ones of one-level optimization. The reason is that the solution space of the up-level optimization is much smaller than that of the one-level optimization. In low-level optimization, APSO is used, which integrates a local acceleration mechanism based on gradient search in iterations, to solve optimal transfer impulse, and it can get faster calculation and maintain global search capability. In up-level optimization, IGA is an improved algorithm that inherits the genetic algorithm and integrates immune algorithm into it. IGA has the ability to generate diverse antibodies, self-regulating mechanisms, and immune memory. These mechanisms can get global search capability without falling into the local optimum. Through antigen memory, the search speed can be accelerated, and search capability of the genetic algorithm can be improved. In this way, IGA is more suitable for sequence optimization search problem of ADR.Figure 2
Two-level optimization.The two-level optimization model can be given as(12)findX=x1,x2,…,xnS=s1,s2,…,snminΔMs.t.xi≤nxi∈1,2,…,nxi≠xji,j∈1,2,…,n,i≠jsi∈0,1i∈1,2,…,nfindT=Δt1,Δt2,…Δtm,tw1,tw2,…,twmminΔMs.t.Δti+twi≤tmaxMfi−ΔMii+1≥0,i=1,2,…,n−1kleft,i≥0.


## 3.2. Optimizer Algorithm Design

### 3.2.1. Low-Level Optimization

The low-level optimization is to determine the optimal transfer time for service satellite, the variable of which isT=Δt1,Δt2,…Δtm,tw1,tw2,…,twm. The objective for optimization is to minimize fuel consumption, shown in equation (5).In this article, the novel accelerated particle swarm optimization (APSO) algorithm is adopted in the low-level optimization [25]. It aims to fully utilize the global exploration of PSO and the rapid local exploitation of the gradient-based method to speed up the global search process. The traditional PSO algorithm has problems such as slow convergence speed, low computational efficiency, and inability to maintain the global search ability in high-dimensional problems. This paper uses PSO as basic framework and integrates a gradient-based search in its iterations, and in this way, APSO algorithm is proposed. When the number of iterations reaches a certain number Tacc, the current position of each particle is used as an initial point for the sequential quadratic programming (SQP) algorithm to exploit its local best and update its position until all the particles in the swarm have been accelerated. The process of APSO is shown in the right part of Figure 3. In this way, APSO can get faster calculation speed when optimizing object-to-object impulse.Figure 3
Flowchart of two-level optimization.


### 3.2.2. Up-Level Optimization

The immune genetic algorithm (IGA) is adopted for up-level optimization. Compared with traditional GA, IGA is an improved genetic algorithm based on biological immune mechanism. The objective function of the problem corresponds to the antigen of the invading organism, and the solution corresponds to antibody produced by the immune system. The algorithm has the following features: (1) diverse antibodies: in IGA algorithm, the ability of the immune system to generate a large number of antibodies to resist various antigens is simulated. This mechanism is used to improve the global search ability of the genetic algorithm rather than falling into the local optimum. (2) Self-regulation mechanism: the immune system has a mechanism to maintain immune balance. By inhibiting and promoting antibodies, it can self-regulate to produce an appropriate amount of antibodies. Using this function, it can improve local search ability of genetic algorithm. (3) Immune memory cell function: some of the cells that produce affinity will be preserved as memory cells. For similar antigens that will invade in the future, the corresponding memory cells will be quickly stimulated to produce a large amount of antibodies. Using this memory function, it can speed up search speed and improve the overall search ability of genetic algorithm.In IGA, antigen stands for objective function, and antibody represents designed variables. The affinity between the antigen and the antibody indicates the matching degree between objective function and optimal variables. Each generation selects excellent antibodies by their affinity as memory cells. The affinity can be calculated as follows:(13)exc=1fit∗p+con∗1−p,0<p<1,where fit and con represent the objective function and antibody concentration, respectively, and p is the evaluation parameter. From equation (13), the calculation for affinity is not only based on the objective function, but also based on the antibody concentration. The high concentration will cause a decrease in affinity. This operation will keep individual varieties and avoid premature convergence.Code, select, crossover, and mutation are used in IGA. Each antibody is coded asC=X,S. The roulette wheel method is used in selection [36]. The one-point method is utilized in crossover. Select one position k in antibody and two different antibodies randomly. The mutation means, select two genes in variable X and exchange them, and then select one gene in variable S and change it to another value. The process of IGA is shown in the left part of Figure 3.Table 1
Parameters for APSO algorithm.
ParametersValueDescriptionsN30Particle sizeTmax100Max number of iterationTacc70The iteration of performing particle accelerationwmax0.8Max inertia weightwmin0.4Min inertia weightc11.5Cognitive factor 1c21.5Cognitive factor 2Table 2
Parameters for IGA algorithm.
ItemvalueDescriptionsNpop50Pop sizeNmem20Memory cells numberTmax100Max generationsPcro0.7Rate of crossoverPmut0.1Rate of mutationp0.95Evaluation parameterIn summary, the procedure of two-level optimization is summarized as follows:Step 1.
:set parameters (shown in Tables1 and 2) to initialize APSO and IGAStep 2.
:calculate the fitness function using APSO in low-level optimizationStep 3.
:calculate the affinity based on the fitness value and antibody concentrationStep 4.
:select the memory cells and operate the process of selection, crossover, and mutation in IGA of up-level optimizationStep 5.
:repeat the Steps2–5 until reaching maximum iterationsThe whole process of two-level optimization is shown in Figure3.


## 3.2.1. Low-Level Optimization

The low-level optimization is to determine the optimal transfer time for service satellite, the variable of which isT=Δt1,Δt2,…Δtm,tw1,tw2,…,twm. The objective for optimization is to minimize fuel consumption, shown in equation (5).In this article, the novel accelerated particle swarm optimization (APSO) algorithm is adopted in the low-level optimization [25]. It aims to fully utilize the global exploration of PSO and the rapid local exploitation of the gradient-based method to speed up the global search process. The traditional PSO algorithm has problems such as slow convergence speed, low computational efficiency, and inability to maintain the global search ability in high-dimensional problems. This paper uses PSO as basic framework and integrates a gradient-based search in its iterations, and in this way, APSO algorithm is proposed. When the number of iterations reaches a certain number Tacc, the current position of each particle is used as an initial point for the sequential quadratic programming (SQP) algorithm to exploit its local best and update its position until all the particles in the swarm have been accelerated. The process of APSO is shown in the right part of Figure 3. In this way, APSO can get faster calculation speed when optimizing object-to-object impulse.Figure 3
Flowchart of two-level optimization.


## 3.2.2. Up-Level Optimization

The immune genetic algorithm (IGA) is adopted for up-level optimization. Compared with traditional GA, IGA is an improved genetic algorithm based on biological immune mechanism. The objective function of the problem corresponds to the antigen of the invading organism, and the solution corresponds to antibody produced by the immune system. The algorithm has the following features: (1) diverse antibodies: in IGA algorithm, the ability of the immune system to generate a large number of antibodies to resist various antigens is simulated. This mechanism is used to improve the global search ability of the genetic algorithm rather than falling into the local optimum. (2) Self-regulation mechanism: the immune system has a mechanism to maintain immune balance. By inhibiting and promoting antibodies, it can self-regulate to produce an appropriate amount of antibodies. Using this function, it can improve local search ability of genetic algorithm. (3) Immune memory cell function: some of the cells that produce affinity will be preserved as memory cells. For similar antigens that will invade in the future, the corresponding memory cells will be quickly stimulated to produce a large amount of antibodies. Using this memory function, it can speed up search speed and improve the overall search ability of genetic algorithm.In IGA, antigen stands for objective function, and antibody represents designed variables. The affinity between the antigen and the antibody indicates the matching degree between objective function and optimal variables. Each generation selects excellent antibodies by their affinity as memory cells. The affinity can be calculated as follows:(13)exc=1fit∗p+con∗1−p,0<p<1,where fit and con represent the objective function and antibody concentration, respectively, and p is the evaluation parameter. From equation (13), the calculation for affinity is not only based on the objective function, but also based on the antibody concentration. The high concentration will cause a decrease in affinity. This operation will keep individual varieties and avoid premature convergence.Code, select, crossover, and mutation are used in IGA. Each antibody is coded asC=X,S. The roulette wheel method is used in selection [36]. The one-point method is utilized in crossover. Select one position k in antibody and two different antibodies randomly. The mutation means, select two genes in variable X and exchange them, and then select one gene in variable S and change it to another value. The process of IGA is shown in the left part of Figure 3.Table 1
Parameters for APSO algorithm.
ParametersValueDescriptionsN30Particle sizeTmax100Max number of iterationTacc70The iteration of performing particle accelerationwmax0.8Max inertia weightwmin0.4Min inertia weightc11.5Cognitive factor 1c21.5Cognitive factor 2Table 2
Parameters for IGA algorithm.
ItemvalueDescriptionsNpop50Pop sizeNmem20Memory cells numberTmax100Max generationsPcro0.7Rate of crossoverPmut0.1Rate of mutationp0.95Evaluation parameterIn summary, the procedure of two-level optimization is summarized as follows:Step 1.
:set parameters (shown in Tables1 and 2) to initialize APSO and IGAStep 2.
:calculate the fitness function using APSO in low-level optimizationStep 3.
:calculate the affinity based on the fitness value and antibody concentrationStep 4.
:select the memory cells and operate the process of selection, crossover, and mutation in IGA of up-level optimizationStep 5.
:repeat the Steps2–5 until reaching maximum iterationsThe whole process of two-level optimization is shown in Figure3.


## 4. Simulation and Analysis

Select 9 pieces of space debris on near GEO orbit. The orbital elements of debris are shown in Table3. To avoid the singularity in classical orbital elements, the eccentricity is initialized as 0.001 as near-circle orbit. The semimajor axis and argument of perigee are initialized as 42166.3 km for GEO and 30°, respectively. It is assumed that both the space depot and service satellite are initially located in the near GEO orbit, where the semimajor axis, eccentricity, inclination, right ascension of ascending node (RAAN), argument of latitude, and true anomaly are 42166.3 km, 0.001, 0.001°, 0°, 0°, and 0°, respectively.Table 3
Initial orbital elements of debris.
Serial numberInclination (°)RAAN (°)True anomaly (°)SD1511015SD2728046SD36.512594SD445138SD5570162SD66.895217SD74145239SD8815277SD94.755334The constants used in this simulation areIsp=300s and g=9.8. Time period for ADR mission is tj=1000s. The refueling time is tf=8000s. The maximum time tmax for each maneuver is 5 orbital periods.The parameters of APSO and IGA are provided in Tables1 and 2.


### 4.1. Case 1

In this case, the efficiency of the two-level optimization is tested by comparing with exhaustive search (ES) method, that is, the enumeration of possible solutions and the choice of the best one. The consumption time for ES grows exponentially with the instance size, so the serial numbers 1–6 of objects in Table2 are selected to be removed. In that way, there are 26⋅6!=46080 different sequences in total, which can be the search space of simulation. The mass of propellant and TDK is 400 kg and 50 kg, respectively, the consumption of which meets the result of the argument in reference [16]. Set the structural mass of service satellite as 400 kg. Each satellite carries 4 TDKs. Therefore, the total mass of service satellite is 1000 kg. The fuel consumption for all sequence is shown in Figure 4. The extreme data with fuel consumption, which is greater than 2000 kg, is removed from figure. Each point in the figure stands for one possible sequence. According to the frequency of replenishment, the mission sequence can be divided into 7 categories, which, respectively, means that the service satellite is replenished 0–6 times. The minimum fuel cost is 443.22 kg, as the optimal mission sequence isX=1,3,6,5,4,2, S=0,0,0,0,1,0. It takes about 49827.5 s for a 3.6 GHZ Intel Core personal computer with a memory of 16 GB to complete this calculation.Figure 4
Fuel consumption of all mission sequence.Using the proposed two-level optimization method, the convergence curve is shown in Figure5. The algorithm converges after about 25 generations. The optimal mission sequence is X=1,3,6,5,4,2, S=0,0,0,0,1,0and the minimum fuel cost is 443.22 kg, which is consistent with result of ES method. However, in two-level optimization, the calculation time in the same computer is 4326.9 s. Clearly, these two methods (ES and two-level optimization method) can both find the optimal solution, but the time consumed by two-level optimization is reduced by 91.3%.Figure 5
Evolutionary curve of two-level optimization.From Figure4, the fuel cost of each category is concentrated in a certain range. For example, 84.65 percent (6616 of 7815 in total) of the data with twice replenishment distribute in the range of fuel consumption 600–1200 kg. Data analysis of the fuel consumption with different sequence category is shown in Figure 6. The data of zero replenishing frequency is removed, because the TDK constraint is not satisfied. In general, when the initial mass is the same, higher replenishing frequency means larger fuel consumption.Figure 6
Fuel consumption with different replenishing frequency.


### 4.2. Case 2

In the revolver mode for debris removing, the quantity of TDKs is an important part for the whole mission. The different TDK number constraints in mission planning are discussed in case 2. Keep structural mass and TDK mass the same as those in case 1. All nine objects in Table3 are selected to be removed. The number of TDKs carried by service satellite is 2, 3, 4, 5, 6, 7, 8, and 9, respectively, and the total mass changes to 900 kg, 950 kg, 1000 kg, 1050 kg, 1100 kg, 1150 kg, 1200 kg, and 1250 kg with the change of the TDKs number.Table4 addresses the optimized mission sequence with different number of TDKs. It is supposed that each debris should be removed by one attached TDK. The minimum fuel consumption for different TDK numbers is 1014.1 kg, 854.4 kg, 904.9 kg, 726.2 kg, 934.1 kg, 981.7 kg, 986.8 kg, and 1066.4 kg, respectively. As the quantity of total pieces of space debris is 9, for the TDK number with 2, 3, and 4, service satellite maneuvers to space depot more than once for TDK replenishment, because the TDKs the service satellite carries are not enough for removing all pieces of debris. This situation increases the fuel consumption as shown in case 1, in which higher frequency of replenishment means larger fuel consumption. For the number of 5, service satellite needs to maneuver to depot once, and the fuel cost is small. When the TDK number increases to 6, 7, 8, and 9, the service satellite maneuvers to space depot twice due to the lack of propellant, and the fuel consumption is large. It can be analyzed that the TDK number influences the replenishing frequency through the constraints of TDKs number and fuel. Fewer TDKs the service satellite carries will cause high replenishing frequency due to constraint of TDK, and larger quantity will have high replenishing frequency as well, due to the constraint of fuel.Table 4
Optimized results with different TDK number.
TDK numberRemoving order XBinary decision variable SReplenishing frequencyTDK redundancyFuel consumption (kg)23-6-4-8-7-1-5-9-20, 0, 1, 0, 1, 0, 1, 0, 1411014.131-3-6-4-5-9-2-7-80, 0, 0, 1, 0, 0, 1, 0, 020854.444-8-9-3-6-5-2-7-10, 0, 0, 1, 0, 0, 1, 0, 023904.952-7-1-3-6-4-8-9-50, 0, 0, 0, 0, 1, 0, 0, 011726.267-1-3-6-4-5-9-2-80, 0, 0, 0, 1, 0, 0, 1, 029934.177-1-3-6-4-5-9-2-80, 0, 0, 0, 1, 0, 0, 1, 0212981.788-9-5-4-2-7-1-3-60, 0, 0, 1, 0, 1, 0, 0, 0215986.898-9-6-3-2-4-5-7-10, 0, 0, 1, 0, 0, 1, 0, 02181066.4In this case, a new concept is introduced: TDK redundancy. The TDK redundancy indicates the total remaining quantity of TDK before replenishment and after the whole mission. This index reflects the usage rate of TDK. Obviously, the initial TDK number and replenishing frequency determine the TDK redundancy. As shown in Table4, for the TDK numbers of 7, 8, and 9, when the replenishing frequency is the same, higher TDK number causes larger TDK redundancy and further leads to larger fuel cost. The reason is that, during the orbital transfer, service satellite will have higher mass (redundant TDKs).To further analyze the influence, the TDK consumption is initialized as1,1,2,1,2,1,1,1,1,. This indicates that the 2 TDKs are consumed for removing debris numbered #3 and #3. The number of TDK carried by service satellite is 2, 3, 4, 5, 6, 7, 8, 9, 10, and 11, respectively, and the total mass changes to 900 kg, 950 kg, 1000 kg, 1050 kg, 1100 kg, 1150 kg, 1200 kg, 1250 kg, 1300 kg, and 1350 kg.Table5 addresses the optimized mission sequence with different number of TDK. The minimum fuel consumption is 861.4 kg with 4 TDK number, twice replenishment, and 1 TDK redundancy. And the suboptimal fuel consumption is 904 kg with 6 TDK number, one replenishment, and 1 TDK redundancy. This indicates that the lower replenishing frequency does not always mean fewer fuel cost, as the mass may be large in orbit transfer. For the TDK numbers as 2 and 3, the replenishing frequency is high, and the fuel consumption is large. For the TDK numbers as 4, 5, 7, 8, 9, 10, and 11, the replenishing frequency is the same, and the larger TDK redundancy leads to more fuel consumption.Table 5
Optimized results with different TDK number and TDK consumption.
TDK numberRemoving order XBinary decision variable SReplenishing FrequencyTDK redundancyFuel consumption (kg)27-1-4-8-9-5-2-3-60, 0, 1, 0, 1, 1, 1, 1, 1531165.135-3-6-4-8-9-2-7-10, 1, 0, 1, 0, 0, 1, 0, 031965.948-9-4-2-5-6-7-1-30, 0, 0, 0, 1, 0, 1, 0, 021861.457-1-2-4-8-9-5-6-30, 0, 0, 1, 0, 0, 1, 0, 024915.663-6-9-5-4-2-7-1-80, 0, 0, 0, 1, 0, 0, 0, 011904.778-9-5-4-2-7-1-3-60, 0, 0, 1, 0, 1, 0, 0, 0210931.383-6-9-5-2-7-1-4-80, 0, 0, 0, 0, 1, 0, 1, 0213971.297-1-3-6-5-2-4-8-90, 0, 0, 0, 1, 0, 1, 0, 0216987.8107-1-6-3-4-5-9-2-80, 0, 0, 0, 1, 0, 0, 1, 02191118.6118-7-1-4-2-3-6-9-50, 0, 0, 1, 0, 1, 0, 0, 02221506.1To sum up, TDK is the main constraint of mission. The TDK number influences replenishing frequency and TDK redundancy and then impacts the fuel consumption. In general, high replenishing frequency means large fuel consumption. But in special conditions, the lowest replenishing frequency may not mean minimum fuel cost, because more TDKs carried by service satellite will lower the frequency of replenishment, while, in turn, they may cause higher mass in orbit transfer. When the replenishing frequency is the same, high TDK redundancy causes more fuel cost. The two indexes (replenishing frequency and TDK redundancy) almost determine the fuel cost in case the target debris is specific. Therefore, for the revolver mode, the TDK number the service satellite carries should be designed beforehand with suitable replenishing frequency and minimum TDK redundancy. This is, in fact, take full advantage of TDK carried by service satellite and reduce orbital maneuver with small satellite mass.


## 4.1. Case 1

In this case, the efficiency of the two-level optimization is tested by comparing with exhaustive search (ES) method, that is, the enumeration of possible solutions and the choice of the best one. The consumption time for ES grows exponentially with the instance size, so the serial numbers 1–6 of objects in Table2 are selected to be removed. In that way, there are 26⋅6!=46080 different sequences in total, which can be the search space of simulation. The mass of propellant and TDK is 400 kg and 50 kg, respectively, the consumption of which meets the result of the argument in reference [16]. Set the structural mass of service satellite as 400 kg. Each satellite carries 4 TDKs. Therefore, the total mass of service satellite is 1000 kg. The fuel consumption for all sequence is shown in Figure 4. The extreme data with fuel consumption, which is greater than 2000 kg, is removed from figure. Each point in the figure stands for one possible sequence. According to the frequency of replenishment, the mission sequence can be divided into 7 categories, which, respectively, means that the service satellite is replenished 0–6 times. The minimum fuel cost is 443.22 kg, as the optimal mission sequence isX=1,3,6,5,4,2, S=0,0,0,0,1,0. It takes about 49827.5 s for a 3.6 GHZ Intel Core personal computer with a memory of 16 GB to complete this calculation.Figure 4
Fuel consumption of all mission sequence.Using the proposed two-level optimization method, the convergence curve is shown in Figure5. The algorithm converges after about 25 generations. The optimal mission sequence is X=1,3,6,5,4,2, S=0,0,0,0,1,0and the minimum fuel cost is 443.22 kg, which is consistent with result of ES method. However, in two-level optimization, the calculation time in the same computer is 4326.9 s. Clearly, these two methods (ES and two-level optimization method) can both find the optimal solution, but the time consumed by two-level optimization is reduced by 91.3%.Figure 5
Evolutionary curve of two-level optimization.From Figure4, the fuel cost of each category is concentrated in a certain range. For example, 84.65 percent (6616 of 7815 in total) of the data with twice replenishment distribute in the range of fuel consumption 600–1200 kg. Data analysis of the fuel consumption with different sequence category is shown in Figure 6. The data of zero replenishing frequency is removed, because the TDK constraint is not satisfied. In general, when the initial mass is the same, higher replenishing frequency means larger fuel consumption.Figure 6
Fuel consumption with different replenishing frequency.


## 4.2. Case 2

In the revolver mode for debris removing, the quantity of TDKs is an important part for the whole mission. The different TDK number constraints in mission planning are discussed in case 2. Keep structural mass and TDK mass the same as those in case 1. All nine objects in Table3 are selected to be removed. The number of TDKs carried by service satellite is 2, 3, 4, 5, 6, 7, 8, and 9, respectively, and the total mass changes to 900 kg, 950 kg, 1000 kg, 1050 kg, 1100 kg, 1150 kg, 1200 kg, and 1250 kg with the change of the TDKs number.Table4 addresses the optimized mission sequence with different number of TDKs. It is supposed that each debris should be removed by one attached TDK. The minimum fuel consumption for different TDK numbers is 1014.1 kg, 854.4 kg, 904.9 kg, 726.2 kg, 934.1 kg, 981.7 kg, 986.8 kg, and 1066.4 kg, respectively. As the quantity of total pieces of space debris is 9, for the TDK number with 2, 3, and 4, service satellite maneuvers to space depot more than once for TDK replenishment, because the TDKs the service satellite carries are not enough for removing all pieces of debris. This situation increases the fuel consumption as shown in case 1, in which higher frequency of replenishment means larger fuel consumption. For the number of 5, service satellite needs to maneuver to depot once, and the fuel cost is small. When the TDK number increases to 6, 7, 8, and 9, the service satellite maneuvers to space depot twice due to the lack of propellant, and the fuel consumption is large. It can be analyzed that the TDK number influences the replenishing frequency through the constraints of TDKs number and fuel. Fewer TDKs the service satellite carries will cause high replenishing frequency due to constraint of TDK, and larger quantity will have high replenishing frequency as well, due to the constraint of fuel.Table 4
Optimized results with different TDK number.
TDK numberRemoving order XBinary decision variable SReplenishing frequencyTDK redundancyFuel consumption (kg)23-6-4-8-7-1-5-9-20, 0, 1, 0, 1, 0, 1, 0, 1411014.131-3-6-4-5-9-2-7-80, 0, 0, 1, 0, 0, 1, 0, 020854.444-8-9-3-6-5-2-7-10, 0, 0, 1, 0, 0, 1, 0, 023904.952-7-1-3-6-4-8-9-50, 0, 0, 0, 0, 1, 0, 0, 011726.267-1-3-6-4-5-9-2-80, 0, 0, 0, 1, 0, 0, 1, 029934.177-1-3-6-4-5-9-2-80, 0, 0, 0, 1, 0, 0, 1, 0212981.788-9-5-4-2-7-1-3-60, 0, 0, 1, 0, 1, 0, 0, 0215986.898-9-6-3-2-4-5-7-10, 0, 0, 1, 0, 0, 1, 0, 02181066.4In this case, a new concept is introduced: TDK redundancy. The TDK redundancy indicates the total remaining quantity of TDK before replenishment and after the whole mission. This index reflects the usage rate of TDK. Obviously, the initial TDK number and replenishing frequency determine the TDK redundancy. As shown in Table4, for the TDK numbers of 7, 8, and 9, when the replenishing frequency is the same, higher TDK number causes larger TDK redundancy and further leads to larger fuel cost. The reason is that, during the orbital transfer, service satellite will have higher mass (redundant TDKs).To further analyze the influence, the TDK consumption is initialized as1,1,2,1,2,1,1,1,1,. This indicates that the 2 TDKs are consumed for removing debris numbered #3 and #3. The number of TDK carried by service satellite is 2, 3, 4, 5, 6, 7, 8, 9, 10, and 11, respectively, and the total mass changes to 900 kg, 950 kg, 1000 kg, 1050 kg, 1100 kg, 1150 kg, 1200 kg, 1250 kg, 1300 kg, and 1350 kg.Table5 addresses the optimized mission sequence with different number of TDK. The minimum fuel consumption is 861.4 kg with 4 TDK number, twice replenishment, and 1 TDK redundancy. And the suboptimal fuel consumption is 904 kg with 6 TDK number, one replenishment, and 1 TDK redundancy. This indicates that the lower replenishing frequency does not always mean fewer fuel cost, as the mass may be large in orbit transfer. For the TDK numbers as 2 and 3, the replenishing frequency is high, and the fuel consumption is large. For the TDK numbers as 4, 5, 7, 8, 9, 10, and 11, the replenishing frequency is the same, and the larger TDK redundancy leads to more fuel consumption.Table 5
Optimized results with different TDK number and TDK consumption.
TDK numberRemoving order XBinary decision variable SReplenishing FrequencyTDK redundancyFuel consumption (kg)27-1-4-8-9-5-2-3-60, 0, 1, 0, 1, 1, 1, 1, 1531165.135-3-6-4-8-9-2-7-10, 1, 0, 1, 0, 0, 1, 0, 031965.948-9-4-2-5-6-7-1-30, 0, 0, 0, 1, 0, 1, 0, 021861.457-1-2-4-8-9-5-6-30, 0, 0, 1, 0, 0, 1, 0, 024915.663-6-9-5-4-2-7-1-80, 0, 0, 0, 1, 0, 0, 0, 011904.778-9-5-4-2-7-1-3-60, 0, 0, 1, 0, 1, 0, 0, 0210931.383-6-9-5-2-7-1-4-80, 0, 0, 0, 0, 1, 0, 1, 0213971.297-1-3-6-5-2-4-8-90, 0, 0, 0, 1, 0, 1, 0, 0216987.8107-1-6-3-4-5-9-2-80, 0, 0, 0, 1, 0, 0, 1, 02191118.6118-7-1-4-2-3-6-9-50, 0, 0, 1, 0, 1, 0, 0, 02221506.1To sum up, TDK is the main constraint of mission. The TDK number influences replenishing frequency and TDK redundancy and then impacts the fuel consumption. In general, high replenishing frequency means large fuel consumption. But in special conditions, the lowest replenishing frequency may not mean minimum fuel cost, because more TDKs carried by service satellite will lower the frequency of replenishment, while, in turn, they may cause higher mass in orbit transfer. When the replenishing frequency is the same, high TDK redundancy causes more fuel cost. The two indexes (replenishing frequency and TDK redundancy) almost determine the fuel cost in case the target debris is specific. Therefore, for the revolver mode, the TDK number the service satellite carries should be designed beforehand with suitable replenishing frequency and minimum TDK redundancy. This is, in fact, take full advantage of TDK carried by service satellite and reduce orbital maneuver with small satellite mass.


## 5. Conclusions

The ADR mission is an effective solution for GEO debris removal. In this paper, the mission planning for ADR mission of revolver mode is studied. A mathematical optimization model for revolver mode is proposed, considering TDK number, mission sequence, and orbital transfer time as constraints of mission. Two-level optimization is organized to solve this problem, including low-level for finding optimal transfer orbit using APSO and up-level for finding best removing sequence using IGA. Compared with ES method, the simulations demonstrate the efficiency of the two-level optimization.Experimental results indicate that the proposed method could solve the mission planning problem in ADR successfully. TDK is the main constraint of mission. The number of TDKs carried by the service satellite influences the fuel consumption by impacting the replenishing frequency and TDK redundancy. More TDKs carried in service satellite will lower the frequency of replenishment while in turn causing large TDK redundancy and higher mass in orbit transfer. Fewer TDKs carried by service satellite will cause high replenishing frequency. When the initial mass is the same, higher frequency of replenishment means larger fuel consumption, and when the replenishing frequency is the same, higher TDK redundancy causes more fuel cost. To reduce fuel consumption, the quantities of TDKs carried by service satellite should be optimized and designed with suitable replenishing frequency and minimum TDK redundancy.



---

*Source: 8284022-2021-12-02.xml*