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.
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 the Kessler 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: Section 2 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 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 provided by the other set of ion engines.
Since it is necessary to keep the Shepherd spacecraft at a constant distance from the debris, the thrust shall be such that the second derivative of the distance between the two spacecraft is null: 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: Under the assumption that the total propulsive power of the IBS spacecraft is constant and that the total propulsive thrust is proportional to it , one can write thus Therefore, the maximum acceleration acting on the IBS-debris combination can be computed as a function of the total available thrust : 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-7 km/s2. If one considers instead the spacecraft alone, the acceleration achievable would be slightly higher, 5·10-7 km/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 a priori and optimised one by one.
Each fetch and de-orbit operation is split in two phases:(i)A de-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)A transfer 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: where is the semimajor axis, the eccentricity, the inclination is the Right Ascension of the Ascending node (RAAN), the argument of perigee, the true anomaly, the semilatus rectus, and the angular momentum; , and 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 components are negligible (since they are multiplied by ). Therefore, a good suboptimal thrust direction can be obtained by imposing as the only non-zero component of the thrust acceleration. Under these assumptions, the variation of Keplerian parameters will be [20] It should be noted that the terms of the variation of and which depend on will also be very small due to the presence of integrated around . Figure 2 visualises the proposed pattern for thrusting arcs.
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: Assuming constant thrust-acceleration in the radial-transverse reference frame: then one can obtain a first-order analytical expansion of the variation of equinoctial elements, parameterised in Longitude and with respect to a reference longitude : where where a are reference values at and 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 be 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 and the elevation (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 where is the apocentre longitude, and are the longitudes at the start and end of thrusting, respectively. is the semiamplitude of the apogee thrusting arc. Note that, given that , there is no variation on and and thus they are omitted. The coasting time is computed from the last of (10) as
Since the thrust magnitude and direction are fixed, the only free control parameter is the semiamplitude 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 over a number of orbits. The nodes 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: and .
In order to evaluate the time and ΔV needed to de-orbit a piece of debris from its initial orbit with semimajor axis , given a set of decision (or control) parameters and , the following procedure was implemented.(1)Compute the set of initial Equinoctial parameters and , where and will be null due to the fact that the initial orbit is circular.(2)Initialise the number of orbits, the total and time of flight to zero: (3)Set and .(4)Initialise the mass of the IBS spacecraft: (5)While is smaller than ,(a).(b)Interpolate the amplitude of the thrusting arc in the current orbit, that is, and compute and .(c)Compute the acceleration acting on the IBS-debris combination from (5).(d)Compute the time of flight spent coasting from to .(e)Compute the equinoctial parameters after the thrusting arc as in (12).(f)Compute the current perigee radius and if this is lower than the threshold km proceed to step 6, otherwise, proceed to step g.(g)Compute the thrusting time from (13) and update the total cost: (h)Update the total time of flight: (i)Update the IBS spacecraft mass: (j)Set and .(6)Back track the value of the longitude for which and compute the related and from (13) and update and accordingly. Compute the equinoctial parameters at from (12).
At this point, one gets the , the time of flight , and the semimajor axis and eccentricity of the final orbit (which are easily computed from ). 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, and , the de-orbit depends exclusively on the mass of the IBS at the beginning of the de-orbit phase and the interpolating values for , that is, and . Therefore, it was decided to precompute the corresponding and for a given set of these three parameters and for each piece of debris (i.e., for each and ). Table 1 reports upper and lower bounds for , and , and the number of samples taken, equally distributed.
Given 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 and 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 and , with respect to a fixed of 300 kg. One can see that the two quantities show opposite trends, the being high when the is low and vice versa. Figures 4(a) and 4(b) show the final semimajor axis and eccentricity, respectively. Note that the minimum transfer corresponds to a quasicircular spiralling trajectory in which the IBS spacecraft is thrusting continuously. On the other hand, the minimum transfer corresponds also to the one with maximum final eccentricity.
(a)
(b)
(a)
(b)
Now, it is desirable that the surrogate model returns the cost as a function of , , , and TOF. From the available data relating the and TOF to the decision variables and , one can derive the functional relationship between and TOF. Given a triplet , each TOF value defines a level curve on the and plane (see Figure 3(a)), which can be mapped into a set of values (see Figure 3(b)). Within this set, one can take the element with minimum . Thus, for each time of flight, between a minimum and a maximum, one can derive the corresponding minimum 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 and, therefore, the final can be computed from the final . In this way one can build the two surrogate models: Figures 5(a) and 5(b) show examples of tridimensional plots ( and , resp.) created by evaluating the surrogated models keeping and 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.
(a)
(b)
3.2. Orbit Transfer Model
According to the scenario presented in Section 3, 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 . 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 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 and RAAN , instead, cannot be ignored without introducing a considerable error in the 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
Thus, in order to account for , the inclination of the initial orbit is fictitiously set to zero, while the final one is set at . 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 () 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, and , the sign of the component of the thrust acceleration (i.e., the sign of ), and the out-of-plane component in the same arcs, and . Define as half the total thrusting arc length and as the ratio of which is devoted to apogee thrusting. In order to have a parameterisation which accounts also for the sign of and , the following one is proposed: with
To define the actual values of and 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, . For and , it is chosen to have a constant value along the entire transfer. The thrusting pattern along each revolution is shown in Figure 6.
Given a set of control parameters , a multirevolution transfer with specified duration , departing from an orbit defined by and targeted to an orbit defined by , is propagated according to the following procedure.(1)Compute the set of initial Equinoctial parameters and . and will be zero since the initial inclination is arbitrarily set to zero.(2)Compute the set of target Equinoctial parameters . Note that and 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: (4)Set and .(5)Initialise the mass of the IBS spacecraft: (6)While .(a)Compute the interpolated values for and . Hence, calculate , and from (21).(b)Compute (c)Compute the current acceleration acting on the spacecraft: (d)Compute the time of flight spent coasting before perigee from to .(e)Compute the Equinoctial parameters after the thrusting perigee arc with an expression analogous to (12).(f)Compute the thrusting time at perigee from (13). If , proceed to step g. Otherwise, break the iterative sequence and go to step 7.(g)Update and TOF: (h)Update the IBS spacecraft mass: (i)Set and .(j)Compute the current acceleration on the spacecraft: (k)Compute the time of flight spent coasting before apogee from to .(l)Compute the equinoctial parameters after the thrusting apogee arc as in (12).(m)Compute the thrusting time at apogee from (13). If , proceed to step n. Otherwise, break the iterative sequence and go to step 7.(n)Update and TOF: (o)Update the IBS spacecraft mass: (p)Set and .(7)Back-track the point at which and compute the corresponding equinoctial parameters and update accordingly.(8)Compute the mismatch between the actual final conditions and the target orbit:
Summarizing, the TPBVP has been reduced to an optimisation problem in the form:
This problem can be solved with a gradient-based optimisation algorithm like MATLAB’s fmincon. Note that the time of flight 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 Table 2. Note that the total plane rotation in this case is 10 degrees. The specified time of flight is 70 days.
First, it is considered the case of a coplanar transfer, that is, will be computed. The optimisation problem was solved with fmincon in 6 iterations and less than 10 seconds, returning a minimum 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 is monotonically increasing while 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 compared to perigee ones (as shown in Figure 8(a)). The azimuth thrust angles (see Figure 8(b)) are both positive since both the perigee and apogee are raised. and are obviously zero because the transfer is coplanar and thus is constantly null.
(a)
(b)
(c)
(a)
(b)
The same problem, but this time with the 10° plane change specified in Table 2, returns a of 1.480 km/s with 1004 revolutions. The high cost of out plane manoeuvres is well exemplified by the fact that the 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 is close to 90° (see Figure 10(b)), meaning that the thrusting action at perigee is mostly devoted to the plane change. In contrast, 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)
(b)
(c)
(d)
(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 Table 3. 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 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 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 4 reports the relative inclination change between the orbit planes of the 5 different objects, as computed from (20).
The de-orbit sequence is defined by the order according to which the five pieces of debris are removed, the time needed to rendezvous with and the time to de-orbit each of them. The order is defined by the integer vector: 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 :The performance of each sequence is assessed according to its total cost and time of flight TOFTot. The latter is computed simply as The total 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 is computed by solving the optimisation problem (33) and the de-orbit cost 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 . Therefore, optimisation problem becomes
The domain 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.
Each 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. 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.
From a visual inspection of the fronts, it is possible to see that sequences starting from debris number 1 seem to present the best combination, since for most of them the 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 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.
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 as the set of the points of the globally optimal Pareto front while is the set of points belonging to the Pareto front corresponding to sequence ord. Define then the ranking parameter of sequence ord as
Conv is given by averaging the distance of each point of from the closest point of . The closest is to and the lower Conv will be. Table 6 reports the ranking of the sequences according to Conv.
As 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 Figure 17. 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°.
Table 7 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 8 shows details about the best solution, with sequence 13452. Note that, in general, the 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.
By 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 . Plane rotation in this case is only 1°.
Table 9 reports details about the minimum 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, 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, °).
A final note is devoted to the assumption mentioned in Section 3.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 and the one given by in-plane phasing . 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: where 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 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 where is the angular velocity of the current orbit, is the one of the arrival orbit and is introducing along the transfer in order to match the phase of the arrival orbit. 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: and can differ by multiples of 2π, therefore, by combining (40):
One can see that, once the transfer type is defined, the left side of (42) is constant and since is an arbitrary integer, one can write and thus the worst-case value for the delay is obtained obviously for . Since we are dealing with a LT transfer in which the semimajor axis is continuously varied, also the angular velocity 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 at the point in which the ratio (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
Note that, given the orbits involved in the transfer, will be generally a fraction of . If a worst-case phase variation is to be achieved, the following simple strategy can be used to estimate the corresponding delay: first, coasting revolutions are performed when the quantity is maximum. In this sense, one can write:
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:
A last coasting revolution is inserted to delete the residual when the semi-major axis which gives the proper angular velocity is reached:
The total delay introduced in the worst case is, therefore, given by the sum of the periods of the coasting revolutions:
By applying the above strategies to the minimum 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 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.
Appendix
In the following, the first-order solution for perturbed Keplerian motion is reported (see also [19, 21]). The equinoctial elements at a longitude with respect to a reference longitude are given by a zero-order term plus a first-order term multiplied by the magnitude ε of the perturbing acceleration: where and are, respectively, thrust azimuth and elevation in the radial-transverse reference frame as in Figure 18.
For , and the zero-order term is simply the value at . For the time instead, this is given by the reference time at plus the variation due to unperturbed Keplerian motion, where is an integral in as reported in (A.5). The first-order terms are The terms are integrals in , as Finally (A.6) report the complete analytical expressions for these integrals