# Preliminary Design of Debris Removal Missions by Means of Simplified Models for Low-Thrust, Many-Revolution Transfers

Authors: Federico Zuiani; Massimiliano Vasile
Journal: International Journal of Aerospace Engineering
 (2012)

---

## Abstract

This paper presents a novel approach for the preliminary design of Low-Thrust, many-revolution transfers. The main feature of the novel approach is a considerable reduction in the control parameters and a consequent gain in computational speed. Each spiral is built by using a predefined pattern for thrust direction and switching structure. The pattern is then optimised to minimise propellant consumption and transfer time. The variation of the orbital elements due to the thrust is computed analytically from a first-order solution of the perturbed Keplerian motion. The proposed approach allows for a realistic estimation ofΔV and time of flight required to transfer a spacecraft between two arbitrary orbits. Eccentricity and plane changes are both accounted for. The novel approach is applied here to the design of missions for the removal of space debris by means of an Ion Beam Shepherd Spacecraft. In particular, two slightly different variants of the proposed low-thrust control model are used for the different phases of the mission. Thanks to their low computational cost they can be included in a multiobjective optimisation problem in which the sequence and timing of the removal of five pieces of debris are optimised to minimise propellant consumption and mission duration.


---

## Body

## 1. Introduction

One of the most critical issues related to the exploitation of space around the Earth is the threat posed by space debris. Since the beginning of the space era in the late 1950s, an increasing number of man-made, inert objects has been orbiting the Earth. Recent statistics revealed around 15000 trackable objects, for a total of some 6000 tons of material. Some of these objects are simply spent upper stages of launch vehicles, some others are satellites which are no longer active due to failures or to having reached their end of life. Others, however, are the results of past collisions. It is easy to imagine that even a single collision between two objects is likely to generate tens of smaller objects as a result. The outcome of a collision in an already crowded environment could generate a cascade of collisions generating an exponentially increasing volume of space debris. In fact, the debris produced by a collision is itself likely to collide with other objects, thereby producing other debris which will generate further collisions, and so on. This chain reaction, known as theKessler Syndrome [1], occurs once the rate of generation of debris due to collisions or simple human-driven additions exceeds the natural debris removal rate. According to Kessler et al., this reaction is likely to be ignited once the object density in a certain orbital band reaches a critical point; once started, it will probably render most spacecraft in that orbital band useless within a matter of months or years.Recent guidelines issued by international spacer regulatory institutions such as the United Nations Committee for the Peaceful Uses of Outer Space (COPUOS) [2] and the Inter-Agency Space Debris Coordination Committee (IADC) [3] prescribe some actions to be followed by national or private agencies putting satellites into orbit in order to mitigate debris growth. For example, it is demanded that every new mission in Low Earth Orbit (LEO) must be planned such that the satellite itself must re-enter in the Earth’s atmosphere within 25 years after the end of the mission. Alternatively, for higher orbits like geostationary orbits, the requirement is for the spacecraft to be placed on a higher graveyard orbit. Measures like these, even if strictly applied (and at the moment compliance with them is on a voluntary basis), are only likely to slow down the accumulation of space debris around the Earth. Therefore, active removal actions will probably be needed in the near future to eliminate at least the most dangerous objects.There have been various proposals on how to remove inert objects from space. They can be generally classified in two major groups: contactless and with direct physical contact. In the latter category, one can find methods based on some form of docking with or capturing the object. Once the removing spacecraft and the piece of debris are attached, the latter is dragged into a re-entry trajectory or to a graveyard orbit. Technical problems related to the attitude state of motion of the piece of debris and the fragility of appendices and cover material (including paint) make this removal solution complicated. A potentially interesting solution is represented by Project ROGER [4], developed by EADS/Astrium with the support of ESA. Among contactless solutions, one can find what is commonly referred to as the space broom [5]. It entails irradiating the target object with a high-power laser which will induce sublimation of the surface material; the ejecta plume will then generate a low-thrusting acceleration which will slowly degrade the debris’ orbit until it reaches an altitude where atmospheric drag will accelerate its re-entry. Such a technique has the advantage that no physical contact is required, on the other hand, current proposals envisage the use of lasers installed on Earth and beaming through the atmosphere. The beam collimation and thrust time is, therefore, limited and this solution is effective for small-sized objects only. Recent proposals have demonstrated that the use of in-space laser systems might be more interesting even to remove larger objects [6]. Other proposals involve, for example, the use of electrodynamic tethers [7], inflatable balloons [8], which are meant to be lightweight and efficient but require, however, the physical attachment of the device to the target object and are, therefore, of difficult application to existing debris.A recent idea simultaneously proposed by Bombardelli and Peláez [9], Ruault et al. [10], and JAXA [11] suggested the use of a collimated beam of ions generated by a spacecraft flying in formation with the piece of debris. In this paper, this concept will be called Ion Beam Shepherd (IBS), using the name introduced by Bombardelli and Peláez. The effect of the ion beam is that of producing a thrusting force, equal in magnitude but opposite in direction, on both the IBS and the piece of debris. This force will induce a thrusting acceleration which can be controlled in order to modify the orbit of the piece of debris. A second ion engine is then fired in a direction opposite to the first one in order to keep the IBS spacecraft at a constant distance from the piece of debris. Among the advantages of this concept is the fact that it employs already existing and proven technologies; it does not require any contact with the target, and the fact that a single spacecraft can be used to fetch and deorbit multiple pieces of debris. In [6], one can find a similar concept that uses concentrated solar light instead of ions to generate a thrust and modify the orbit of debris.Assuming a scenario in which a single IBS needs to de-orbit multiple pieces of debris, one would need to solve an interesting mission design problem: the optimisation of the de-orbit sequence and trajectories for multiple target objects in minimum time and with minimum propellant. In the hypothetical mission scenario which is analysed in this work, it is assumed that a number of pieces of debris have been shortlisted as priority targets due to the threat they pose to satellites operating in LEO. For example, Liou and Johnson [12] propose some criteria to choose the object whose removal will be most effective to mitigate the risk of collisions. They underline that an effective removal strategy must be targeted first to large objects in crowded orbits up to 1500 km. Thus, a removal mission by means of an IBS spacecraft is planned to be launched from the Earth. Its task is that of removing five objects lying on different low Earth orbits. The design of such a mission is a complex optimisation problem, because it requires the computation of multiple low-thrust, many-revolution transfers. Therefore, this work proposes an approach to the fast estimation and optimisation of the cost and time duration of the fetch and de-orbit sequences. In past works, other authors have already proposed approaches to the design of Low-Thrust (LT), many-revolution transfers, based on analytical solutions to an optimal orbit raising problem under the assumption of small eccentricity [13–15] or on averaging techniques [16–18]. This work, proposes a different approach, based on a first-order solution of perturbed Keplerian motion. The approach in this paper aims at capturing the definition and optimisation of the thrusted arcs for each orbit without sacrificing computational speed. The approach can be classified as direct method for trajectory design as it does not derive the necessary conditions for optimality but translates the initial optimal control problem into an NLP problem.The paper is organised as follows: Section2 will briefly outline the IBS concept and in particular will outline how to compute the thrusting acceleration generated on a given target object; Section 3 will analyse an hypothetical mission profile for the removal mission and most important, Sections 3.1 and 3.2 will present in detail the proposed trajectory models. Section 4 will then show how the mission design problem can be then translated into as a series of multiobjective optimisation problems which are solved with a stochastic optimiser. The results are then presented and discussed.


## 2. The Ion Beam Shepherd Spacecraft

As shown by Bombardelli and Peláez [9], the concept behind the Ion Beam Shepherd is relatively simple and envisions employing a spacecraft provided with two sets of ion engines mounted along the same axis but in opposite directions (see Figure 1). The jet from one of the sets will be directed towards the piece of debris and will exert a thrusting force Fp1 on it. Due to Newton’s third law, an opposite force of same magnitude will also act on the spacecraft itself, but this component will be balanced by the thrust Fp2 provided by the other set of ion engines.Figure 1
Ion Beam Shepherd spacecraft.Since it is necessary to keep the Shepherd spacecraft at a constant distance from the debris, the thrustFp2 shall be such that the second derivative of the distance d between the two spacecraft is null:
(1)d̈=Fp2-Fp1mIBS-Fp1md=0.
Note that in (1) the acceleration terms due to the gravity of the central body have been neglected since it is assumed that the debris and the Shepherd are in close proximity and arranged in a leader-follower configuration. A more accurate and detailed modelling of the proximal motion dynamics of these two bodies is beyond the scope of this study. Thus, in the following sections, the IBS-debris combination will be treated as a point mass, in order to apply two-body dynamics. By rearranging the terms in  (1), one obtains:
(2)Fp2=Fp1(1+mIBSmd).
Under the assumption that the total propulsive power of the IBS spacecraft Ptot is constant and that the total propulsive thrust is proportional to it Ftot, one can write
(3)Fp1+Fp2=Ftot∝Ptot,
thus
(4)Ftot=Fp1(2+mIBSmd).
Therefore, the maximum acceleration acting on the IBS-debris combination can be computed as a function of the total available thrust Ftot:
(5)εIBS-debr=Fp1md=Ftot2md+mIBS.
It is assumed here to have a 1000 kg IBS spacecraft with a total available thrust of 0.5 N. Such a high thrust would correspond to a substantial power and propulsion system mass, however, this is deemed realistic if one considers that the payload of the IBS spacecraft is in fact its propulsion and power systems. Hence, the propulsion and power systems might be oversized compared to other applications in which ion engines are used for propulsion only. Note that the validity of the methodology proposed in this paper would not be affected even if lower thrust levels were considered. Thus, in this case, considering, for example, an 800 kg debris, the magnitude of the acceleration would be 1.923·10-7km/s2. If one considers instead the spacecraft alone, the acceleration achievable would be slightly higher, 5·10-7km/s2. Given this order of magnitude, the thrust acceleration can be considered as a perturbative force compared to the Earth’s gravitational force, and therefore, the analytical approach to the propagation of the low-thrust motion described in [19] can be applied.


## 3. Mission Profile

The objective of this study is that of optimising the performance and cost of a debris de-orbiting mission performed by a single spacecraft. As mentioned in the introduction, it is assumed that there are five pieces of debris of different masses and lying in circular orbits with different radii and orientations. It is assumed that the IBS spacecraft departs from a low-Earth parking orbit, rendezvous with the first object, transfers it to an elliptical re-entry orbit, rendezvous with the second object, transfers it to a second elliptical re-entry orbit, and so forth until all five pieces of debris are removed. One important issue is defining in which order the pieces of debris need to be de-orbited. In the following, all possible sequences are generated apriori and optimised one by one.Each fetch and de-orbit operation is split in two phases:(i)
Ade-orbit phase, in which the perigee of the orbit of the piece of debris is lowered such that the orbit will decay naturally in a relatively short time. In this study, it is assumed that this condition is met if the perigee altitude of the debris’ orbit is equal or lower than 300 km.(ii)
Atransfer phase, in which the IBS spacecraft rendezvouses with the next piece of debris (which lies on a circular orbit), after having abandoned the current piece of debris on an orbit with a 300 km perigee altitude.Given the magnitude of the available thrust acceleration, both phases require a spiral orbit transfer. If a direct transcription approach is used to optimise each spiral the number of parameters that needs to be defined is very high leading to high computational times. The latter fact would make the solution of a multiobjective optimisation of all possible de-orbiting sequences computationally intractable. Thus, in this paper, a simplified, highly efficient, trajectory model is proposed for each one of the two phases.


### 3.1. De-Orbiting Trajectory Model

The objective of the de-orbiting phase is that of lowering an initial circular orbit such that its perigee is equal or below 300 km, which basically translates into a perigee lowering manoeuvre. Therefore, it is appropriate to assume that in general, as soon as the initial circular orbit becomes slightly eccentric, one keeps thrusting around the apogee in order to lower the perigee. The thrust level will also be kept at its maximum in order to minimize gravity losses. Moreover, since the de-orbit condition is independent of the final orbit’s orientation, one can reasonably assume that the perigee lowering will be performed in-plane. In this sense, the only Keplerian parameters which need to be altered are the semimajor axis and eccentricity. By analysing the structure of Gauss’ variational equations:(6)dadt=2a2h(esin⁡θar+praθ),dedt=1h{psin⁡θar+[(p+r)cos⁡θ+re]aθ},didt=rcos⁡θhah,dΩdt=rsin⁡θhsin⁡iah,dωdt=1he[-pcos⁡θar+(p+r)sin⁡θaθ]-rsin⁡θcos⁡ihsin⁡iah,dθdt=hr2+1eh[pcos⁡θar-(p+r)sin⁡θaθ],
where a is the semimajor axis, e the eccentricity, i the inclination Ω is the Right Ascension of the Ascending node (RAAN), ω the argument of perigee, θ the true anomaly, p the semilatus rectus, and h the angular momentum; ar,aθ, and ah are the radial, transversal, and out-of-plane components of the thrust acceleration. If one considers the case of thrusting with maximum acceleration along arcs which are symmetrical around apogee (θ=π), one can see that the contributions to semimajor axis and eccentricity variations given by the ar components are negligible (since they are multiplied by sinθ). Therefore, a good suboptimal thrust direction can be obtained by imposing aθ as the only non-zero component of the thrust acceleration. Under these assumptions, the variation of Keplerian parameters will be [20](7)dadt=2a2hpraθ,dedt=1h[(p+r)cos⁡θ+re]aθ,didt=0,dΩdt=0,dωdt=(p+r)hesin⁡θaθ,dθdt=hr2+(p+r)sin⁡θaθ.
It should be noted that the terms of the variation of ω and θ which depend on aθ will also be very small due to the presence of sinθ integrated around θ=π. Figure 2 visualises the proposed pattern for thrusting arcs.Figure 2
Thrusting arc around apogee with thrust directed along transverse direction.In order to obtain a fast propagation of the thrusting arcs, the analytical propagation of perturbed motion with finite perturbative elements in time (FPET) derived in [19, 21] will be used. In order to employ FPET, one has also to assume that the thrust acceleration is constant around each thrusting arc, which is reasonable given the low propellant consumption per arc. This assumption ensures that, if the engine thrust is constant, the resulting acceleration can also be considered constant over short thrusting arcs.Motion propagation with FPET is based on a first-order analytical solution of perturbed Keplerian motion. In this formulation, the state is expressed in nonsingular equinoctial elements:(8)a,P1=esin⁡(Ω+ω),P2=ecos⁡(Ω+ω),Q1=tani2sin⁡Ω,Q2=tani2cos⁡Ω,L=Ω+ω+θ.
Assuming constant thrust-acceleration in the radial-transverse reference frame:
(9)a=[araθah]T=ε[cos⁡⁡αcos⁡⁡βsin⁡⁡αcos⁡⁡βsin⁡⁡β]T,
then one can obtain a first-order analytical expansion of the variation of equinoctial elements, parameterised in Longitude L and with respect to a reference longitude L0:
(10)a(L)=a0(L0)+εa1(L0,ΔL,α,β),P1(L)=P10(L0)+εP11(L0,ΔL,α,β),P2(L)=P20(L0)+εP21(L0,ΔL,α,β),Q1(L)=Q10(L0)+εQ11(L0,ΔL,α,β),Q2(L)=Q20(L0)+εQ21(L0,ΔL,α,β),t(L)=t00(L0,ΔL)+εt11(L0,ΔL,α,β),
where(11)L=L0+ΔL,
where a a0P10,P20,Q10,Q20 are reference values at L0 and a1,P11,P21,Q11,Q21 are first-order terms as reported in [19, 21]. In [19],it has also been shown that this analytical propagation scheme provides good accuracy along relatively long trajectory arcs.As explained above, the only nonzero component of the acceleration will beaθ and since the aim is obtaining a decrease of the orbit energy, it will also be in the negative direction. Therefore, the acceleration azimuth will be α=-π/2 and the elevation β=0 (since, as already mentioned, the motion will be within the initial orbit plane). The variation of equinoctial elements after an apogee thrusting arc will be given by
(12)E+=[aP1P2]L+=La+ΔLa=[aP1P2]L-=La-ΔLa+ε[a1(2ΔLa,π2,0)P11(2ΔLa,π2,0)P21(2ΔLa,π2,0)]=E-+εf(2ΔLa,π2,0),
where La is the apocentre longitude, L- and L+ are the longitudes at the start and end of thrusting, respectively. ΔLa is the semiamplitude of the apogee thrusting arc. Note that, given that β=0, there is no variation on Q1 and Q2 and thus they are omitted. The coasting time is computed from the last of (10) as(13)tthrust=t00(La-ΔLa,2ΔLa)+εt11(2ΔLa,-π2,0).Since the thrust magnitude and direction are fixed, the only free control parameter is the semiamplitudeΔLa for each orbit. In order to keep the number of decision variables to a minimum, the semiamplitude for each orbit is computed from a piece-wise linear polynomial interpolating a limited number of ΔLa,i over a number of orbits. The nodes ΔLa,i are equally distributed between orbit 1 and an arbitrary number of orbits (in this paper, 1200 was found to be adequate). In this paper, the number of interpolating nodes was limited to 2: ΔLa1 and ΔLaf.In order to evaluate the time and ΔV needed to de-orbit a piece of debris from its initial orbit with semimajor axis adebr0, given a set of decision (or control) parameters ΔLa1 and ΔLaf, the following procedure was implemented.(1)
Compute the set of initial Equinoctial parametersL0 and E0=[adebr0P10P20]T, where P10 and P20 will be null due to the fact that the initial orbit is circular.(2)
Initialise the number of orbits, the totalΔV and time of flight to zero:
(14)Norbit=0,ΔV=0,TOF=0.(3)
SetE-=E0 and Lcoast=L0.(4)
Initialise the mass of the IBS spacecraft:(15)mIBS=mIBS0.(5)
WhileNorbit is smaller than Norbitsmax,
(a)
Norbit=Norbit+1.(b)
Interpolate the amplitude of the thrusting arc in the current orbit, that is,ΔLa and compute L-=La-ΔLa and L+=La+ΔLa.(c)
Compute the accelerationεIBS-debr acting on the IBS-debris combination from (5).(d)
Compute the time of flighttcoast spent coasting from Lcoast to L-.(e)
Compute the equinoctial parameters after the thrusting arcE+ as in (12).(f)
Compute the current perigee radiusrp and if this is lower than the threshold r̅p=300 km proceed to step 6, otherwise, proceed to step g.(g)
Compute the thrusting timetthrust from (13) and update the total ΔV cost:
(16)ΔV=ΔV+εIBS-debrtthrust.(h)
Update the total time of flight:(17)TOF=TOF+tcoast+tthrust.(i)
Update the IBS spacecraft mass:(18)mIBS=(mIBS+2mdebr)exp⁡(-εIBS-debrtthrustIspg0)-2mdebr.(j)
SetE-=E+ and Lcoast=L+.(6)
Back track the value of the longitudeLf for which rp=r̅p and compute the related and tthrust from (13) and update TOF and ΔV accordingly. Compute the equinoctial parameters Ef at Lf from (12).At this point, one gets theΔV, the time of flight TOF, and the semimajor axis and eccentricity of the final orbit (which are easily computed from Ef). It is important to note that, given the simplifications introduced, once one sets the initial mass and orbit of the piece of debris, and the characteristics of the IBS propulsion system, that is, Ftot and Isp, the de-orbit depends exclusively on the mass of the IBS mIBS0 at the beginning of the de-orbit phase and the interpolating values for ΔLa, that is, ΔLa1 and ΔLaf. Therefore, it was decided to precompute the corresponding ΔV and TOF for a given set of these three parameters and for each piece of debris (i.e., for each mdebr and adebr0). Table 1 reports upper and lower bounds for mIBS0, ΔLa1 and ΔLaf, and the number of samples taken, equally distributed.Table 1
Bounds and number of samples for the de-orbit parameters.
mIBS0ΔLa1ΔLafLower boundMdry+100=350 kg00Upper boundMlaunch=1000 kgππSamples85050Given the limited number of decision variables, for each piece of debris, one has 20000 de-orbit instances to propagate. Since each instance requires typically 1·10−2 s of CPU time, with a code implemented in MATLAB and running on a 3.16 GHz, 4 GB desktop PC running Windows 7, the whole computation can be completed in roughly five minutes. The set of de-orbit ΔV and TOF is then used to build a response surface, or surrogate model, of the de-orbiting process. Figures 3(a) and 3(b) show examples of two-dimensional surface, respectively, for ΔV and TOF, with respect to a fixed mIBS0 of 300 kg. One can see that the two quantities show opposite trends, the ΔV being high when the TOF is low and vice versa. Figures 4(a) and 4(b) show the final semimajor axis and eccentricity, respectively. Note that the minimum TOF transfer corresponds to a quasicircular spiralling trajectory in which the IBS spacecraft is thrusting continuously. On the other hand, the minimum ΔV transfer corresponds also to the one with maximum final eccentricity.(a)ΔV and (b) TOF surfaces with respect to ΔLa1 and ΔLaf for mIBS0= 300 kg, adebr0 = 7128 km and mdebr = 120 kg.
(a)(b)(a) final semi-major axis and (b) eccentricity after de-orbit with respect toΔLa1 and ΔLaf for mIBS0 = 300 kg, adebr0 = 7128 km and mdebr = 120 kg.
(a)(b)Now, it is desirable that the surrogate model returns theΔV cost as a function of mIBS0, mdebr, adebr0, and TOF. From the available data relating the ΔV and TOF to the decision variables ΔLa1 and ΔLaf, one can derive the functional relationship between ΔV and TOF. Given a triplet mIBS0,mdebr,adebr0, each TOF value defines a level curve on the ΔLa1 and ΔLaf plane (see Figure 3(a)), which can be mapped into a set of ΔV values (see Figure 3(b)). Within this set, one can take the element with minimum ΔV. Thus, for each time of flight, between a minimum and a maximum, one can derive the corresponding minimum ΔV cost. A similar procedure is followed to find the functional relationship between the final semimajor axis and the TOF. Note that there is no need to do the same for the eccentricity given the fact that the final perigee radius is fixed at r̅p and, therefore, the final e can be computed from the final a. In this way one can build the two surrogate models: (19)ΔV=fΔV,interp(TOF,mIBS0,mdebr,adebr0),af=faf,interp(TOF,mIBS0,mdebr,adebr0).
Figures 5(a) and 5(b) show examples of tridimensional plots (mIBS0-TOF-ΔV and mIBS0-TOF-af, resp.) created by evaluating the surrogated models keeping adebr0 and mdebr fixed. In Figure 5(a), one can see that there is a large plateau region corresponding to large time of flights and a smaller region close to the minimum TOF where the de-orbit cost increases very steeply and the final semimajor axis in Figure 5(b) similarly decreases. The complete procedure for the creation of the interpolated de-orbit cost models requires few minutes of CPU time and once completed allows for a very fast estimation of the de-orbit cost. The surrogated models will be extremely useful in the multiobjective optimisation of debris removal sequences as it will be shown in the following sections.3D Plot of surrogate models foradebr0=7128 km and mdebr=120 kg: (a) ΔV; (b) af.
(a)(b)


### 3.2. Orbit Transfer Model

According to the scenario presented in Section3, after having left the debris on a re-entry orbit, the IBS will have to transfer to the orbit of the next debris and rendezvous with it. The design of such a transfer arc would normally require the solution of a fixed-time two point boundary value problem (TPBVP) which would be computationally very expensive given the high number of control parameters and constraints involved. A second simplified model was then created to quickly estimate the cost of a low-thrust multirevolution orbit transfer with boundary constraints. The approach and assumptions presented in this section are similar to those already introduced for the de-orbit model.First, given the limited acceleration provided by low-thrust propulsion systems, one should consider that the orbit transfer will require a high number of multiple revolutions around the Earth, typical in the range of hundreds to few thousands. In this sense, it is possible to argue that achieving the proper phasing to transfer from the initial to final orbit would not be a major issue. Even a small variation ofω and θ per revolution would be sufficient to attain the required orientation to rendezvous with the piece of debris. Moreover, it is important to bear in mind that, in order to de-orbit the previous debris in the sequence, the IBS spacecraft, started from a circular orbit which was subsequently modified into an elliptical one with perigee r̅p. Thus, it would be also possible to conveniently adjust the start point of the de-orbit procedure from the circular orbit in order to obtain the proper phasing once this is completed. For all these reasons, it is assumed that in this particular case, the phasing problem will have a negligible effect on the ΔV and time required to rendezvous with the next piece of debris in the sequence. Therefore, in the following, it is assumed that it is not necessary to match the arrival ω and θ computed with the simplified model with those of the target object. Matching the target inclination i and RAAN Ω, instead, cannot be ignored without introducing a considerable error in the ΔV cost. In order to match the inclination and RAAN difference, one need to take into account only the geometric angle between planes of the initial and final orbits, which is given by
(20)Δi=arccos⁡(-cos⁡i0cos⁡(π-if)+sin⁡i0sin⁡(π-if)cos⁡(Ωf-Ω0)).Thus, in order to account forΔi, the inclination of the initial orbit is fictitiously set to zero, while the final one is set at Δi. The matching of the RAAN is assured by performing the circularisation properly. The assumption is that the deorbiting of one piece of debris starts at a true anomaly such that the resulting elliptical orbit has the line of apses perpendicular with the line of the nodes of the following piece of debris. Since the orbits of the debris are assumed to be circular, it is always possible to start the deorbiting at the right true anomaly with minimum delay. This hypothesis will be discussed in more detail with some numerical examples in Section 4.With these assumptions, the main issue in designing the multirevolution transfer will be that of achieving the required change in the apogee and perigee radii in order to match those of the final orbit, and to achieve the required rotation of the orbit plane.The control of eccentricity and semimajor axis, required to match the target perigee and apogee altitudes, can be obtained by inserting two thrust arcs per revolution, one around the apogee and one around the perigee. This methodology is analogous to what was done in the previous section for the perigee lowering. In the same way, the radial component of the thrust acceleration is set to zero. The transverse component this time can have either positive or negative sign (α=±π/2) depending whether the perigee (or apogee) needs to be raised or lowered.Since a plane change is required, the out-of-plane component of the thrust acceleration is nonzero. Thanks to this, the control parameters can be reduced to the semiamplitude of the apogee and perigee thrusting arcs,ΔLa and ΔLp, the sign of the θ component of the thrust acceleration (i.e., the sign of αa,αp=±π/2), and the out-of-plane component in the same arcs, βa and βp. Define ΔLthrust as half the total thrusting arc length and rt as the ratio of ΔLthrust which is devoted to apogee thrusting. In order to have a parameterisation which accounts also for the sign of αa and αp, the following one is proposed:
(21)αa={π2ΔLt≥0,-π2ΔLt<0,αp={αa0≤rt≤1,-αa1<rt≤2,ΔLa={rt|ΔLthrust|0≤rt≤1,(2-rt)|ΔLthrust|1<rt≤2,ΔLp=|ΔLthrust|-ΔLa,
with
(22)ΔLthrust∈[-ππ],rt∈[02].To define the actual values ofΔLthrust and rt in each revolution, an interpolating strategy from a set of nodal values similar to the one used for the de-orbit model is adopted. In this case, however, the interpolated values will not be computed with respect to the current revolution number but with respect to the time. Again the number of interpolating nodes can be chosen arbitrarily and is set to 2 in this case, ΔLt1,ΔLtf,rt1,rtf. For βa and βp, it is chosen to have a constant value along the entire transfer. The thrusting pattern along each revolution is shown in Figure 6.Figure 6
Thrusting arcs around apogee and perigee.Given a set of control parameters[ΔLt1ΔLtfrt1rtfβaβp], a multirevolution transfer with specified duration TOF¯, departing from an orbit defined by [a0e00]T and targeted to an orbit defined by [afefΔi]T, is propagated according to the following procedure.(1)
Compute the set of initial Equinoctial parametersL0 and E0=[a0P10P20Q10Q20]T. Q10 and Q20 will be zero since the initial inclination is arbitrarily set to zero.(2)
Compute the set of target Equinoctial parametersE̅f=[a̅fP̅1fP̅2fQ̅1fQ̅2f]T. Note that P̅1f and P̅2f will be zero since in this case the target orbit is a circular one.(3)
Initialise the total ΔV and Time of flight to zero:
(23)ΔV=0,TOF=0.(4)
SetEp-=E0 and Lcoast,a=L0.(5)
Initialise the mass of the IBS spacecraft:(24)mIBS=mIBS0.(6)
WhileTOF<TOF¯.
(a)
Compute the interpolated values forΔLt and rt. Hence, calculate αa,αp,ΔLa, and ΔLp from (21).(b)
Compute(25)La-=La-ΔLaLa+=La+ΔLa,Lp-=Lp-ΔLpLp+=Lp+ΔLp.(c)
Compute the current acceleration acting on the spacecraft:(26)εIBS=FtotmIBS.(d)
Compute the time of flighttcoast,p spent coasting before perigee from Lcoast,p to Lp-.(e)
Compute the Equinoctial parameters after the thrusting perigee arcEp+ with an expression analogous to (12).(f)
Compute the thrusting time at perigeetthrust,p from (13). If (TOF¯-TOF)<tthrust,p, proceed to step g. Otherwise, break the iterative sequence and go to step 7.(g)
UpdateΔV and TOF:
(27)ΔV=ΔV+εIBStthrust,p,TOF=TOF+tcoast,p+tthrust,p.(h)
Update the IBS spacecraft mass:(28)mIBS=mIBSexp⁡(-εIBStthrust,pIspg0).(i)
SetEa-=Ep+ and Lcoast,a=Lp+.(j)
Compute the current acceleration on the spacecraft:(29)εIBS=FtotmIBS.(k)
Compute the time of flighttcoast,a spent coasting before apogee from Lcoast,a to La-.(l)
Compute the equinoctial parameters after the thrusting apogee arcEa+ as in (12).(m)
Compute the thrusting time at apogeetthrust,a from (13). If (TOF¯-TOF)<tthrust,a, proceed to step n. Otherwise, break the iterative sequence and go to step 7.(n)
UpdateΔV and TOF:
(30)ΔV=ΔV+εIBStthrust,a,TOF=TOF+tcoast,a+tthrust,a.(o)
Update the IBS spacecraft mass:(31)mIBS=mIBSexp⁡(-εIBStthrust,aIspg0).(p)
SetEp-=Ea+ and Lcoast,p=La+.(7)
Back-track the point at whichTOF=TOF¯ and compute the corresponding equinoctial parameters Ef=[afP1fP2fQ1fQ2f]T and update ΔV accordingly.(8)
Compute the mismatch between the actual final conditions and the target orbit:(32)Ceq=[af-a̅fef-e̅f=P1f2+P2f2-P̅1f2+P̅2f2if-i̅f=2(arctanQ1f2+Q2f2-arctanQ̅1f2+Q̅2f2)].Summarizing, the TPBVP has been reduced to an optimisation problem in the form:(33)min⁡xΔVs.t.Ceq=0withx=[ΔLt1,ΔLtf,rt1,rtf,βa,βp].This problem can be solved with a gradient-based optimisation algorithm like MATLAB’sfmincon. Note that the time of flight TOF¯ is specified a priori and, therefore, it might occur that this duration is too short as to obtain the change in the orbital parameters specified by the boundary constraints. In this case, the problem is infeasible and the optimisation is terminated after a maximum of 50 iterations if the constraints are not satisfied.In the following, an example of transfer from an elliptical orbit with 300 km perigee altitude and eccentricity 0.031 (corresponding to the final orbit of a de-orbiting strategy) to a circular orbit of 1100 km altitude (corresponding to the orbit of the next debris in an hypothetical removal sequence). Parameters of the two orbits are reported in Table2. Note that the total plane rotation Δi in this case is 10 degrees. The specified time of flight is 70 days.Table 2
Parameters of departure and arrival orbits.
a (km)ei (deg)Departure6892.240.0310Arrival7478.16010First, it is considered the case of a coplanar transfer, that is,Δi=0 will be computed. The optimisation problem was solved with fmincon in 6 iterations and less than 10 seconds, returning a minimum ΔV cost of 0.301 km/s, with 1001 revolutions. Figures 7(a)–7(c) report, respectively, the variation of semimajor axis, eccentricity, apogee and perigee radii. One can see that a is monotonically increasing while e on the other hand is monotonically decreasing to zero. In order to reach the desired circular orbit, the perigee had to be raised by almost 700 km while the apogee had to be raised by some 400 km. This higher effort needed to raise the perigee explains the larger amplitude of apogee thrusting arcs ΔLa compared to perigee ones ΔLp (as shown in Figure 8(a)). The azimuth thrust angles αp,αa (see Figure 8(b)) are both positive since both the perigee and apogee are raised. βp and βa are obviously zero because the transfer is coplanar and thus Δi is constantly null.(a) Variation of semi-major axis, (b) eccentricity, (c) perigee, and apogee radiuses for multi-revolution orbital transfer (coplanar case).
(a)(b)(c)Control parameters for multi-revolution orbital transfer (coplanar case): (a) thrust arc length; (b) azimuth and elevation.
(a)(b)The same problem, but this time with the 10° plane change specified in Table2, returns a ΔV of 1.480 km/s with 1004 revolutions. The high cost of out plane manoeuvres is well exemplified by the fact that the ΔV required is more than four times larger than a coplanar transfer. As can be seen in Figures 9(a), 9(b), and 9(d), semimajor axis, eccentricity, apogee, and perigee radii show a similar behaviour to the coplanar case while this time also the inclination (as in Figure 9(c)) increases monotonically to 10 degrees. By analysing the control parameters in Figure 10(a) one can see that this time the amplitude of the perigee arcs is in general larger than the apogee ones, even if, like in the coplanar case, the increase in perigee is much larger than that of the apogee. This fact is explained by the fact that the out-of-plane component at perigee βp is close to 90° (see Figure 10(b)), meaning that the thrusting action at perigee is mostly devoted to the plane change. In contrast, βa is smaller in magnitude, around −70° (the opposite sign is due to the fact that it is advantageous to invert the out-of-plane component twice per revolution), therefore, with a higher in plane component devoted to perigee raising.(a) variation of semi-major axis, (b) eccentricity, (c) plane change, (d) perigee, and apogee radiuses for multirevolution orbital transfer (10° plane change).
(a)(b)(c)(d)Control parameters for multirevolution orbital transfer (10° plane change): (a) thrust arc length; (b) azimuth and elevation.
(a)(b)


## 3.1. De-Orbiting Trajectory Model

The objective of the de-orbiting phase is that of lowering an initial circular orbit such that its perigee is equal or below 300 km, which basically translates into a perigee lowering manoeuvre. Therefore, it is appropriate to assume that in general, as soon as the initial circular orbit becomes slightly eccentric, one keeps thrusting around the apogee in order to lower the perigee. The thrust level will also be kept at its maximum in order to minimize gravity losses. Moreover, since the de-orbit condition is independent of the final orbit’s orientation, one can reasonably assume that the perigee lowering will be performed in-plane. In this sense, the only Keplerian parameters which need to be altered are the semimajor axis and eccentricity. By analysing the structure of Gauss’ variational equations:(6)dadt=2a2h(esin⁡θar+praθ),dedt=1h{psin⁡θar+[(p+r)cos⁡θ+re]aθ},didt=rcos⁡θhah,dΩdt=rsin⁡θhsin⁡iah,dωdt=1he[-pcos⁡θar+(p+r)sin⁡θaθ]-rsin⁡θcos⁡ihsin⁡iah,dθdt=hr2+1eh[pcos⁡θar-(p+r)sin⁡θaθ],
where a is the semimajor axis, e the eccentricity, i the inclination Ω is the Right Ascension of the Ascending node (RAAN), ω the argument of perigee, θ the true anomaly, p the semilatus rectus, and h the angular momentum; ar,aθ, and ah are the radial, transversal, and out-of-plane components of the thrust acceleration. If one considers the case of thrusting with maximum acceleration along arcs which are symmetrical around apogee (θ=π), one can see that the contributions to semimajor axis and eccentricity variations given by the ar components are negligible (since they are multiplied by sinθ). Therefore, a good suboptimal thrust direction can be obtained by imposing aθ as the only non-zero component of the thrust acceleration. Under these assumptions, the variation of Keplerian parameters will be [20](7)dadt=2a2hpraθ,dedt=1h[(p+r)cos⁡θ+re]aθ,didt=0,dΩdt=0,dωdt=(p+r)hesin⁡θaθ,dθdt=hr2+(p+r)sin⁡θaθ.
It should be noted that the terms of the variation of ω and θ which depend on aθ will also be very small due to the presence of sinθ integrated around θ=π. Figure 2 visualises the proposed pattern for thrusting arcs.Figure 2
Thrusting arc around apogee with thrust directed along transverse direction.In order to obtain a fast propagation of the thrusting arcs, the analytical propagation of perturbed motion with finite perturbative elements in time (FPET) derived in [19, 21] will be used. In order to employ FPET, one has also to assume that the thrust acceleration is constant around each thrusting arc, which is reasonable given the low propellant consumption per arc. This assumption ensures that, if the engine thrust is constant, the resulting acceleration can also be considered constant over short thrusting arcs.Motion propagation with FPET is based on a first-order analytical solution of perturbed Keplerian motion. In this formulation, the state is expressed in nonsingular equinoctial elements:(8)a,P1=esin⁡(Ω+ω),P2=ecos⁡(Ω+ω),Q1=tani2sin⁡Ω,Q2=tani2cos⁡Ω,L=Ω+ω+θ.
Assuming constant thrust-acceleration in the radial-transverse reference frame:
(9)a=[araθah]T=ε[cos⁡⁡αcos⁡⁡βsin⁡⁡αcos⁡⁡βsin⁡⁡β]T,
then one can obtain a first-order analytical expansion of the variation of equinoctial elements, parameterised in Longitude L and with respect to a reference longitude L0:
(10)a(L)=a0(L0)+εa1(L0,ΔL,α,β),P1(L)=P10(L0)+εP11(L0,ΔL,α,β),P2(L)=P20(L0)+εP21(L0,ΔL,α,β),Q1(L)=Q10(L0)+εQ11(L0,ΔL,α,β),Q2(L)=Q20(L0)+εQ21(L0,ΔL,α,β),t(L)=t00(L0,ΔL)+εt11(L0,ΔL,α,β),
where(11)L=L0+ΔL,
where a a0P10,P20,Q10,Q20 are reference values at L0 and a1,P11,P21,Q11,Q21 are first-order terms as reported in [19, 21]. In [19],it has also been shown that this analytical propagation scheme provides good accuracy along relatively long trajectory arcs.As explained above, the only nonzero component of the acceleration will beaθ and since the aim is obtaining a decrease of the orbit energy, it will also be in the negative direction. Therefore, the acceleration azimuth will be α=-π/2 and the elevation β=0 (since, as already mentioned, the motion will be within the initial orbit plane). The variation of equinoctial elements after an apogee thrusting arc will be given by
(12)E+=[aP1P2]L+=La+ΔLa=[aP1P2]L-=La-ΔLa+ε[a1(2ΔLa,π2,0)P11(2ΔLa,π2,0)P21(2ΔLa,π2,0)]=E-+εf(2ΔLa,π2,0),
where La is the apocentre longitude, L- and L+ are the longitudes at the start and end of thrusting, respectively. ΔLa is the semiamplitude of the apogee thrusting arc. Note that, given that β=0, there is no variation on Q1 and Q2 and thus they are omitted. The coasting time is computed from the last of (10) as(13)tthrust=t00(La-ΔLa,2ΔLa)+εt11(2ΔLa,-π2,0).Since the thrust magnitude and direction are fixed, the only free control parameter is the semiamplitudeΔLa for each orbit. In order to keep the number of decision variables to a minimum, the semiamplitude for each orbit is computed from a piece-wise linear polynomial interpolating a limited number of ΔLa,i over a number of orbits. The nodes ΔLa,i are equally distributed between orbit 1 and an arbitrary number of orbits (in this paper, 1200 was found to be adequate). In this paper, the number of interpolating nodes was limited to 2: ΔLa1 and ΔLaf.In order to evaluate the time and ΔV needed to de-orbit a piece of debris from its initial orbit with semimajor axis adebr0, given a set of decision (or control) parameters ΔLa1 and ΔLaf, the following procedure was implemented.(1)
Compute the set of initial Equinoctial parametersL0 and E0=[adebr0P10P20]T, where P10 and P20 will be null due to the fact that the initial orbit is circular.(2)
Initialise the number of orbits, the totalΔV and time of flight to zero:
(14)Norbit=0,ΔV=0,TOF=0.(3)
SetE-=E0 and Lcoast=L0.(4)
Initialise the mass of the IBS spacecraft:(15)mIBS=mIBS0.(5)
WhileNorbit is smaller than Norbitsmax,
(a)
Norbit=Norbit+1.(b)
Interpolate the amplitude of the thrusting arc in the current orbit, that is,ΔLa and compute L-=La-ΔLa and L+=La+ΔLa.(c)
Compute the accelerationεIBS-debr acting on the IBS-debris combination from (5).(d)
Compute the time of flighttcoast spent coasting from Lcoast to L-.(e)
Compute the equinoctial parameters after the thrusting arcE+ as in (12).(f)
Compute the current perigee radiusrp and if this is lower than the threshold r̅p=300 km proceed to step 6, otherwise, proceed to step g.(g)
Compute the thrusting timetthrust from (13) and update the total ΔV cost:
(16)ΔV=ΔV+εIBS-debrtthrust.(h)
Update the total time of flight:(17)TOF=TOF+tcoast+tthrust.(i)
Update the IBS spacecraft mass:(18)mIBS=(mIBS+2mdebr)exp⁡(-εIBS-debrtthrustIspg0)-2mdebr.(j)
SetE-=E+ and Lcoast=L+.(6)
Back track the value of the longitudeLf for which rp=r̅p and compute the related and tthrust from (13) and update TOF and ΔV accordingly. Compute the equinoctial parameters Ef at Lf from (12).At this point, one gets theΔV, the time of flight TOF, and the semimajor axis and eccentricity of the final orbit (which are easily computed from Ef). It is important to note that, given the simplifications introduced, once one sets the initial mass and orbit of the piece of debris, and the characteristics of the IBS propulsion system, that is, Ftot and Isp, the de-orbit depends exclusively on the mass of the IBS mIBS0 at the beginning of the de-orbit phase and the interpolating values for ΔLa, that is, ΔLa1 and ΔLaf. Therefore, it was decided to precompute the corresponding ΔV and TOF for a given set of these three parameters and for each piece of debris (i.e., for each mdebr and adebr0). Table 1 reports upper and lower bounds for mIBS0, ΔLa1 and ΔLaf, and the number of samples taken, equally distributed.Table 1
Bounds and number of samples for the de-orbit parameters.
mIBS0ΔLa1ΔLafLower boundMdry+100=350 kg00Upper boundMlaunch=1000 kgππSamples85050Given the limited number of decision variables, for each piece of debris, one has 20000 de-orbit instances to propagate. Since each instance requires typically 1·10−2 s of CPU time, with a code implemented in MATLAB and running on a 3.16 GHz, 4 GB desktop PC running Windows 7, the whole computation can be completed in roughly five minutes. The set of de-orbit ΔV and TOF is then used to build a response surface, or surrogate model, of the de-orbiting process. Figures 3(a) and 3(b) show examples of two-dimensional surface, respectively, for ΔV and TOF, with respect to a fixed mIBS0 of 300 kg. One can see that the two quantities show opposite trends, the ΔV being high when the TOF is low and vice versa. Figures 4(a) and 4(b) show the final semimajor axis and eccentricity, respectively. Note that the minimum TOF transfer corresponds to a quasicircular spiralling trajectory in which the IBS spacecraft is thrusting continuously. On the other hand, the minimum ΔV transfer corresponds also to the one with maximum final eccentricity.(a)ΔV and (b) TOF surfaces with respect to ΔLa1 and ΔLaf for mIBS0= 300 kg, adebr0 = 7128 km and mdebr = 120 kg.
(a)(b)(a) final semi-major axis and (b) eccentricity after de-orbit with respect toΔLa1 and ΔLaf for mIBS0 = 300 kg, adebr0 = 7128 km and mdebr = 120 kg.
(a)(b)Now, it is desirable that the surrogate model returns theΔV cost as a function of mIBS0, mdebr, adebr0, and TOF. From the available data relating the ΔV and TOF to the decision variables ΔLa1 and ΔLaf, one can derive the functional relationship between ΔV and TOF. Given a triplet mIBS0,mdebr,adebr0, each TOF value defines a level curve on the ΔLa1 and ΔLaf plane (see Figure 3(a)), which can be mapped into a set of ΔV values (see Figure 3(b)). Within this set, one can take the element with minimum ΔV. Thus, for each time of flight, between a minimum and a maximum, one can derive the corresponding minimum ΔV cost. A similar procedure is followed to find the functional relationship between the final semimajor axis and the TOF. Note that there is no need to do the same for the eccentricity given the fact that the final perigee radius is fixed at r̅p and, therefore, the final e can be computed from the final a. In this way one can build the two surrogate models: (19)ΔV=fΔV,interp(TOF,mIBS0,mdebr,adebr0),af=faf,interp(TOF,mIBS0,mdebr,adebr0).
Figures 5(a) and 5(b) show examples of tridimensional plots (mIBS0-TOF-ΔV and mIBS0-TOF-af, resp.) created by evaluating the surrogated models keeping adebr0 and mdebr fixed. In Figure 5(a), one can see that there is a large plateau region corresponding to large time of flights and a smaller region close to the minimum TOF where the de-orbit cost increases very steeply and the final semimajor axis in Figure 5(b) similarly decreases. The complete procedure for the creation of the interpolated de-orbit cost models requires few minutes of CPU time and once completed allows for a very fast estimation of the de-orbit cost. The surrogated models will be extremely useful in the multiobjective optimisation of debris removal sequences as it will be shown in the following sections.3D Plot of surrogate models foradebr0=7128 km and mdebr=120 kg: (a) ΔV; (b) af.
(a)(b)


## 3.2. Orbit Transfer Model

According to the scenario presented in Section3, after having left the debris on a re-entry orbit, the IBS will have to transfer to the orbit of the next debris and rendezvous with it. The design of such a transfer arc would normally require the solution of a fixed-time two point boundary value problem (TPBVP) which would be computationally very expensive given the high number of control parameters and constraints involved. A second simplified model was then created to quickly estimate the cost of a low-thrust multirevolution orbit transfer with boundary constraints. The approach and assumptions presented in this section are similar to those already introduced for the de-orbit model.First, given the limited acceleration provided by low-thrust propulsion systems, one should consider that the orbit transfer will require a high number of multiple revolutions around the Earth, typical in the range of hundreds to few thousands. In this sense, it is possible to argue that achieving the proper phasing to transfer from the initial to final orbit would not be a major issue. Even a small variation ofω and θ per revolution would be sufficient to attain the required orientation to rendezvous with the piece of debris. Moreover, it is important to bear in mind that, in order to de-orbit the previous debris in the sequence, the IBS spacecraft, started from a circular orbit which was subsequently modified into an elliptical one with perigee r̅p. Thus, it would be also possible to conveniently adjust the start point of the de-orbit procedure from the circular orbit in order to obtain the proper phasing once this is completed. For all these reasons, it is assumed that in this particular case, the phasing problem will have a negligible effect on the ΔV and time required to rendezvous with the next piece of debris in the sequence. Therefore, in the following, it is assumed that it is not necessary to match the arrival ω and θ computed with the simplified model with those of the target object. Matching the target inclination i and RAAN Ω, instead, cannot be ignored without introducing a considerable error in the ΔV cost. In order to match the inclination and RAAN difference, one need to take into account only the geometric angle between planes of the initial and final orbits, which is given by
(20)Δi=arccos⁡(-cos⁡i0cos⁡(π-if)+sin⁡i0sin⁡(π-if)cos⁡(Ωf-Ω0)).Thus, in order to account forΔi, the inclination of the initial orbit is fictitiously set to zero, while the final one is set at Δi. The matching of the RAAN is assured by performing the circularisation properly. The assumption is that the deorbiting of one piece of debris starts at a true anomaly such that the resulting elliptical orbit has the line of apses perpendicular with the line of the nodes of the following piece of debris. Since the orbits of the debris are assumed to be circular, it is always possible to start the deorbiting at the right true anomaly with minimum delay. This hypothesis will be discussed in more detail with some numerical examples in Section 4.With these assumptions, the main issue in designing the multirevolution transfer will be that of achieving the required change in the apogee and perigee radii in order to match those of the final orbit, and to achieve the required rotation of the orbit plane.The control of eccentricity and semimajor axis, required to match the target perigee and apogee altitudes, can be obtained by inserting two thrust arcs per revolution, one around the apogee and one around the perigee. This methodology is analogous to what was done in the previous section for the perigee lowering. In the same way, the radial component of the thrust acceleration is set to zero. The transverse component this time can have either positive or negative sign (α=±π/2) depending whether the perigee (or apogee) needs to be raised or lowered.Since a plane change is required, the out-of-plane component of the thrust acceleration is nonzero. Thanks to this, the control parameters can be reduced to the semiamplitude of the apogee and perigee thrusting arcs,ΔLa and ΔLp, the sign of the θ component of the thrust acceleration (i.e., the sign of αa,αp=±π/2), and the out-of-plane component in the same arcs, βa and βp. Define ΔLthrust as half the total thrusting arc length and rt as the ratio of ΔLthrust which is devoted to apogee thrusting. In order to have a parameterisation which accounts also for the sign of αa and αp, the following one is proposed:
(21)αa={π2ΔLt≥0,-π2ΔLt<0,αp={αa0≤rt≤1,-αa1<rt≤2,ΔLa={rt|ΔLthrust|0≤rt≤1,(2-rt)|ΔLthrust|1<rt≤2,ΔLp=|ΔLthrust|-ΔLa,
with
(22)ΔLthrust∈[-ππ],rt∈[02].To define the actual values ofΔLthrust and rt in each revolution, an interpolating strategy from a set of nodal values similar to the one used for the de-orbit model is adopted. In this case, however, the interpolated values will not be computed with respect to the current revolution number but with respect to the time. Again the number of interpolating nodes can be chosen arbitrarily and is set to 2 in this case, ΔLt1,ΔLtf,rt1,rtf. For βa and βp, it is chosen to have a constant value along the entire transfer. The thrusting pattern along each revolution is shown in Figure 6.Figure 6
Thrusting arcs around apogee and perigee.Given a set of control parameters[ΔLt1ΔLtfrt1rtfβaβp], a multirevolution transfer with specified duration TOF¯, departing from an orbit defined by [a0e00]T and targeted to an orbit defined by [afefΔi]T, is propagated according to the following procedure.(1)
Compute the set of initial Equinoctial parametersL0 and E0=[a0P10P20Q10Q20]T. Q10 and Q20 will be zero since the initial inclination is arbitrarily set to zero.(2)
Compute the set of target Equinoctial parametersE̅f=[a̅fP̅1fP̅2fQ̅1fQ̅2f]T. Note that P̅1f and P̅2f will be zero since in this case the target orbit is a circular one.(3)
Initialise the total ΔV and Time of flight to zero:
(23)ΔV=0,TOF=0.(4)
SetEp-=E0 and Lcoast,a=L0.(5)
Initialise the mass of the IBS spacecraft:(24)mIBS=mIBS0.(6)
WhileTOF<TOF¯.
(a)
Compute the interpolated values forΔLt and rt. Hence, calculate αa,αp,ΔLa, and ΔLp from (21).(b)
Compute(25)La-=La-ΔLaLa+=La+ΔLa,Lp-=Lp-ΔLpLp+=Lp+ΔLp.(c)
Compute the current acceleration acting on the spacecraft:(26)εIBS=FtotmIBS.(d)
Compute the time of flighttcoast,p spent coasting before perigee from Lcoast,p to Lp-.(e)
Compute the Equinoctial parameters after the thrusting perigee arcEp+ with an expression analogous to (12).(f)
Compute the thrusting time at perigeetthrust,p from (13). If (TOF¯-TOF)<tthrust,p, proceed to step g. Otherwise, break the iterative sequence and go to step 7.(g)
UpdateΔV and TOF:
(27)ΔV=ΔV+εIBStthrust,p,TOF=TOF+tcoast,p+tthrust,p.(h)
Update the IBS spacecraft mass:(28)mIBS=mIBSexp⁡(-εIBStthrust,pIspg0).(i)
SetEa-=Ep+ and Lcoast,a=Lp+.(j)
Compute the current acceleration on the spacecraft:(29)εIBS=FtotmIBS.(k)
Compute the time of flighttcoast,a spent coasting before apogee from Lcoast,a to La-.(l)
Compute the equinoctial parameters after the thrusting apogee arcEa+ as in (12).(m)
Compute the thrusting time at apogeetthrust,a from (13). If (TOF¯-TOF)<tthrust,a, proceed to step n. Otherwise, break the iterative sequence and go to step 7.(n)
UpdateΔV and TOF:
(30)ΔV=ΔV+εIBStthrust,a,TOF=TOF+tcoast,a+tthrust,a.(o)
Update the IBS spacecraft mass:(31)mIBS=mIBSexp⁡(-εIBStthrust,aIspg0).(p)
SetEp-=Ea+ and Lcoast,p=La+.(7)
Back-track the point at whichTOF=TOF¯ and compute the corresponding equinoctial parameters Ef=[afP1fP2fQ1fQ2f]T and update ΔV accordingly.(8)
Compute the mismatch between the actual final conditions and the target orbit:(32)Ceq=[af-a̅fef-e̅f=P1f2+P2f2-P̅1f2+P̅2f2if-i̅f=2(arctanQ1f2+Q2f2-arctanQ̅1f2+Q̅2f2)].Summarizing, the TPBVP has been reduced to an optimisation problem in the form:(33)min⁡xΔVs.t.Ceq=0withx=[ΔLt1,ΔLtf,rt1,rtf,βa,βp].This problem can be solved with a gradient-based optimisation algorithm like MATLAB’sfmincon. Note that the time of flight TOF¯ is specified a priori and, therefore, it might occur that this duration is too short as to obtain the change in the orbital parameters specified by the boundary constraints. In this case, the problem is infeasible and the optimisation is terminated after a maximum of 50 iterations if the constraints are not satisfied.In the following, an example of transfer from an elliptical orbit with 300 km perigee altitude and eccentricity 0.031 (corresponding to the final orbit of a de-orbiting strategy) to a circular orbit of 1100 km altitude (corresponding to the orbit of the next debris in an hypothetical removal sequence). Parameters of the two orbits are reported in Table2. Note that the total plane rotation Δi in this case is 10 degrees. The specified time of flight is 70 days.Table 2
Parameters of departure and arrival orbits.
a (km)ei (deg)Departure6892.240.0310Arrival7478.16010First, it is considered the case of a coplanar transfer, that is,Δi=0 will be computed. The optimisation problem was solved with fmincon in 6 iterations and less than 10 seconds, returning a minimum ΔV cost of 0.301 km/s, with 1001 revolutions. Figures 7(a)–7(c) report, respectively, the variation of semimajor axis, eccentricity, apogee and perigee radii. One can see that a is monotonically increasing while e on the other hand is monotonically decreasing to zero. In order to reach the desired circular orbit, the perigee had to be raised by almost 700 km while the apogee had to be raised by some 400 km. This higher effort needed to raise the perigee explains the larger amplitude of apogee thrusting arcs ΔLa compared to perigee ones ΔLp (as shown in Figure 8(a)). The azimuth thrust angles αp,αa (see Figure 8(b)) are both positive since both the perigee and apogee are raised. βp and βa are obviously zero because the transfer is coplanar and thus Δi is constantly null.(a) Variation of semi-major axis, (b) eccentricity, (c) perigee, and apogee radiuses for multi-revolution orbital transfer (coplanar case).
(a)(b)(c)Control parameters for multi-revolution orbital transfer (coplanar case): (a) thrust arc length; (b) azimuth and elevation.
(a)(b)The same problem, but this time with the 10° plane change specified in Table2, returns a ΔV of 1.480 km/s with 1004 revolutions. The high cost of out plane manoeuvres is well exemplified by the fact that the ΔV required is more than four times larger than a coplanar transfer. As can be seen in Figures 9(a), 9(b), and 9(d), semimajor axis, eccentricity, apogee, and perigee radii show a similar behaviour to the coplanar case while this time also the inclination (as in Figure 9(c)) increases monotonically to 10 degrees. By analysing the control parameters in Figure 10(a) one can see that this time the amplitude of the perigee arcs is in general larger than the apogee ones, even if, like in the coplanar case, the increase in perigee is much larger than that of the apogee. This fact is explained by the fact that the out-of-plane component at perigee βp is close to 90° (see Figure 10(b)), meaning that the thrusting action at perigee is mostly devoted to the plane change. In contrast, βa is smaller in magnitude, around −70° (the opposite sign is due to the fact that it is advantageous to invert the out-of-plane component twice per revolution), therefore, with a higher in plane component devoted to perigee raising.(a) variation of semi-major axis, (b) eccentricity, (c) plane change, (d) perigee, and apogee radiuses for multirevolution orbital transfer (10° plane change).
(a)(b)(c)(d)Control parameters for multirevolution orbital transfer (10° plane change): (a) thrust arc length; (b) azimuth and elevation.
(a)(b)


## 4. Multiobjective Optimisation

The aim is now that of optimising the timing and sequence of a removal mission by means of a single IBS spacecraft. It is assumed that the spacecraft departs from an LEO with a 250 km semi-major axis altitude and coplanar with respect to the first piece of debris in the sequence. The five target objects have the orbital parameters and mass reported in Table3. The mass and orbital parameters have been chosen arbitrarily while adhering to the observations in [9, 12] that the most dangerous debris are located in LEO and generally weigh a few hundred kilos. Different values for i and Ω are also taken in order to consider the fact that the pieces of debris, in principle, will be orbiting on different planes. Note that TDo,min has been computed with the procedure detailed in Section 3.1, and therefore, depends on the characteristics of the IBS spacecraft. Moreover, it is also important to remark that these are only best case figures values which were computed with a minimum hypothetical wet mass of 350 kg (much lower than the actual launch mass of 1000 kg). The surrogate models in (19) can in general consider wet masses between 350 kg and 1000 kg, as shown, for example in Figures 5(a) and 5(b).Table 3
Mass, initial orbit parameters, and minimum de-orbit time of the debris.
Debris nr.mass (kg)a (km)ei (deg)Ω (deg)TDO,min (days)15006828.1601652.6721207128.16021503.3633006978.160−22003.6844007478.160−19011.1258007178.16004512.25Table4 reports the relative inclination change between the orbit planes of the 5 different objects, as computed from (20).Table 4
Relative inclination change|Δi| (deg) between orbit planes of the debris.
Debris number234512.161.471.9512—3.632.6523——2.5224———1The de-orbit sequence is defined by the order according to which the five pieces of debris are removed, the time needed to rendezvous withTRV and the time to de-orbit TDO each of them. The order is defined by the integer vector: (34)ord=[i1i2i3i4i5],
which collects the indexes of the objects in the a single debris removal sequence. Since there are five objects, there are 120 possible de-orbit sequences. The other parameters are contained in the vector x:(35)x=[TRV,i1TDO,i1TRV,i2TDO,i2TRV,i3TDO,i3TRV,i4TDO,i4TRV,i5TDO,i5],The performance of each sequence is assessed according to its total ΔVTot cost and time of flight TOFTot. The latter is computed simply as(36)TOFTot=∑x,
The total ΔV cost is calculated sequentially by adding up the costs of each of the ten phases (rendezsvous and de-orbit for each debris). In particular, the cost of the rendezvous ΔVRV is computed by solving the optimisation problem (33) and the de-orbit cost ΔVDO is calculated from the surrogated model in (19). The final conditions after de-orbit are also computed from (19) since they will be the departure conditions for the following rendezvous step. The propellant mass consumption is also taken into account and updated throughout the entire sequence computation. In order to have only a real-valued optimisation problem, it is chosen here to treat each of the 120 sequences as a biobjective optimisation problem with ord fixed and ten design variables defined in x. Therefore, optimisation problem becomes (37)min⁡x∈D[TOFTot(x)ΔVTot(x)].The domainD is defined by the upper and lower boundaries defined in Table 5. Note that the lower boundaries for de-orbit time are set according to the sequence and the minimum times reported in Table 3.Table 5
Optimisation boundaries.
ParameterTRV,i1TDO,i1TRV,i2TDO,i2TRV,i3TDO,i3TRV,i4TDO,i4TRV,i5TDO,i5Lower Bound5TDO,min,i15TDO,min,i25TDO,min,i35TDO,min,i45TDO,min,i5Upper Bound1005010050100501005010050Each biobjective optimisation problem is solved with MACS, a hybrid-memetic stochastic optimisation algorithm [22]. MACS was run for 40000 function evaluations with 30 agents. Each of the 120 optimisation instances required roughly 6 days of computational time to complete. The outputs are represented by the Pareto optimal solutions w.r.t. ΔVTot and ToFTot. Figure 11 to Figure 15 collect the Pareto fronts according to the number of the first object in the sequence, that is, the first index in the vector ord, as introduced in (34). In each figure, each colour represents the Pareto front corresponding to one of the 24 debris removal sequences starting with the same object. For example, Figure 11 includes the Pareto fronts of sequences 12345, 13245, 14235, 15234, 12435, and so forth.Figure 11
Pareto fronts for sequences starting with debris nr. 1.Figure 12
Pareto fronts for sequences starting with debris nr. 2.Figure 13
Pareto fronts for sequences starting with debris nr. 3.Figure 14
Pareto fronts for sequences starting with debris nr. 4.Figure 15
Pareto fronts for sequences starting with debris nr. 5.From a visual inspection of the fronts, it is possible to see that sequences starting from debris number 1 seem to present the bestΔVTot-TOFTot combination, since for most of them the ΔV cost is comprised between 2 and 2.5 km/s. The corresponding times of flight are comprised roughly between 100 and 500 days. The sequences starting with debris nr. 3 and nr. 2 also have a good ΔV while those starting with nr. 4 and nr. 5 appear to be worst. By combining all the partial Pareto fronts, one obtains the globally optimal solutions, as reported in Figure 16.Figure 16
Global Pareto front.One can see that the global Pareto front is composed by individual solutions belonging exclusively to sequence 13452, which is, therefore, globally dominant. In order to rank the degree of optimality of each sequence, it is proposed to use an approach inspired by the performance metrics for optimisation algorithms proposed in [23]. Define PFg as the set of the points of the globally optimal Pareto front while PFord is the set of points belonging to the Pareto front corresponding to sequence ord. Define then the ranking parameter of sequence ord as (38)Conv(ord)=100∑i=1Nordmin⁡gj∈PFg⁡((fi-gj)/δ)Nordfi∈PFord,∀i=1,…,Nord.Conv is given by averaging the distance of each point ofPFord from the closest point of PFg. The closest PFord is to PFg and the lower Conv will be. Table 6 reports the ranking of the sequences according to Conv.Table 6
Ranking of the de-orbit sequences.
RankordConv(ord)RankordConv(ord)RankordConv(ord)1134520414251321.18815214331.432135425.14421523421.26823214531.623135246.61433245121.46835412331.724124536.78445213421.46845413231.835125437.25453452121.52854213532.436315429.41463514221.79865231432.707314529.85473521421.99874253133.0583451211.59483425122.02882143533.9692451312.15495243122.04895423134.05101524312.16504513222.13902314534.31111253412.33515431223.39912351434.56123125412.37522154323.60925342134.67131543213.24532431523.62932534134.71143512413.87544135223.81941432534.91151325414.22554315223.90954125335.20163152414.36561243524.40963251435.42171534214.48573412524.53971423535.65181342516.30581532424.89983254136.42192453116.53595314224.90995123436.81201452316.65602315425.611004215336.91211435216.69615312425.671015142338.10223415217.16625124325.801025432138.25232513417.17634351225.831034523140.16241235417.47643142525.951045143240.98251425317.63651234525.961054152341.91263124517.81662145326.011064532144.72271542317.85675241326.091073241545.05285134217.88685132426.561084235145.38291453218.05693524126.681094352145.43302541318.07702514326.771105321445.72315421318.17712415326.931114325145.87323541218.48723421527.521122354146.11332134519.32732153428.001135234146.89342543119.43743215429.651144132547.14351324519.55754312530.171154153247.50363542119.56762345130.401165324148.31372531419.97772435130.971172341548.84384521319.98784512331.071184231548.85394531220.07792413531.181194321552.91402135420.07805341231.221204123565.42As one would expect, sequence 13452 has the lowest Conv since it coincides with part of the global Pareto front. Sequences 13524, 13542 and 12543 have also a low Conv index and thus they are quite close to the globally optimal solution, as shown in Figure17. In general, as already noted before, there is a strong dependence of the quality of the sequence from its first element. One can see that the first ranks are occupied mostly by sequences starting with debris nr. 1 and 3, while those with nr. 4 and 5 have highest Conv and, therefore, occupy predominantly the worst ones. Those starting with nr. 2 are somewhat in the middle. The fact that solutions with nr. 1 and 3 are privileged as first elements in the sequence might be explained from the fact that they lie in the two lowest orbits (see Table 3), and therefore, are easier to reach (please keep in mind that for the rendezvous with the first debris there is no plane change since it is assumed to depart from a coplanar orbit). Another interesting observation is that the best sequences tend to avoid the largest plane changes. For example, in 13452 the plane changes are 1.47°, 2.52°, 1°, and 2°. On the contrary, in the worst one according to Conv, that is, 32415, they are 3.63°, 2.65°, 1.95°, and 1°.Figure 17
Pareto fronts corresponding to the four best sequences according to Conv.Table7 reports the minimum values for the performance parameters associated to each sequence, that is, the extreme points of the Pareto fronts. Similar considerations to those made previously also apply to this case, with best values given by sequences starting with nr. 1 and 3 and the worst ones with nr. 4 and 5.Table 7
BestΔVTot and ToFTot for each sequence. Best values are in bold. Worst values are underlined.
ordmin(ΔVTot) (km/s)min(ToFTot) (days)ordmin(ΔVTot) (km/s)min(ToFTot) (days)ordmin(ΔVTot) (km/s)min(TOFTot) (days)123452.30108.17245132.17106.64423152.53116.99123542.26107.81245312.26105.94423512.63114.24124352.27106.25251342.24102.96425132.28104.48124532.13100.63251432.41109.53425312.36109.50125342.18105.81253142.30107.66431252.42109.26125432.13103.26253412.49107.22431522.34107.36132452.22102.72254132.26104.36432152.67116.73132542.11103.03254312.27107.97432512.56113.93134252.15103.12312452.20100.73435122.43108.12134521.9896.35312542.10103.46435212.53111.63135242.07101.03314252.30106.79451232.46106.61135422.02100.08314522.1297.810451322.33102.12142352.45115.10315242.15104.32452132.25102.99142532.21104.45315422.12100.52452312.42110.85143252.42112.87321452.44111.30453122.27101.57143522.21105.30321542.35107.73453212.55111.39145232.25107.07324152.51115.31512342.49110.19145322.27105.43324512.33107.85512432.38107.21152342.29107.13325142.38107.70513242.36106.79152432.14102.56325412.40109.14513422.25103.38153242.27106.71341252.42109.68514232.53113.27153422.17102.45341522.33104.78514322.55112.77154232.26109.45342152.36112.19521342.29106.29154322.24106.63342512.30107.17521432.44108.36213452.31103.98345122.18101.86523142.46116.72213542.24103.07345212.24104.80523412.61112.97214352.58115.15351242.27103.81524132.30106.23214532.38106.26351422.30105.62524312.37108.24215342.40113.76352142.32109.19531242.29103.32215432.32110.97352412.37111.19531422.36108.06231452.47113.45354122.28101.50532142.60114.17231542.36107.94354212.29108.91532412.62116.59234152.63114.64412352.70116.91534122.45106.79234512.48111.27412532.47107.83534212.46112.91235142.55114.91413252.55113.15541232.49112.69235412.54111.22413522.37108.25541322.38105.83241352.42107.59415232.54111.57542132.24104.67241532.43108.50415322.57112.14542312.42115.17243152.38110.73421352.44115.59543122.22107.06243512.42108.28421532.47108.75543212.50116.42Table8 shows details about the best VTot solution, with sequence 13452. Note that, in general, the ΔV cost of each phase is relatively low, thus leading to the minimum total cost of 1.98 km/s. Correspondingly, their duration is long, meaning that slow but more efficient transfers are preferred. This behaviour is also confirmed by the fact that the de-orbit conditions have nonnegligible eccentricities, which means also that the amplitude of the apogee thrusting arcs during de-orbit (see Figure 2) is kept to a minimum. In this way, propellant is devoted to lowering the perigee only with minimum variation of the apogee altitude.Table 8
Debris removal sequence and timing for minimumΔVTot.
PhaseFinal Keplerian elementsDuration (days)ΔV (km/s)mass (kg)a (km)ei (deg)Ω (deg)Departure6628.160.010165——1000Nr. 1 reached6828.16016550.115996.11Nr. 1 de-orbited6752.690.01116522.060.043993.21Nr. 3 reached6978.160−220088.100.239985.17Nr. 3 de-orbited6826.440.022−220025.960.084980.63Nr. 4 reached7478.160−19066.710.476964.88Nr. 4 de-orbited7055.540.053−19034.330.221951.69Nr. 5 reached7178.16004555.890.241943.91Nr. 5 de-orbited6912.180.03404530.770.144931.48Nr. 2 reached7128.160215056.980.297922.12Nr. 2 de-orbited6901.390.032215033.990.124917.24By analysing in more detail the ΔV cost breakdown, one can see, for example, that the highest figures 0.476 km/s are given by the rendezvous with debris nr. 4 from the de-orbit conditions of debris nr. 3. This high value is justified by the fact that reaching the final orbit radius of 7478.16 requires an apogee raise of 501 km from 6977 km and a perigee raise of 800 km from 6678 km. At the same time, there is also a rotation of the orbit plane of 2.52°. By comparison, the rendezvous with nr. 5 after the de-orbit of nr. 4 is comparatively cheaper even if the radius of the target orbit is still high. In this case, the perigee raise is 500 km while the apogee on the other hand needs to be lowered by 252 km from 7430 km since piece of debris nr. 4 is released on a relatively eccentric orbit with e=0.053. Plane rotation in this case is only 1°.Table9 reports details about the minimum ToFTot solution, again with sequence 13452. In contrast to what has been remarked for the previous case, here obviously the duration of each phase is kept to a minimum. For example, values for de-orbit times are very close to the minima reported in Table 3. Conversely, ΔV costs are higher than those in Table 8. Moreover, one can see that the de-orbit trajectories are quasicircular, which suggests that the thrusting arcs are not restricted to apogee passages but cover almost entirely each revolution (i.e., with reference to Figure 2, La≈180°).Table 9
Debris removal sequence and timing for minimum TOFTot.
PhaseFinal Keplerian elementsDuration (days)ΔV (km/s)Mass  (kg)a (km)ei (deg)Ω (deg)Departure6628.160.010165——1000Nr. 1 reached6828.16016550.115996.11Nr. 1 de-orbited6685.240.0011654.040.081990.61Nr. 3 reached6978.160−22008.590.312980.14Nr. 3 de-orbited6701.870.004−22006.290.154971.87Nr. 4 reached7478.160−19014.790.664950.17Nr. 4 de-orbited6789.060.016−19017.130.362928.76Nr. 5 reached7178.1600457.990.281919.92Nr. 5 de-orbited6715.720.00604515.90.252898.39Nr. 2 reached7128.16021509.870.466884.28Nr. 2 de-orbited6725.900.00721506.750.221875.87A final note is devoted to the assumption mentioned in Section3.2 that the delay due to phasing will be relatively negligible compared to the total transfer time. First of all, one has to consider that each de-orbit-rendezvous couplet is actually a transfer between two circular orbits with different altitude, phasing and orbit plane. In this sense, the related transfer strategy first lowers the perigee down to 300 km; then, in the second phase the apogee, and perigee altitudes are adjusted to match those of the target orbit and at the same time the orbit plane is rotated around the line of nodes. In order to obtain a worst-case estimation of the delay, it is chosen to decompose the latter into the contribution determined by the inclination change twait,Δiand the one given by in-plane phasing twait,Δϕ. The former stems from the assumption made in Section 3.2 that the perigee lowering phase from the initial circular orbit is started such that the lines of apses is perpendicular to the line of nodes defined by the intersection of the orbit planes of the current piece of debris and the next one in the sequence. The maximum wait time is obtained when the line of nodes is aligned with the line of apses and is, therefore, given by half the orbit period of the departure circular orbit:
(39)max⁡(twait,Δi)=πn0,
where n0 is the angular velocity of the initial circular orbit. After the line of apses is properly aligned in order to reach the target orbit plane, there remains, however, the problem of in-plane phasing. As a first step, the case of a quasicircular transfer is considered, noting that this is actually the case for minimum ΔV sequences as the one reported in Table 8. If one considers the case of a transfer between two circular coplanar orbits, the phasing of the departure and arrival ones can by expressed as
(40)ϕ0=n0twait,Δϕ+Δϕtransf(ttransf),ϕf=nf(twait,Δϕ+ttransf)+Δϕ0,
where n is the angular velocity of the current orbit, nf is the one of the arrival orbit and is introducing along the transfer in order to match the phase of the arrival orbit. Δϕ0 is the nominal phase difference between the two orbits at time of departure, computed simply from the initial and final argument of perigee and true anomaly:
(41)Δϕ0=(ωf-ω0)+(θf-θ0),ϕ0 and ϕf can differ by multiples of 2π, therefore, by combining (40):
(42)twait,Δϕ(n-nf)=Δϕ0+Δϕtransf+nfttransf+2kπ,k∈Z.One can see that, once the transfer type is defined, the left side of (42) is constant and since k is an arbitrary integer, one can write
(43)twait,Δϕ=ΔϕTot|n-nf|ΔϕTot∈[02π],
and thus the worst-case value for the delay twait is obtained obviously for ΔϕTot=2π. Since we are dealing with a LT transfer in which the semimajor axis is continuously varied, also the angular velocity n at a certain point of the transfer is varying accordingly. Also, since it is assumed that the transfer is quasicircular, one can insert a coasting arc of duration twait,Δϕ at the point in which the ratio 1/|n-nf| (which depends on the radii of the current and target orbits) is at its lowest. This condition typically occurs when the end of the de-orbit phase is reached.If the transfer type is not quasicircular but involves spirals with nonnegligible eccentricity, then an arbitrary delay cannot be introduced without altering the position of the lines of nodes. However, it is still possible to introduce an arbitrary number of coasting arcs of duration equal to the orbital period of the osculating orbit, that is, one full revolution. The phase variation obtained by one such revolution is(44)Δϕ2π(n)=2π|n-nf|n.Note that, given the orbits involved in the transfer,Δϕ2π will be generally a fraction of 2π. If a worst-case phase variation ΔϕTot=2π is to be achieved, the following simple strategy can be used to estimate the corresponding delay: first, k coasting revolutions are performed when the quantity |n-nf| is maximum. In this sense, one can write:
(45)k=⌊Δϕ2π(nk)2π⌋,nk=argmin⁡nΔϕ2π(n).This will bring the phase difference to a quantity which is lower than the maximum phase variation per revolution achievable, leaving a residual phase difference:(46)Δϕres=2π-kΔϕ2π(nk).A last coasting revolution is inserted to delete the residual when the semi-major axis which gives the proper angular velocitynres is reached:
(47)nres=arg⁡(Δϕ2π(n)=Δϕres).The total delay introduced in the worst case is, therefore, given by the sum of the periods of the coasting revolutions:(48)max⁡(twait,Δϕ)=2π(knk+1nres).By applying the above strategies to the minimumΔV and minimum time of flight sequences, we can obtain a worst-case estimation of the additional time introduced by phasing. The maximum delay introduced by the apses alignment in both cases would be 0.14 days. For the minimum time of flight case in Table 9 (i.e., quasicircular sequence), the worst-case delay due to Δϕ is 2.68 days, leading to a total delay of 2.82 days. This value equates to a 2.93% increase compared to the nominal time of flight of 96.35 days, which can be considered acceptable for a preliminary study. On the contrary, in the case of minimum time of flight sequence as in Table 8, the delay due to Δϕ would be 4.58 days, and the total delay 4.72 days, corresponding to a 1.12% increase on the nominal time of flight of 419.79 days. For these reasons, neglecting the phasing appears to be an acceptable approximation in this preliminary study.


## 5. Conclusions

This work presented a novel computational approach for the preliminary design of multispiral trajectories. The approach was applied to the design of an orbit debris removal mission by means of an IBS spacecraft. The models proposed here for the computation of low-thrust many-revolution transfers, allowed for a considerable reduction in control parameters and at the same time for a fast propagation of low-thrust motions thanks to the analytical propagation with FPET. Thanks to the reduced computational cost for the evaluation of a single fetch and deorbit operation, a multiobjective optimisation problem could be solved in which thousands of different debris removal sequences were examined to find the optimal ones with respect toΔV cost and total removal duration. As a result, a considerable number of optimal candidate solutions were found. Analysis of the results showed that the particular removal sequence 13452 is globally optimal. A ranking criterion was proposed to grade all the candidate sequences and identify those that are suboptimal. From the analysis, it was found that there is a dependency of the quality of the sequence on the first target object. Among the open issues for future developments, there is, for example, the possibility of integrating the problem of the sequence choice directly into a single multi-objective optimisation instance, thus obtaining a mixed continuous and discrete optimisation problem. This can be crucial when missions with more than 5–10 debris are considered since the decomposition in fixed-sequence continuous optimisation problems becomes less computationally tractable.Future work will deal with comparing the proposed approach with similar methods like orbit averaging. In addition, although the proposed method has been applied to the special case of a debris removal mission, it is suitable to be extended and applied to more general trajectory design problems which involve many-revolution transfers from elliptical to circular or from elliptical to elliptical orbits. Current developments are incorporating gravity perturbations in the analytical solution to allow the computation of more accurate solutions.



---

*Source: 836250-2012-05-10.xml*