A Novel Two-Level Optimization Strategy for Multi-Debris Active Removal Mission in LEO

: Recent studies of the space debris environment in Low Earth Orbit (LEO) have shown that the critical density of space debris has been reached in certain regions. The Active Debris Removal (ADR) mission, to mitigate the space debris density and stabilize the space debris environment, has been considered as a most effective method. In this paper, a novel two-level optimization strategy for multi-debris removal mission in LEO is proposed, which includes the low-level and high-level optimization process. To improve the overall performance of the multi-debris active removal mission and obtain multiple Pareto-optimal solutions, the ADR mission is seen as a Time-Dependant Traveling Salesman Problem (TDTSP) with two objective functions to minimize the total mission duration and the total propellant consumption. The problem includes the sequence optimization to determine the sequence of removal of space debris and the transferring optimization to define the orbital maneuvers. Two optimization models for the two-level optimization strategy are built in solving the multi-debris removal mission, and the optimal Pareto solution is successfully obtained by using the non-dominated sorting genetic algorithm II (NSGA-II). Two test cases are presented, which show that the low level optimization strategy can successfully obtain the optimal sequences and the initial solution of the ADR mission and the high level optimization strategy can efficiently and robustly find the feasible optimal solution for long duration perturbed rendezvous problem.


Introduction
The space debris is turning into one of the critical problems which strongly affect the design and the operation of current and future space missions. As the number of uncontrolled fragments increases, such as the rocket upper stage, the non-functional satellites, etc., the mission operations of currently in-orbit satellites are forced to perform collision avoidance maneuver regularly [Schaub, Jasper, Anderson et al. (2015); Felicetti and Emami (2016)]. If the maneuver of the in-orbit satellite is not timely, the collision will lead to serious results, such as failure and disintegration of the satellite. In 2009, the unintentional collision between the non-functional Russian satellite Kosmos 2251 and the operational US satellite Iridium 33 which occurred at an altitude of 789 km [Jankovic, Paul and Kirchner (2016)]. The collision has led to the production of a huge amount of new cataloged fragments in space. In fact, the fragments of collision event among the objects already in space have driven the evolution of the environment over several decades and result in an exponential increase of the cataloged fragments. A self-sustaining, cascading process is known as the Kessler syndrome [Lion andJohnson (2008, 2009); Lion (2011)]. The Kessler syndrome is first predicted by Kessler et al. in 1978 [Kessler andCour-Palais (1978)], and this phenomenon eventually become the most important long-term source of the space debris. To mitigate this phenomenon, space debris mitigation guidelines (SDMG) have been proposed by the IADC. Nevertheless, recent studies have shown that the space debris mitigation measures alone are not sufficient to guarantee the long-term utilization of some important orbital regions and to stabilize the current space debris environment [Lion, Johnson and Hill (2010)]. The ADR mission, to mitigate the space debris density and stabilize the space debris environment, has been considered as a most effective method. According to the recent studies [Lion andJohnson (2008, 2009); Lion, Johnson and Hill (2010)], the 5-10 space debris in the crowded altitudes and the inclination bands are removed per year to stabilize the debris population. In the ADR mission, the Orbital Transfer Spacecraft (OTS) must execute several major mission phases, such as launch, phasing, far range rendezvous, close range rendezvous, capturing and de-orbiting. The orbital rendezvous between the space debris and the OTS is a necessary process for realizing the ADR mission, and the design of a phasing strategy is very important for a rendezvous mission. A good phasing strategy can save propellant and improve the safety of rendezvous missions [Fehse (2003)]. For reasons of launch cost, the number of debris removed in the same ADR mission should be set as high as possible, while the propellant of the OTS is limited. The OTS repeatedly rendezvous with the debris, service it and then drag into the graveyard orbit to release it until all debris are removed. The design of phasing strategy for a limited-time, multi-debris ADR mission leads to a hybrid optimization problem, which includes the sequence optimization to determine the removal sequence of the space debris and the transferring optimization to define the phasing maneuvers. The optimization of multi-debris ADR mission in LEO has been the subject of extensive research. Zuiani et al. [Zuiani and Vasile (2012)] have solved a 5 debris ADR mission for an OTS using electric propulsion for both orbital transfer and debris processing. Bonnal et al. [Bonnal, Ruault and Desjean (2013)] developed four strategies disposing the debris in the sun-synchronous orbit (SSO). Chemical and electric propulsion systems were analyzed and evaluated while mission planning was not a focus of their research. Missel et al. [Missel and Mortari (2013)] provided a path optimization strategy for ADR mission, focusing on the proposed Space Sweeper with Sling-Sat (4S) mission. Madakat et al. [Madakat, Morio and Vanderpooten (2013)] proposed a bi-objective time dependent travelling salesman problem model for the problem of optimally removing debris and use a branch and bound approach to deal with the ADR mission. Cerf employed branch and bound algorithm to optimize the debris selection and the trajectories of space debris collecting mission with the impulsive maneuver [Cerf (2013)] and the low thrust maneuver [Cerf (2015)]. The selection of debris for removal in SSO region and guidance to the selected debris using low-thrust propulsion is studied by Olympio et al. [Olympio and Frouvelle (2014)]. The mission planning of LEO active debris removal problem considering some complex constraints (communication time window constraint, terminal state constraint and the time distribution constraint) is studied by Yu et al. [Yu, Chen and Chen (2015)]. A bi-objective Debris Active Removal mission was optimized using a Branch & Bound method with a dedicated algorithm by Berend et al. [Berend and Olive (2016)]. The mission is seen as a Time-Dependant Traveling Salesman Problem (TDTSP) with two objective functions to minimize the total mission duration and the total propellant consumption. Yang et al. [Yang, Tang and Jiang (2018)] uses the "Branch & Bound" method for solving the complex GTOC8 problem. As can be seen from the above discussions, the existing studies have two obvious deficiencies. First, the sequence optimization problem to determine the removal sequence of the space debris is generally used by the enumeration algorithm, either explicit or implicit. The explicit enumeration algorithm (brute force approach) is used by Zuiani et al. [Zuiani and Vasile (2012) ;Bonnal, Ruault and Desjean (2013) ;Yu, Chen and Chen (2015)]. The implicit enumeration algorithm (the branch and bound method, the branch and prune algorithm) is used by Madakat et al. [Madakat, Morio and Vanderpooten (2013); Cerf (2013); Olympio and Frouvelle (2014); Berend and Olive (2016)]. For the long duration of one year in the ADR mission, the search spaces of the maneuver impulses and burn times are very large, so the sequence optimal strategy consumes a large amount of the computing time. The total number of possible debris sequence is n!, if the search domain around good solutions is considered as the focus, the Branch & Bound algorithm becomes less efficient [Berend and Olive (2016)]. Yang et al. [Yang, Li and Baoyin (2015)] introduces a method for gravity-assisted low-thrust asteroid mission planning in a two-level approach. When the simplified optimization model is used, the terminal conditions for the rendezvous cannot be satisfied. Second, the transfer strategy adopted obscures the relationships between the two optimization objectives, and the characteristics of the multi-debris active removal problem are not clearly shown either. Therefore, an appropriate optimization strategy and effective optimization algorithm are eagerly sought to successfully solve this complex problem. In this paper, a two-level optimization strategy for multi-debris removal mission in LEO is proposed, which includes the low-level and high-level optimization process. To improve the overall performance of the multi-debris active removal mission and obtain multiple Pareto-optimal solutions, the ADR mission is seen as a Time-Dependant Traveling Salesman Problem (TDTSP) with two objective functions to minimize the total mission duration and the total propellant consumption. The ADR missions in the different regions of LEO are analyzed, in which the orbital inclination is approximately 99° in one region and 82° in the other region. The remainder of this paper is organized as follows. In Section 2, the multi-debris removal mission scenario and the space debris of interest are given. Section 3 introduced the dynamics modeling of the space debris and the OTS and then the multi-objective multi-debris rendezvous two-level optimization strategy is built in Section 4. After that, the model is solved by a multi-objective genetic algorithm (MOGA). The optimal Pareto solution is successfully obtained by using the non-dominated sorting genetic algorithm II (NSGA-II) for the selected removal sequence. Two representative numerical problems are used to validate the proposed method, to analyze the characteristics of the multi-debris rendezvous missions in Section 5. Finally, the conclusions are given in Section 6.

Problem description 2.1 Space debris of interest
According to the studied results of Bastida et al. [Bastida and Krag (2009);Castronuovo (2011)], there are three most catastrophic collision orbit regions in LEO, which have been identified as most critical because of the large concentration of mass. In these regions, the Kessler syndrome is most likely to be triggered [Lion andJohnson (2008, 2009);Lion (2011)]. Therefore, an efficient ADR mission is critical to remove the debris in these regions. The semi-major axis and inclination ranges of the debris in these regions have the characteristics as following: (1) a=7378±100 km, i=82±1°. In this region, the highest number of catastrophic collisions will occur in the next 200 years [Bastida and Krag (2009);Castronuovo (2011)].
(2) a=7178±100 km, i=99±1°. Region (2) corresponds to sun-synchronous conditions, which is very valuable for the commercial and scientific application. (3) a=7228±100 km, i=71±1°. In this region, the debris mainly is the spent rocket bodies with large mass and cross-section. Due to the highest number of catastrophic collisions in region (1) and the high commercial interest and scientific application in SSO and their intensive use in region (2), especially for remote sensing missions, the debris to be removed in these regions are considered as the priority targets in the ADR mission in this paper.

Multi-debris active removal mission
In the ADR mission, the OTS executes phasing maneuvers to reduce the difference of the semi major axis, the inclination and the right ascension of the ascending node (RAAN) for the different debris. Considering the n debris in the most catastrophic collision orbital region and locating on different orbits in LEO, the OTS is required to remove the n debris one by one until the n selected debris has been removed. In the process, the tasks are executed by the OTS which includes: (a) Performs phasing maneuver and comes to a hold position for the selected debris. (b) Inspects the debris and estimates the pose of the debris. (c) Deploys a de-orbiting kit in the direction of the debris and ensures the de-orbiting kit completes the close approximation, the attitude synchronization maneuver, adhesion, de-tumbling and de-orbiting. When the de-orbiting kit fails to attach to the debris, the other de-orbiting kit is deployed by the OTS. (d) Transfers to the next selected debris. The orbital elements of debris are denoted by D i={ai, ei, ii, Ωi, ωi, θi : i=1, 2, 3,…, n} and the initial orbital elements of the OTS are denoted by S0={a0, e0, i0, Ω0, ω0, θ0}, where a is the semi-major axis, e is the eccentricity, i is the orbital inclination, Ω is the RAAN, ω is the argument of perigee and θ is the true anomaly, as shown in Fig. 1. The whole operation process of ADR mission can be divided into two parts: transferring and servicing. Phase (a) is defined as transferring and phases (b) and (c) can be grouped into servicing. Phase (d) is seen as the transferring part of the next debris. The specific ADR process is shown in Fig. 2. The OTS is launched into the initial orbit at t0 and the servicing of the n-th debris is completed at tf, then the ADR mission ends. According to the time, the whole ADR mission can be divided into n phases. For the removal mission of the k-th debris, the time of the mission is defined as follows: t2k-1 refers to the time when the OTS leaves the (k-1)-th debris, t2k refers to the time when the OTS arrives and begins to serve the k-th debris, t2k+1 indicates the service ending time for the k-th debris. When k=n, there is t2n+1=tf. The transferring time can be expressed as t2k-t2k-1 for the phasing maneuver time of the OTS, and the servicing time can be expressed as t2k+1-t2k. Due to consider scenarios with de-orbiting kit, the de-orbiting kit is performed attaching to the debris and provide the necessary velocity increment and drag the debris into the final disposal orbit [Covello (2012)]. Since the large range phasing maneuver only needs to be executed by the OTS, the de-orbiting kit strategy is more economical to the propellant consumption of the OTS [Berend and Olive (2016)]. Summing up, the goal of the mission is to find the optimal service sequence and rendezvous path to visit all debris such that the total propellant consumption and the total mission duration are optimal, and the terminal constraint conditions should be satisfied meanwhile. So, the ADR mission can be seen as a Time-Dependant Traveling Salesman Problem (TDTSP) with two objective functions to minimize the total mission duration and the total propellant consumption.

OTS Initial Orbit
De-orbiting

Servicing phase
For the cooperative target, the de-orbiting kit must execute many complicated operation in the space rendezvous and docking mission, such as the close proximity, relative navigation, station maintenance. Zhang et al. [Zhang and Park (2013)] have studied the fuel consumption and the time duration of the multiphase orbital rendezvous mission and obtained the Pareto-optimal solutions satisfying engineering constraints. Because the debris is the non-cooperative targets, the complicated operation, such as attitude synchronization, debris capture and the stability of the combination, must be also completed by the de-orbiting kit. It is difficult to deal with the tumbling debris in the ADR mission by the de-orbiting kit, which can even lead to a failure of the mission. When the de-orbiting kit fails to attach to the debris, the OTS can designate and deploy the other de-orbiting kit for the same debris. The time consumption of servicing the non cooperative targets will be much longer. In the studies of Cerf et al. [Cerf (2013); Berend and Olive (2016)], the duration time is identical, but the characteristic parameter α is introduced in this paper to describe the effects of the servicing time for the tumbling space debris. The servicing time of the k-th space debris can be expressed as: where αk is the scale coefficient which is related to the mass mk, the moment of inertia matrix Jk and the angular velocity ωk of the space debris. ∆Tc is the reference time which is obtained by Zhang et al. [Zhang and Park (2013)]. When the service time is assumed to be same, all αk are a constant. When the service time is very small compared with the total mission duration, the service time can be ignored, ∆Tk=0.

Dynamic model including J2 perturbation
Due to the debris uncontrolled and subjected to perturbations (Earth gravitational perturbations, Sun and Moon attraction, Solar pressure radiation, geomagnetic field, etc.), it is valid for a long-duration, non-coplanar rendezvous ADR mission, in which perturbations such as J2 should be considered. In two different optimization strategies, the two different dynamic models are used, which are the osculating orbital dynamic model and the average dynamic model.

Osculating orbital dynamical model
The osculating orbital dynamical model in Cartesian coordinate system is widely used in the spacecraft orbit maneuver problem (such as orbital optimization design) with the position vector r and the velocity vector v as the state variable. But the position and the velocity in the Cartesian coordinates are the fast variables, for the long duration, the smaller integral step size is needed to be chosen in the orbital integral, which greatly increases the computation time. Under the effect of the J2 perturbation, the perturbation acceleration is small and consequently the solution of the motion can be described in terms of "almost constant" orbital elements. If the Gauss orbital equation is used to describe the motion of the space debris and the OTS, the classical orbital elements (a, e, i, Ω, ω, θ) are used. However, these elements exhibit singularities for e=0 and i=0 0 , 180 0 .
Therefore, the equinoctial orbital elements (p, f, g, h, k, L) that avoid the singularities in the classical orbital elements have been introduced by Broucke et al. [Broucke and Cefola (1972); Walker, Ireland and Owens (1985)]. At the same time, only L is the fast variable, so the dynamic equations are used to solve a low-thrust earth orbit transfer problem as described by Betts et al. [Betts and Huffman (1993)]. The equinoctial orbital elements to formulate the osculating orbital dynamical equations in the matrix form as follows: where y= [p, f, g, h, k, L] T is the equinoctial orbital element matrix. M is a 6×3 matrix and D is a 6×1 matrix, which are expressed as follows: The relationship between the equinoctial orbital elements and the classic orbital elements are given as: The equinoctial orbital elements y are related to the Cartesian state (r, v) according to the expressions The J2 perturbation acceleration af by the equinoctial orbital elements y can be expressed as: In the ADR mission, the transfer time of the OTS is far longer duration. In order to fully understand the relationship between the two objectives, the multiple impulses and multiple revolutions Lambert rendezvous algorithm considering the effect of J2 perturbation is used. The solution of Eq.
(2) based on the numerical integration method can be expressed as: The position vector r and velocity vector v are obtained by Eqs. (6) and (7), which can be given as: In the every impulse position of the transferring phase, an impulsive ∆vkl is applied, where k is the k-th transferring phase, l is the l-th impulse position, and k≤n, l≤N. n and N is the debris number and the impulse number of every transferring phase, respectively. When the superscript "-" indicates the state before an impulse and the superscript "+" indicates the state after an impulse, in the instantaneous t, the position vector and the velocity vector satisfy the following as: For the intermediate impulse, according to Eq. (10), the following equations must be satisfied For the first impulse in every transferring phase, the equations are given by If the initial coast time is not considered, the motion state of the OTS is consistent with that of the (k-1)-th debris at the initial time of the transferring phase. The equations at the final conditions for every transferring phase are given by If the finial coast time is not considered, the motion state of the OTS is consistent with that of the k-th debris at the finial time of the transferring phase.
For the multiple revolutions Lambert rendezvous problem, there exist 2Nmax+1 mathematical solution for Nmax revolutions [Shen and Tsiotras (2003)]. However, not all of these solutions are feasible. In actual engineering, the solutions require the perigee height to be greater than a minimum value hmin to avoid the influence of the earth's atmosphere, such as hmin>200 km, and the apogee height to be lower than a maximum value hmax to avoid expensive changes in large eccentricity [Yang, Luo, Zhang et al. (2015) ].
The feasible number of revolutions satisfies the following equation where N(amin) and N(amax) are the numbers of revolutions corresponding to the minimum and maximum semi-major axis.

Averaged dynamic model
For the space debris and the OTS in the near-circular LEO, the average effect of the J2 perturbations on the orbit semi-major axis, the eccentricity and the inclination is very small in the long duration [Cerf (2013); Howard (2005); Vallado (2007)]. The space debris in the LEO is thus still moving on the near-circular orbits, at very close altitudes and inclinations. Considering the long-duration secular effect of the J2 perturbation, the averaged dynamic equations of the space debris and the OTS can be expressed as follows: where µ is the gravitational constant and µ=3.986×10 14 m 3 /s 2 , J2 is the first zonal coefficient and J2=1.086×10 -3 , Re is the mean equatorial radius of the Earth and Re=6378.137 km [Howard (2005)]. It is very difficult to optimize the mixing problem by using the osculating orbital dynamical model for two reasons. First, the search spaces of the optimal variable, such as impulse amplitude, burn time and the transferring time, are very large because of the long duration. Second, the total number of possible debris sequence is n!. If the search domain around good solutions is considered as the focus, the Branch & Bound algorithm becomes less efficient [Berend and Olive (2016)]. Therefore, a two-level optimization strategy are eagerly sought to successfully solve the complex mixing problem, which includes the sequence optimization to determine the removal sequence of the debris and the transferring optimization to define the orbital maneuvers.

Two-level optimization strategy
The optimal transfer strategy depends on the initial and final orbits, the mission constraints and the OTS capabilities. In the general case, finding the optimal control law is a challenging task. In order to fully understand the relationship between the total fuel consumption and the time duration, two different type problems have to be made. First, the removal sequence of the debris should be optimized. Second, each transferring stage should also be optimized. It should be noted that there exists an infinity of possible transfer stages for a given trajectory linking two given debris. Due to the bi-objective optimization, a Pareto frontier can be obtained, which gives an interesting insight on the trade-off between the time duration and the fuel consumption. In the two-level optimization strategy, the mixing optimization problem is divided into the low-level optimal process and the high-level optimal process. The low-level optimal process determines several optimal sequences and the time interval, as the initial guess of the high-level optimal process. The high-level optimal process finds the set of the optimal control input that could trigger the optimal trajectory with the Pareto-optimal solution, and the rendezvous boundary conditions can be satisfied. The specific solution procedure is denoted as follows.

Low-level optimization process
In the low-level optimization strategy, the main objective is to determine several optimal sequences and the time interval as the initial input of the high-level optimal process. In order to reduce the computation time of the feasible solution and obtain several optimal debris removal sequences as fast as possible, several hypotheses are introduced: (a) The impulse maneuver is used by the OTS, (b) The circular drifting orbit strategy is used for the correction of the RAAN, (c) The correction of the true anomaly and the perigee argument in the orbit transfer are negligible. During the whole removal process of the debris, corresponding to the removal sequence D=(D1, D2, …, Dn), the debris one by one until the n-th debris is removed by the OTS. For the removal of the k-th debris, the time sequence of the OTS is tk=[t2k-1, tkd1, tkd2, t2k, t2k+1], where tkd1 is the start time of the drifting phase, tkd2 is the end time of the drifting phase. Due to the whole ADR mission starts at t1 and ends at tf, there is t2n+1=tf. In the interval of time t2k and t2k+1, the de-orbiting kit completes operations including the close proximity, synchronization, capture and the stability of the combination. When the de-orbiting kit completes the service of the k-th debris at the time t2k+1, the transferring phase starts for the (k+1)-th debris. Repeating these operations, the OTS would rendezvous with debris which are selected and finish the whole mission in tf. According to the hypothesis, the states of the OTS and the k-th debris at time t can be describe as:

t e t i t t t a t e t i t t
The orbital elements of the drifting orbit for the k-th space debris removal can be denoted as: In the transferring stages of the OTS, the Hohmann transfer strategy is used [Cerf (2013); Berend and Olive (2016)]. In the first and second transfer stage, the velocity increments are [∆vk1, ∆vk2] and [∆vk3, ∆vk4]. In the case an inclination change is performed simultaneously with a shape change as shown in Fig. 3, the plane change accompanied by velocity change is the most efficient, so the impulsive maneuver is assessed as the norm of the vectors difference In order to reduce the fuel consumption of the OTS, the drifting time ∆tkd is calculated so that the RAAN of the drift orbit and the debris orbit are equal at the end of the drifting stage. As a time constraint, it is calculated using the following formula The servicing time of the k-th debris can be expressed as: (25) In this time interval, the de-orbiting kit completes the near range approximation, attitude synchronization, capture, stability of the combination and other operations. The time interval is assumed to be related to the mass, the moment of inertia, the relative velocity and the tumbling angular velocity of the tumbling debris. For different debris, the servicing time interval is different, but it is a constant ∆Tk in the servicing phase of the debris.
The main objective of the low-level optimal process is to determine several optimal sequences and the time interval as the initial guess of the high-level optimal process. For the ADR mission, the optimal variable vector of the mixing problem in low-level optimization process is expressed as: where Dj is the j-th optimal removal sequences and [Ed]j=[E1d, E2d,…, End]j is the drifting orbit parameters vector of the j-th optimal sequences. The low-level optimization process of the ADR mission can be modeled as: where m is the number of the optimal removal sequences that satisfies the minimum fuel consumption in one year, ∆Tmax is the maximum removal time in the ADR mission, ∆vmax is the maximum characteristic velocity of the OTS. In order to prevent the optimal solution is lost before the high-level optimization process and reduce the computing time of the high-level optimization process, m is set to five for n=5 in this paper.

High-level optimization process
According to the results obtained in the low-level optimization process, the sequence D of the space debris removal is fixed in the high-level optimization process, and the optimization problem is a piecewise optimal control problem with some state constraints. Along with the time axis, the generic continuous-time state of the OTS for the k-th space debris could be abstractly depicted as Fig. 4.

Design variables
The design variables include the impulse number N, the impulse time interval ∆tkmi and the velocity increment ∆vkmi. As shown in Fig. 4, there are n×(4N-5) design variables for an N-impulse maneuver plan in the n-debris active removal mission: where tkmi is the i-th impulse time in the k-th transferring phase. ∆vkmi is the i-th velocity increment vector in the k-th transferring phase, Nk is the impulse number of the k-th transferring phase. In order to save the calculation time, the impulse number can be set to a constant, such as Nk=2, 3 or 4.

Objective function
In this paper, the ADR mission contains two kinds of objective functions, which are the total time duration and the total propellant consumption, respectively. In the transferring stage, when the OTS transfers from the i-th debris to the j-th debris, the fuel consumption ∆vij and duration ∆Tij are required. Based on the definition of the time in the ADR mission, as shown in Fig. 2 and Fig. 4, the transferring time can be expressed as: (29) The first objective is the total time duration, which can be expressed as: where eij is a binary selection variable for the switching from the Di to Dj, and eij equals to 0 or 1. ej is a binary selection variable for the service to the j-th debris, and ej equals to 0 or 1. When the service time is considered, the ej is equal to 1, otherwise, ej is equal to 0. When the service time of all space debris is not considered, the Eq. (30) can be simplified as: The second objective is the total propellant consumption of the ADR mission, which can be expressed as: For the different reasons, the total time duration and the total fuel consumption are constrained. For the total duration, it determines the mean number of space debris that is removed each year. According to the recent studied [Lion andJohnson (2008, 2009); Lion (2011)], the 5-10 space debris are removed per year to stabilize the debris population.
Therefore, the maximum value of ∆T is one year in this paper. For the total propellant consumption, it determines the burn duration of the OTS's rocket engine or the security margin of the OTS. The maximum velocity increment ∆vmax is provided by the OTS, so the total velocity increment must be less than ∆vmax, i.e., ∆v≤∆vmax. In the ADR mission, it is interesting to seek the shortest mission time and the lowest propellant consumption. Indeed, the shorter mission time is used, the more space debris is removed in one year. The lower fuel is consumed, the longer time is used by the OTS in ADR mission.

Constraints
The duration between two adjacent maneuvers should be long enough that the trajectories of the spacecraft can be properly measured, so as: where ∆tmin is the needed duration between two adjacent maneuvers. When the OTS arrives to the k-th debris at t2k, the OTS's state must coincide with the k-th debris's state. So, the position and velocity of the OTS must satisfy a set of linkage conditions Associated with any pair (Dj, Dk), the transferring time should satisfy where tkmin describes the required minimum time for transferring, tkmax gives the permitted maximum duration for the k-th debris. These values can be set reasonably according to the time values determined in the low-level optimization process. To deal with the nonlinear equality constraints presented in Eq. (34), the perturbed Lambert algorithm is used, so the last two impulses can be obtained as follows:

Optimization model
Overall, the high-level optimization process of the ADR mission can be rewritten as follows: ;

Optimization algorithm
To solve the multi-objective optimization problem, the existing methods can be roughly categorized into two classes: aggregative methods and Pareto-based methods. The aggregative method is to convert multiple objectives into single objective based on the assumption [Youssef and Eimaraghy (2008)]. On the contrary, the Pareto-based method is to provide a list of interesting trade-offs between the objectives rather than a lone solution supposing [Deb, Pratap, Agarwal et al. (2002)]. Since the two objectives under consideration are in conflict, the Pareto-based method is adopted. In the field of spacecraft trajectory optimization, NSGA-II is widely used in the low thrust orbit optimization design of the interplanetary spacecraft [Hartmann, Coverstone-Carroll and Williams (1998); Coverstone-Carroll, Hartmann and Mason (2000)], low thrust optimal transfer orbit [Lee, Allmen, Fink et al. (2005)], multi-objective optimal rendezvous [Luo, Lei and Tang (2007); Luo, Tang and Parks (2008)] and RLV optimal re-entry trajectory [Braun, Lupken, Flegel et al. (2013)]. Considering the successful application of NSGA-II for related problems, the NSGA-II is also employed to solve the mixing optimal problem in this paper. The NSGA-II is conceived as one of the famous Pareto-based multi-objective evolutionary algorithms (MOEAs). The main advantages of the NSGA-II approach compared with the other MOEAs are: (1) A fast non-dominated sorting approach ranks the solutions of a population by layers of non-dominated solution; (2) A crowding distance-based comparison operator is utilized to select solution for diversity preservation; (3) An elitism selection procedure is used to identify the best solutions from the individuals combining the parent and offspring populations with respect to fitness and spread. For further details about NSGA-II, the interested reader can refer to Deb et al. [Deb, Paratap, Agarwal et al. (2002)].

Solution for the two-level optimization strategy
In the two-level optimization strategy, the mixing optimization problem is divided into the low-level and high-level optimal process. The low-level optimal process determines several optimal sequences and the time interval, as the initial guess of the high-level optimal process. The high-level optimal process finds the set of the optimal control input that could trigger the optimal trajectory with the Pareto-optimal solution. The averaged dynamic model is used in the low-level optimization problem, which can only obtain an approximate solution, but cannot attain the desired state in the absolute dynamics system, and the rendezvous boundary conditions cannot be satisfied. To deal with this problem, the osculating orbital dynamical model is proposed to refine the approximate solution to a precise solution. The solution process of the two-level optimization strategy is summarized as follows: Step I. l=0, l is number of the remove sequence and l≤n!, n is the number of the debris, set the initial orbital elements of both the OTS and the debris, Es (t0) and Ek (t0), as well as the maximum removal time ∆Tmax, the maximum characteristic velocity ∆vmax, boundary values of design variables Ed, the number of the optimal removal sequences m, and the common parameters for the NSGA-II.
Step II. The Pareto frontier of each sequence is obtained using Eq. (27) according to the NSGA-II method. By comparison, the m optimal sequences Dj and the corresponding design variables (Ed)j are obtained.
Step III. j=1, set the needed duration between two adjacent maneuvers ∆tmin, the common parameters for the NSGA-II.
Step IV. The minimum duration time for transferring tkmin and the permitted maximum duration time describes tkmax for the k-th debris based on the estimate value (Ed)j of the low level optimization problem.
Step V. Solve the high level optimization problem of Eq. (37) using the NSGA-II, obtain the Pareto frontier of the j-th sequence, the absolute trajectories and the maneuver plan of the OTS.
Step VI. If j>m, the calculation is terminated and the optimal sequence Dopt and the corresponding design variables x is the final solution; else go to Step IV.

Numerical results
The ADR mission is a most effective method to mitigate the space debris density and stabilize the space debris environment. A self-sustaining, cascading process is known as the Kessler syndrome [Lion andJohnson (2008, 2009); Lion (2011)], and this phenomenon will eventually become the most important long-term source of the space debris. According to the recent studies [Lion andJohnson (2008, 2009); Lion (2011)], the 5-10 space debris in the crowded altitudes and the inclination bands are removed per year to stabilize the debris population. The five space debris with high priority has been extracted from an initial list of 477 candidate debris in LEO in the 500-1200 km altitude range. Due to the plane change maneuver are very costly in term of propellant consumption, the space debris with the close inclination values are been choose in the removal list. The debris selection criterions have been detailed in Braun et al. [Braun, Lupken, Flegel et al. (2013)]. In this study, the two test cases given here are the debris in two different regions for ADR mission. In the first region, the semi-major axis and inclination ranges of space debris have the characteristics as following: a=7178±100 km, i=99±1°. The region corresponds to sun-synchronous conditions, which is very valuable for the commercial and scientific application [Castronuovo (2011)]. In the second region, the semi-major axis and inclination ranges of space debris have the characteristics as following: a=7378±100 km, i=82±1°. The highest number of catastrophic collisions will occur in the next 200 years [Bastida and Krag (2009);Castronuovo (2011)].
The constants used in this study are Re=6378.137 km, J2=1.082626×10 -3 and µ=3.9860044×10 14 m 3 /s 2 . The removal duration is less than 730 days (two years) in the low level optimization process and 365 days (one year) in the high level optimization process, where t1=0, tf≤730 (365) days, and ∆T=tf−t0. Tabs. 1 and 2 list the initial orbital elements of the debris for the two different regions. The first test case was used by Cerf [Cerf (2015)], and the second test case was used by Yu et al. [Yu, Chen and Chen (2015)]. The first test case has differences of 200 km in a, 20 deg in Ω and 7 deg in θ; the second test case has differences of 23 km in a, 2 deg in Ω and 260 deg in θ. The scale coefficient αk is assumed to be related to the mass, the moment of inertia matrix and the tumbling which is different and the reference time ∆Tc=14400 s (4 hours) is obtained by Zhang et al. [Zhang and Park (2013)]. The bound of the service time is [50,840] hours for the total ADR mission in the Yu et al. [Yu, Chen and Chen (2015)]. In this study, the classical orbital elements and the Hohmann transfer strategy are used in the low-level optimization process. The equinoctial orbital elements and the perturbed Lambert algorithm are used in the high-level optimization process.  If the out-of-plane maneuver is used directly to correct the large initial RNNA difference, the propellant consumption would be prohibitively very high. For the OTS, this is obviously not appropriate. Instead, the natural orbital precession rate due to the Earth's J2 perturbation can be used to modify the RNNA difference of the orbits for the OTS and the space debris. This means that the RNNA difference can be indirectly corrected via an in-plane maneuver combined with long-duration orbital drifting under J2 perturbation [Cerf (2015); Olympio and frouvelle (2014); Yu, Chen and Chen (2015); Berend and Olive (2016)]. For the LEO short-duration near-coplanar rendezvous, the orbital drifting time would be very long-duration, so the RNNA difference is corrected using the out-of-plane maneuvers.
According to the characteristic of the ADR mission in the low-level optimization process, the first Hohmann transfer needs to be performed at the beginning time to drifting orbit so as to correct the RNNA difference, and the second Hohmann transfer needs to be performed at the rendezvous time to correct the semi-major axis difference between the OTS and the space debris. In actual engineering, the drifting orbit require the semi-major axis to be greater than a minimum value to avoid the influence of the earth's atmosphere for the OTS, and to be lower than a maximum value to avoid expensive orbital maneuver in large difference of semi-major axis. Due to the plane change maneuver are very costly in term of propellant consumption, it can be estimated that 100 m/s of cross-track maneuver to correct 1 deg of inclination difference, so the inclination of the OTS with the close inclination values are been choose in the removal process. Based on the above analysis, the boundary values of the design variables in the low-level optimization process are set as in Tab. 3. The low-level optimal process determines the optimal sequences, the time interval and the increment of velocity, which are the initial guess of the high-level optimization process. The boundary values of the design variables in the high-level optimization process are determined by initial guess values. The NSGA-II is used for optimization to calculate the Pareto frontier of the ADR mission, with the other algorithm parameter values being the same as those given in Tab. 4. Taking into account the randomness of NSGA-II, each case is run 10 times. By removing the repeated and dominated solutions, the Pareto optimal solutions of the problem are obtained.

Optimization results of the low-level optimization process
To validate the performance of the proposed low-level optimization model, the first test case is solved for four impulse maneuvers. The success rate for each configuration is 100%, which demonstrates that the low-level optimization model successfully locates feasible solutions in the experiments. In addition, the Pareto frontier of the bi-objective optimization problem shows that the NSGA-II is efficient and robust in obtaining the near-optimal solutions. In order to compare the result of the Cerf [Cerf (2015)], the Pareto-fronts for the different sequences (24135 and 24153) of the ADR mission using the low-level optimization process are illustrated in Fig. 5, in which the mission duration does not include the Servicing time. From Fig. 5, the results can be seen: (1) in general, with the increase of the time duration, the velocity increment of the OTS in the ADR mission would decrease; (2) with the increase of the time duration, the changes of the velocity increment for the different sequences of the ADR mission would be different; (3) for the different ADR sequences, such as 24135 and 24153, when the ADR time duration less than 336 days, the velocity increment of the sequence 24153 is less than that of the sequence 24135, the results are consistent with that in the Cerf [Cerf (2015)].
In order to obtain more clearly the effect of different ADR sequences on the velocity increments of the OTS, Fig. 6 shows the Pareto optimal solutions in the five different ADR sequences. In addition to the results obtained in Fig. 5, it is also showed that the velocity increments of the five different ADR sequences are bounding together in the time duration [335,375] days, that is, the velocity increments are almost equal. The result is mainly due to the true anomaly of the debris is relatively concentrated, which provides a good reference for the removal of concentrated debris. In order to further verify the effectiveness of the low-level optimization model, the second test case is solved for four impulse maneuvers. The Pareto-fronts for the five different sequences of the ADR mission using the low-level optimization process are illustrated in Fig. 7. From Fig. 7, the results can be seen: (1) with the increase of the time duration, the velocity increment of the OTS in the ADR mission would decrease, but the velocity increment of some sequences decreases very slowly, such as 45123; (2) the velocity increments of the five different ADR sequences are relatively scattered in the time duration [100,500] days, the result is mainly due to the true anomaly of the debris is relatively dispersed. As discussed earlier, the low-level optimal process determines several optimal sequences and the time interval, as the initial guess of the high-level optimal process. In the next section, the high-level optimal process finds the set of the optimal control input that could trigger the optimal trajectory with the Pareto-optimal solution, and the rendezvous boundary conditions can be satisfied.

Optimization results of the high-level optimization process
To solve the multi-objective optimization problem, the Pareto-based method is to provide a list of interesting trade-offs between the objectives rather than a lone solution supposing. In this paper, the Pareto optimal solutions have been calculated using the NSGA-II, which is conceived as one of the famous Pareto-based multi-objective evolutionary algorithms (MOEAs) [Deb, Pratap, Agarwal et al. (2002)]. Based on the results provided by the low level optimization process, the Pareto optimal solutions of the ADR mission is obtained by NSGA-II using the high level optimization process. For the second case, the Pareto-front of the three optimal sequences (15423, 54123 and 14523) is shown in Figs. 8-10, respectively. From these Figures, the results can be seen: (1) the optimal velocity increments are calculated, which is 186.94 m/s, 214.82 m/s and 226.47 m/s, respectively; (2) when the time duration is larger than a certain threshold, such as 154.33 day, 238.01 day and 175.64 day for three sequences respectively, the velocity increment decreases very slowly while the time duration increases; (3) for the impulse maneuver, the change of time duration will significantly increase the velocity increment of the OTS, so the points of the Pareto front are sparse and discontinuous. At the discontinuity of the Pareto front, the velocity increment of the OTS cannot be significantly reduced through the time duration increase.  (14523) It is important to note that the optimal solutions for the three different sequences shown in Figs. 8-10 are the subset of the solution spaces for the different sequences of the ADR mission that have been explored, so the Pareto-front actually depends on the multi-objective optimization method and the parameter setting of the given optimization method. Therefore, the Pareto-front is obtained using the multi-objective optimization methods, which is dominated by the true Pareto-front of the problem. The true Pareto front is unknown but can be approached using the better optimization algorithm and reasonable parameter settings. For the optimal solutions of the cases using the two-level optimization strategy, three different reasons guarantee the sub-optimal solutions approaching the true Pareto-front of the ADR mission. The first reason is the choice of the low-level optimal process, which uses a simplified optimal model, providing the several optimal sequences and the time interval as the initial guess of the high-level optimal process. The low-level optimal process is the key of a fast computation. The second reason is the high-level optimal process, which can trigger the optimal trajectory with the Pareto-optimal solution because of the more impulses are adopted, the J2 perturbation effects are consider and the rendezvous boundary conditions can be satisfied. The third reason is the choice of the NSGA-II, which have the main advantages compared with other multi-objective optimization algorithms, such as a fast non-dominated sorting, a crowding distance-based comparison operator and an elitism selection procedure.

Conclusion
A two-level optimization strategy for multi-debris removal mission in LEO is proposed, which includes the low-level optimization process and the high-level optimization process.
To improve the overall performance of the multi-debris active removal mission and obtain multiple Pareto-optimal solutions, the ADR mission is seen as a Time-Dependant Traveling Salesman Problem (TDTSP) with two objective functions to minimize the total mission duration and the total propellant consumption. The problem mixes the sequence optimization to determine the sequence of removal of space debris and the transferring optimization to define the orbital maneuvers. Two optimization models for the low-level optimization process and the high-level optimization process are built in solving the multi-debris removal mission, and the NSGA-II is employed to obtain the optimal Pareto solution. Numerical results show the following.
1) The two level optimization strategy and solution proposed in this paper can effectively address the ADR mission planning problem.
2) The low level optimization strategy can successfully obtain the optimal sequences and the initial solution for the ADR mission.
3) The high level optimization strategy can efficiently and robustly find the feasible optimal solution for long duration perturbed rendezvous problem. In future work, the following studies will be considered. First, the low-thrust transfer strategy and hybrid-thrust transfer strategy will be considered. Second, how to select the priority space debris which can minimize the collision risk and the velocity increment of the OTS. Third, the actual characteristic of the space debris will be considered.