Robust Model Predictive Control for an Ion Beam Shepherd in a large-debris removal mission

Contents lists available at ScienceDirect Acta Astronautica journal homepage: www.elsevier.com/locate/actaastro Research paper Robust Model Predictive Control for an Ion Beam Shepherd in a large-debris removal mission Javier Urriosa,∗, Rafael Vazqueza, Francisco Gavilana, Ignacio Alvaradob aDepartment of Aerospace Engineering, Universidad de Sevilla, Sevilla, 41092, Spain bDepartment of Systems Engineering and Automation, Universidad de Sevilla, Sevilla, 41092, Spain ARTICLE INFO Keywords: Active debris removal Ion Beam Shepherd Model Predictive Control Robust control Ray-marching ABSTRACT The increasing accumulation of space debris poses significant risks to spacecraft, making the development of effective debris mitigation technologies necessary. This paper explores the Ion Beam Shepherd (IBS) method as a potential contactless solution for deorbiting large debris objects. The IBS system concept involves a spacecraft equipped with an ion thruster to direct a controlled ion beam at the debris, generating a small force that gradually lowers its orbit. A proposed configuration of the chaser’s actuator system integrates radial and out-of-plane cold-gas thrusters along with in-track ion thrusters to enhance control and safety while maintaining low mission costs. A robust Model Predictive Control (MPC) strategy is implemented, using the theory of MPC for Tracking to ensure accurate positioning and effective deorbiting. This theoretical approach addresses uncertainties and perturbations to robustly guarantee safe distances between the chaser and the debris. Additionally, a new ray-marching-based algorithm is introduced to estimate the force and torque exerted by the ion beam on the target, considered as a 6 degrees of freedom object, improving simulation accuracy and control performance assessment. A comprehensive simulation of the deorbit of a large debris object is performed, demonstrating the potential of the IBS technology for future large-debris removal missions. This research advances the conceptual framework and control techniques for the IBS technology, advancing towards its future implementation in space debris mitigation. 1. Introduction The proliferation of space debris has become a critical issue in contemporary space exploration and satellite operations. According to [1], there are more than 30 000 tracked debris objects orbiting Earth, with millions of smaller fragments posing additional risks. These debris objects, traveling at high velocities, can cause significant damage to operational spacecraft, endangering space missions and the lives of astronauts aboard the International Space Station. The accumulation of debris also raises concerns about the potential onset of Kessler syndrome [2], a scenario where cascading collisions among space objects lead to a self-sustaining chain reaction of debris generation, further aggravating the problem. To mitigate this risk, the space industry has proposed various debris removal strategies [3], which can be broadly categorized as contact and non-contact methods. Contact methods often require precise maneuvering and docking capabilities, whereas contactless methods can operate at a distance, but may require high power or advanced technology to achieve the desired effects. Examples of contact methods are net captures, harpoons, or robotic arms, while the main proposed ∗Corresponding author. E-mail addresses: [email protected] (J. Urrios), [email protected] (R. Vazquez), [email protected] (F. Gavilan), [email protected] (I. Alvarado). technologies for contactless removal are laser ablation and ion beams. Ongoing efforts are focused on translating these strategies into practical applications. For example, Astroscale, a Japanese company, recently demonstrated its capabilities by capturing images of a rocket upper stage from a distance of several hundred meters using the ADRAS-J satellite. This mission aims to remove the debris using a robotic arm. Moreover, NASA’s Office of Technology, Policy, and Strategy recently published an analysis on the costs and benefits of contact methods for space debris remediation [4]. While contact methods are effective, they require direct contact with the debris, making them less suitable for faster spinning objects. The Ion Beam Shepherd (IBS) technology (which was nearly simultaneously proposed by [5–7]) represents a promising alternative for space debris mitigation. The IBS system uses a spacecraft equipped with an ion thruster to project a controlled stream of charged particles (ions) towards a target debris object. This ion beam generates a small force, allowing the debris to be gradually deorbited into Earth’s atmosphere for disintegration. Usually, a target with an almost circular orbit is chosen (although some studies consider higher eccentricities, such https://doi.org/10.1016/j.actaastro.2024.09.054 Received 1 August 2024; Received in revised form 17 September 2024; Accepted 22 September 2024 Acta Astronautica 225 (2024) 689–706 Available online 26 September 2024 0094-5765/© 2024 The Authors. Published by Elsevier Ltd on behalf of IAA. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ). J. Urrios et al. Fig. 1. Ion Beam Shepherd concept. as [8] or [9]), and the shepherd (or chaser) is set to travel the same orbit but a bit ahead. The chaser then points its ion beam backwards to slow down the debris (see Fig. 1). This technology offers several advantages: (i) it avoids direct contact with debris, increasing the safety of operations; (ii) it works independently of the target shape, rotational state or material; (iii) it is reusable until fuel depletion, potentially enabling multi-target Active Debris Removal (ADR) missions; (iv) taking into account that the momentum exchange is performed remotely, it is relatively efficient given its high specific impulse and that ions are concentrated in a narrow cone; and (v) it is based on mature and space-tested electric propulsion technology. Moreover, the IBS is well suited to tackle large space debris objects, which are of special interest. Indeed, according to [4], the deorbiting of large debris objects is highly beneficial, because it not only eliminates the necessity of spacecraft to constantly maneuver to avoid collisions, but also mitigates the generation of smaller debris resulting from surface degradation. Additionally, the deorbiting of these large objects effectively neutralizes their high potential for further debris creation. However, the ion beam also poses demanding challenges; the relative position between the chaser and the target must be accurately kept for the method to be effective, and formation flying with the target at such a small distance requires a robust control that guarantees no collision and deals with uncertainty and perturbations. For the IBS method to be effective, the ion beam should be accurately and efficiently characterized, in order to properly compute its interaction with the debris. Efforts towards these tasks were made in [10–12], and even some in-house software was developed to study such interaction [13]. The next task was to build a subsystem capable of generating such a beam and optimize it to improve the beam convergence [12] and minimize power consumption [14]. Guidelines and constraints to be taken into account in the preliminary design of an IBS mission were outlined in [15]. Some algorithms for on-board force computation (although not torque) have also been developed [16,17], the latter leveraging deep learning techniques. Furthermore, additional research has investigated more advanced consequences of the ion beam–target interaction such as surface sputtering and backscattering contamination, through numerical simulations [18] and experiments [12,19], with the aim of further improving the models. In terms of the relative positioning and control of the chaser with respect to the target, different configurations have been developed. Conventionally, the shepherd satellite consists of the main ion beam thruster, usually called the Impulse Transfer Thruster (ITT), which is the one that points to the space debris object (target), and another Impulse Compensation Thruster (ICT) in the opposite direction, to counteract the reaction force generated by the ITT [7]. In order to extract energy as efficiently as possible, the ion beam should act in the opposite direction of the target’s velocity, that is, it should be positioned in front of the debris. Initially, as part of a preliminary study, in [20] it was assumed that the chaser had control thrusters capable of continuous control in all directions, assumed no measurement errors and a spherical target (and thus no considerations in terms of attitude). Then, [21] extended the formulation to include state errors while preserving robustness, and made use of Pulse-Width Modulated hydrazine thrusters, while maintaining the spherical form hypothesis. Later, the attitude evolution of the target was thoroughly studied [22,23], along with possible strategies to control it [9,24]. The efficiency of the relative position control presented in [20] can be further enhanced by considering alternate control strategies. For instance, control in the in-track direction could be achieved by modulating the force exerted by the ICT. This approach was investigated in [25], where the problem was simplified to a planar case, and the system was controlled solely through ICT modulation, without any control in the radial direction. This control strategy demonstrated the potential for reduced fuel consumption. To achieve this, and 𝐻∞based controller was implemented for a time-varying and parametric uncertain plant. However, the chaser remained subject to periodic relative position errors due to the low controllability of the system and the small but non-zero eccentricity of the orbits. Later, more eccentric orbits were considered in [8], where a periodic Linear Quadratic Regulator (LQR) with a periodic reference input was implemented to minimize periodic errors due to the time-varying plant. These relative position errors were reduced to around 1 m for a LEO orbit with 0.1 eccentricity. Moreover, the work in [26] extended the formulation in [25] to the three-dimensional case, and discovered that it was also feasible to control the system in the out-of-plane direction by only performing yaw attitude deviations, which reduces fuel consumption since there is no need for an additional thruster. More recently, a new work proposes multiple ITTs to control the relative position of the chaser and the attitude of the target simultaneously, while deorbiting the target faster [27]. Although the developments mentioned above are interesting and crucial, some of them may not be suitable for implementation in the near future given the low Technology Readiness Level (TRL) of the IBS technology. For example, even if the yaw control in [26] works in simulations, it might not be safe to send such a satellite without a highly responsive controller in the out-of-plane direction. Similarly, although it is technically possible to control in-plane dynamics just with modulation of the ICT [25], adding an actuator in the radial direction is a reasonable way to reduce the probability of impact between the chaser and the target. Lastly, as powerful as the idea of using multiple beams for simultaneous attitude and relative position control [27] sounds, for an initial mission, it is advisable to aim for a minimal number of ion beams to mitigate the risk of component or system malfunction. By reducing the complexity of the propulsion system, the overall reliability of the spacecraft can be enhanced, thereby increasing the likelihood of mission success. This paper proposes a new configuration of the chaser actuator system, where the ITT and ICT in the in-track directions are maintained, but on–off cold-gas thrusters are added in the radial and out-of-plane directions (see Fig. 2). Cold-gas thrusters are simple and reliable actuators, especially the ones that operate in a binary on–off mode, and the addition of them in the radial and out-of-plane directions ensures controllability and security, while not increasing the mission budget too much. For the control, a Robust Model Predicitive Controller for Tracking is implemented. Model Predictive Control (MPC) [28] is an advanced control strategy used in a variety of industries, including aerospace. Unlike traditional control methods, MPC looks into the future by predicting the system’s behavior over a specified time horizon and optimizing control actions accordingly. It utilizes a dynamic model of the system, which is usually Linear Time-Invariant (LTI), to simulate future behavior and selects control inputs that minimize a defined cost function, accounting for constraints on states and inputs. MPC has already been used in similar use cases, such as rendezvous [29], also in the presence of on–off thrusters [30]. There are important aspects that must be taken into account in the design of an MPC controller, mainly feasibility, stability, and robustness [31]. Feasibility refers to the ability of the MPC controller to find Acta Astronautica 225 (2024) 689–706 690 J. Urrios et al. a solution that satisfies all the constraints imposed on the system while achieving the control objectives. Stability is the ability of the closedloop system to remain bounded and converge to a desired operating point or trajectory over time. Lastly, robustness indicates the ability of the controller to maintain satisfactory performance in the presence of uncertainties, disturbances, and model mismatches. There are several formulations in which these aspects can be tackled. In this work, a Model Predictive Controller based on the formulations developed in [32,33] is implemented, which combines MPC for Tracking and Robust MPC. Robust MPC has been merged with MPC for Tracking before (see [34]), employing the notion of tubes to ensure robustness. Here, however, to ensure robustness, the constraint tightening approach in [32] is taken instead, and thus the final controller is different. MPC for Tracking [33] was initially developed to enable systems to follow piece-wise changing references without losing the feasibility of the controller. This is not of much use in this work, since the chaser should regulate itself to a nominal nonchanging relative position with respect to the target. However, another consequence of implementing the MPC for Tracking scheme is that the region of attraction of the controller (region in which the controller successfully operates) is expanded, enhancing its capabilities and, in this case, further reducing fuel consumption. Robust MPC (for example, the formulation in [32]) was developed to account for uncertainties in the dynamic model given to the MPC. These uncertainties could be caused by perturbations, linearization errors, or others. Assuming the uncertainties are unknown but bounded, the goal of the Robust MPC scheme is to find a solution that minimizes the cost function while ensuring the constraints are met for any value of the uncertainties, and guaranteeing recursive feasibility and stability. This formulation is extremely useful for the IBS deorbiting scenario, since robust constraints on the position of the chaser can be set so that the ion beam always impacts the target efficiently for any target’s attitude state. The contribution of this work lies in the conceptual design of a robust control system for a realistic IBS mission scenario, utilizing a new actuator configuration. Following the literature [25,27], the deorbiting of a large debris object has been considered. However, it must be acknowledged that, once in very low orbits, the process of reentry from such large objects would be random with associated unacceptable risks on ground and thus a secondary mechanism to ensure a safe reentry would be required, which is beyond the scope of this paper. The mission aims to deorbit the current most statistically concerning debris object according to [35], the Zenit-2 second stage rocket body (NORAD ID: 28353), down to an orbit of 340 km, although in principle this framework can be used to deorbit many other debris objects. Perturbations and unmodeled dynamic effects are bounded and fed to the controller to ensure the chaser maintains a safe distance from the target, while preserving a fairly large region of attraction and not impacting fuel consumption. This proposed configuration incorporates cold-gas thrusters in the radial and out-of-plane directions, alongside the existing in-track ion thrusters (ITT and ICT). This design concept aims to enhance controllability and safety without significantly increasing the mission budget. Additionally, this paper introduces a novel method for computing the force and torque exerted by the ion beam on the target using a ray-marching algorithm that is integrated into the mission simulator. In this way, the performance of the controller can be measured more accurately. By addressing these aspects, this work advances the theoretical understanding and control algorithms of the IBS technology for large space debris mitigation, providing a potential solution for future implementation. This paper is organized as follows. In Section 2, the problem statement is presented and the equations of motion are introduced, together with the modeled perturbations. Section 3includes an exhaustive explanation of how the ion beam is modeled, introducing a new algorithm based on ray-marching to compute the force and torque exerted on the Fig. 2. Proposed IBS configuration. target. Next, Section 4tackles all aspects related to control, introducing a linear model of the system, the constraints to meet by the controller, and the controller itself. Disturbance bounding is also performed in this section, and considerations regarding recursive feasibility, convergence, and the tuning of the controller are present. In Section 5, a simulation of the whole deorbiting mission is performed, and relevant results are shown and discussed. Lastly, the conclusions are summarized in Section 6, along with possible future work. 2. Problem statement and equations of motion The deorbiting of a large space debris object using IBS technology is considered in this paper. Depending on the dimensions of the target body, and leaving security margins to account for differences between the target’s Center of mass (COM) and geometric center, the distance between the chaser’s COM and the ion beam source, and some extra space in case of malfunctioning, the chaser will be positioned a certain distance 𝑑𝐶𝑇 in front of the target.1The chaser has to maintain this distant position accurately, because due to the separation between the two bodies, small deviations in the chaser position means the beam hits the target only partially, notably reducing efficiency. To perform the deorbit, the chaser must maintain a precise relative position with the target to correctly transfer the momentum of the ion beam. Thus, the spacecraft must have the proper actuators, and a suitable controller must be designed. The spacecraft is equipped with an ITT that maintains its thrust constant with a value of 𝐹𝐼𝑇 𝑇 and an ICT in the opposite direction that can vary its thrust 𝐹𝐼𝐶𝑇 for control purposes. In the perpendicular directions, the chaser is equipped with cold-gas thrusters. It is assumed that the chaser possesses an attitude control system capable of always orienting the ITT in the negative in-track direction, which in perfect circular orbits coincides exactly with the target direction. The chaser is assumed to have sensors capable of estimating the relative state of the target and also some orbit parameters used by the controller such as the mean motion 𝑛of the orbit. These parameters are estimated with a certain error that is bounded. In the subsequent sections, the following reference frames are used: •Earth-Centered-Inertial (ECI) 𝑂𝐼𝑥𝐼𝑦𝐼𝑧𝐼, which is an inertial reference frame centered on the Earth, coincident with J2000. •Local-Vertical Local-Horizontal frame (LVLH) 𝑂𝐿𝑥𝐿𝑦𝐿𝑧𝐿. Although this frame is not used by itself, it helps us introduce later frames. Generically, an LVLH frame is a frame centered on an object that orbits a certain body, and their axes are defined as follows. The 𝑥𝐿axis points towards the object position vector as measured from the orbited body, 𝑧𝐿points towards the angular momentum vector of the object, and 𝑦𝐿completes the Cartesian right-handed coordinate system. •Chaser-centered LVLH frame (CLVLH) 𝑂𝐶𝑥𝐶𝑦𝐶𝑧𝐶. It is the LVLH frame centered on the chaser’s COM and orbiting the Earth (see Fig. 3). 1See Nomenclature at end of paper Acta Astronautica 225 (2024) 689–706 691 J. Urrios et al. Fig. 3. Chaser-centered LVLH frame definition. Fig. 4. Block diagram of the problem. •Target-centered LVLH frame (TLVLH) 𝑂𝑇𝑥𝑇𝑦𝑇𝑧𝑇. Analogous to the CLVLH, but the origin is now the target’s COM, which is assumed to coincide with its geometric center. •Target Body Axes (TBA) 𝑂𝐵𝑥𝐵𝑦𝐵𝑧𝐵. With its origin at the target’s COM, it is fixed to the target, with its 𝑧𝐵axis coinciding with the cylinder revolution axis. The state variables of the problem are the position and velocity of both the chaser and the target in the ECI frame, 𝐫𝐼 𝐶,𝐯𝐼 𝐶,𝐫𝐼 𝑇,𝐯𝐼 𝑇; the attitude quaternion of the target with respect to ECI 𝐪𝐵 𝐼; the angular velocity of the target in the body axes 𝝎𝐵 𝑇; and the mass of the chaser 𝑚𝐶. A schematic view of the system is included in Fig. 4. 2.1. Equations of translational motion The target and chaser dynamics are given by the restricted two-body problem with added perturbations. More specifically, the equations of motion of the COMs of the chaser and target in the ECI frame are  𝐫𝐼 𝐶= − 𝜇𝐸 𝑟3 𝐶 𝐫𝐼 𝐶+𝐚𝐼 𝐶,𝐵 +𝐚𝐼 𝐶,𝐶 +𝐚𝐼 𝐶,𝐷 +𝐚𝐼 𝐶,𝐽2(1)  𝐫𝐼 𝑇= − 𝜇𝐸 𝑟3 𝑇 𝐫𝐼 𝑇+𝐚𝐼 𝑇 ,𝐵 +𝐚𝐼 𝑇 ,𝐷 +𝐚𝐼 𝑇 ,𝐽2,(2) where 𝜇𝐸= 3.9860044 ⋅105km3/s2is the geocentric gravitational constant, and the subscripts 𝐶,𝑇represent the chaser and target respectively. The first term of both equations represents the unperturbed dynamics of the two-body problem, while the rest are accelerations due to the ion beam (𝐚𝐼 𝐶,𝐵 and 𝐚𝐼 𝑇 ,𝐵,𝐵subscript from beam), control (𝐚𝐼 𝐶,𝐶 , 𝐶subscript from control), drag (𝐚𝐼 𝐶,𝐷 and 𝐚𝐼 𝑇 ,𝐷,𝐷subscript from drag) and J2 perturbation (𝐚𝐼 𝐶,𝐽2and 𝐚𝐼 𝑇 ,𝐽2). Drag and J2 perturbations are chosen because they are the most prominent for a satellite in a Low Earth Orbit (LEO). 2.1.1. Ion beam acceleration The force 𝐅𝐶 𝐵exerted by the ion beam on the target is calculated using the method explained in Section 3, and it is expressed in the CLVLH frame. Therefore, the acceleration is obtained as 𝐚𝐼 𝑇 ,𝐵 =𝑅𝐼 𝐶 𝐅𝐶 𝐵 𝑚𝑇 ,(3) where 𝑅𝐼 𝐶is the DCM that transforms from CLVLH to ECI. The equilibrium state of the system is chosen so that the ion beam fully impacts the target, maximizing the force transmitted to it (𝐹𝐶,𝑚𝑎𝑥 𝐵=𝐹𝐼𝑇 𝑇 ) and thus imparting a maximum acceleration of 𝐚𝐼,𝑚𝑎𝑥 𝐵,𝑇 . For the system to remain in this equilibrium (ignoring perturbations and unmodeled dynamics), the nominal value of the ICT must be such that the acceleration of the chaser and the target due to the beam are equal 𝐚𝐼 𝐶,𝐵 =𝐚𝐼,𝑚𝑎𝑥 𝑇 ,𝐵 .(4) Note that during the mission, the ion beam may impact the target only partially. The assumption in (4) is only considered to linearize the system along that configuration. 2.1.2. Control acceleration There is a component of the acceleration that the chaser experiments that comes from the action of controlling the system, 𝐚𝐼 𝐶,𝐶 . The design of the controller as well as the actuators’ behavior are explained in Section 4. Logically, there is no such acceleration term for the target. 2.1.3. Drag perturbation To model drag, the next simple model is used 𝐚𝑖,𝐷 = − 1 2𝜌𝑎(ℎ𝑖)𝑣𝑟𝑖𝐯𝑟𝑖 𝐴𝑓,𝑖 𝑚𝑖 𝑐𝐷,𝑖,(5) where the subindex 𝑖stands for chaser 𝐶and target 𝑇respectively; 𝜌𝑎(ℎ𝑖)is the atmosphere density, which is computed using the International Standard Atmosphere (ISA) model knowing the height ℎ𝑖;𝐯𝑟𝑖 is the relative velocity between the object and the air (the latter is assumed to rotate in solidarity with the Earth); 𝐴𝑓,𝑖 is the frontal area; and 𝑐𝐷,𝑖 is the drag coefficient. The chaser frontal area is assumed to be constant since it always points in the same direction, while the target frontal area depends on its attitude. Knowing the attitude with respect to the TLVLH frame 𝑅𝐵 𝑇, the relative velocity in that same frame 𝐯𝑟𝑇 𝑇, and that the cylinder symmetry axis goes along 𝑧𝐵, the frontal area can be computed as 𝛿= min (arccos ([001]𝑅𝐵 𝑇 𝐯𝑟𝑇 𝑇 𝑣𝑟𝑇), arccos ([0 0 −1]𝑅𝐵 𝑇 𝐯𝑟𝑇 𝑇 𝑣𝑟𝑇)),(6) 𝐴𝑓,𝑇 = 2𝑅𝑐𝐻𝑐sin 𝛿+𝜋𝑅2 𝑐cos 𝛿, (7) where 𝑅𝑐is the cylinder radius, 𝐻𝑐its height, and 𝛿is the angle between the cylinder axis and the in-track direction (which is between 0 and 90 degrees). 2.1.4. J2 perturbation J2 perturbation, with its origin in Earth’s oblateness, is also included. This kind of harmonic perturbations, which appear due to the fact that the Earth is not a perfect sphere, are computed in the Earth-Centered Earth-Fixed frame (ECEF) for obvious reasons. However, taking into account that the 𝑧axes of both ECI and ECEF are almost coincident, and that J2 only depends on the latitude of the satellite (and not longitude), the J2 perturbation can be expressed directly in the ECI frame almost without losing accuracy. The acceleration due to J2 is 𝐚𝐼 𝑖,𝐽2=3 2𝜇𝐸𝐽2 𝐫𝐼 𝑖 𝑟3 𝑖(𝑅𝐸 𝑟𝑖)2 ◦⎛⎜⎜⎝ 5(𝑧𝐼 𝑖 𝑟𝑖)2⎡⎢⎢⎣ 1 1 1⎤⎥⎥⎦ −⎡⎢⎢⎣ 1 1 3⎤⎥⎥⎦⎞⎟⎟⎠ ,(8) where again the subscript 𝑖stands for chaser 𝐶and target 𝑇respectively, 𝐽2= 0.0010826269 and 𝑅𝐸= 6378.1366 km is the radius of the Earth. The operator ◦denotes the Hadamard product or element-wise product of two vectors. Acta Astronautica 225 (2024) 689–706 692 J. Urrios et al. 2.2. Equations of rotational motion When the beam impacts the space debris object, it exerts not only a force but also a torque, which is also computed by the algorithm presented in Section 3. This torque changes the angular velocity of the debris, and it could potentially accelerate it indefinitely, which is not desirable since phenomena such as centrifugal fragmentation could occur. To ensure this does not happen, and to simulate a more realistic scenario (the frontal area of the debris changes depending on its attitude), the rotational dynamics are also included in the simulation:  𝐪𝐵 𝐼=1 2𝐪𝐵 𝐼∗[0 𝝎𝐵 𝑇],(9)  𝝎𝐵 𝑇=−1 (𝐌𝐵 𝐵−𝝎𝐵 𝑇×𝝎𝐵 𝑇),(10) where 𝐪𝐵 𝐼is the scalar-first attitude quaternion that represents the attitude of the target with respect to ECI, ∗represents quaternion multiplication, 𝝎𝐵 𝑇is the target angular velocity, 𝐌𝐵 𝐵is the torque exerted by the beam, and is the inertia tensor of the target, with the latter three quantities expressed in the body frame. In this work, the target is modeled as a cylinder with constant density, which means that the COM and the geometric center coincide. With these assumptions, the inertia tensor is computed as 𝑥=𝑦=1 12 𝑚𝑇𝐻2 𝑐+1 4𝑚𝑇𝑅2 𝑐,(11) 𝑧=1 2𝑚𝑇𝑅2 𝑐,(12) 𝑥𝑦 =𝑦𝑧 =𝑥𝑧 = 0,(13) where 𝑚𝑇is the target mass, and 𝑅𝑐and 𝐻𝑐are the cylinder radius and height respectively. Note that the only torque source for the rotational dynamics is assumed to be the ion beam. There could be more, such as gravity gradient or drag, but they are small and slow, and have been neglected for simplicity since they would have little impact on rotational motion. 2.3. Fuel consumption As the chaser is controlled with both cold-gas thrusters and ion beams, it uses fuel and thus loses mass. This mass variation is important, as it changes the force that the actuators need to exert to produce the same acceleration. Note that the rules that are commanded to the actuators are computed at the sampling time, even if they are applied some time later. The forces actuating on the chaser are the ITT and the ICT in the in-track direction, and the cold-gas thrusters in the radial and out-of-plane directions. Recall from Section 2.1.1 that ignoring perturbations and unmodeled dynamics, and assuming that the beam completely impacts the target, the ICT force must be such that the system remains in equilibrium with no control (𝐚𝐼 𝐶,𝐶 = 0, 𝐚𝐼 𝐶,𝐶 is now referred to as  𝐮for clarity). In that scenario, at a generic time instant, the in-track accelerations of the target and the chaser are 𝑎𝑇 ,𝑖𝑡 = −𝐹𝐼𝑇 𝑇 ∕𝑚𝑇and 𝑎𝐶,𝑖𝑡 = (𝐹𝐼𝑇 𝑇 −𝐹𝐼𝐶𝑇 )∕𝑚− 𝐶respectively, where 𝑚− 𝐶 represents the chaser mass at the most recent sampling time. Therefore, for a generic acceleration 𝑢𝑦commanded in the in-track direction, the ICT must be such that 𝑢𝑦=𝑎𝑖𝑡,𝐶 −𝑎𝑖𝑡,𝑇 , that is 𝐹𝐼𝐶𝑇 =𝐹𝐼𝑇 𝑇 (1 + 𝑚− 𝐶 𝑚𝑇)−𝑚− 𝐶𝑢𝑦.(14) Thus, the complete fuel consumption equation is given by 𝑚𝐶= − (2 + 𝑚− 𝐶 𝑚𝑇)𝐹𝐼𝑇 𝑇 −𝑚− 𝐶𝑢𝑦 𝑔0𝐼𝑠𝑝,𝐼𝐵 −𝑚− 𝐶(|𝑢𝑥|+|𝑢𝑧|) 𝑔0𝐼𝑠𝑝,𝐶𝐺 ,(15) where 𝐼𝑠𝑝,𝐼𝐵 and 𝐼𝑠𝑝,𝐶𝐺 are the specific impulses of the ion beam and the cold-gas thruster respectively, while 𝑔0= 9.80665 m∕s2is the standard gravity. Fig. 5. Depiction of the impact of the ion beam in the target, with a generic relative position and attitude of the later. 2.4. System uncertainties Within the system, uncertainties related to sensoring and actuation are included. These uncertainties are modeled as uniform distributions so that they have a bounded maximum value. This is necessary to ensure the robustness of the controller in Section 4. If uncertainties were to be modeled as normal distributions or similar, a chance-based approach can be used, in which it is possible to give a probability (which is now not 100%) that the constraints will be met (see, for example, [29]). The spacecraft has sensors that estimate the relative position and velocity of the target with respect to the chaser, and these measurements are assumed to have some error. If some advanced sensoring system based on lidar/radar is considered, bounding values of the relative position errors are of the order of ±1 cm for the relative position and ±5 cm/min for the relative velocity. Moreover, the ion beam force is assumed to be not completely stationary, slightly changing the divergence angle of the beam. This effect is modeled by considering that the force exerted on the target can vary as much as ±5% of the nominal value. 3. Ion beam modeling. Force and torque computation Modeling the ion beam is crucial for an accurate simulation of the controlled system. However, computing the force and torque that the beam exerts on the target is not trivial. With an arbitrary position and attitude of the target, perhaps only part of the beam impacts the object and the rest does not (see Fig. 5), and the points where it impacts are crucial for force and torque calculations. To solve this task, some groups discretize the target surface as a triangle mesh and locally compute the force and torque on each one [22]. Others use properties of the assumed ion beam model to project the force and compute it in two dimensions [16]. Another group has even developed internal software for this task [13], based on finite element methods and a shadowing algorithm. In this work, another approach based on ray-marching is chosen, and it is explained in this section. 3.1. Ray-marching basics Ray-marching is a rendering technique used in computer graphics to generate images by casting rays into a scene and incrementally advancing along each ray to find intersections with objects [36]. In this paper, no images are generated using this technique, but rather the technique is used to compute the distance from a point to the first intersecting object in some given direction, that is, from the ion beam cone vertex to the target in some given direction. The idea is that one can find the distance by iteratively advancing along the ray a distance equal to the shortest distance to the closest object, as Fig. 6 shows. Acta Astronautica 225 (2024) 689–706 693 J. Urrios et al. Fig. 6. Ray-marching depiction. Fig. 7. Relationship between a target differential surface element d𝑆and a beam differential surface element d𝑆′. If that shortest-distance computation is cheap (which is the case for simple geometrical shapes, such as a cylinder, as shown in Section 3.3), then the algorithm converges rapidly to the desired distance. Raymarching may not be suited for more complex geometries, but to the authors’ knowledge, no complex geometry for the space debris object has been used yet for the study of a potentially viable IBS deorbit mission, and even then, more complex shapes generated from Boolean operations applied on simple mathematical shapes are also suited for ray-marching. 3.2. Beam-focused discretization of the force Under some typical simplifying assumptions [20], the differential force that the ion beam exerts in a differential surface element of the target is given by d𝐅𝐵=𝑚𝐼𝑛𝐼𝐔(−𝐍⋅𝐔)d𝑆, (16) where 𝑚𝐼is the ion mass, 𝑛𝐼is the number density of ions, 𝐔is the ion velocity vector, and 𝐍is the external unit normal vector to the elementary area d𝑆. Under the same assumptions, a self-similar model of the fluid can be developed to determine the ion number density and velocity at any point far from the beam source, which is the region of interest. This model assumes the beam expands in a conical shape, disregarding factors like internal ambipolar electric fields [10]. However, this simplification is common in similar studies [16,20,25,27], as it has limited impact on the overall behavior of the beam. For the ray-marching-based algorithm implementation, however, it is more interesting to discretize the force in terms of a differential surface element of the beam, not the target. To do this, observe Fig. 7, where a generic differential surface element of the target d𝑆is shown along with a differential surface element of the beam d𝑆′. Note that, since the ions move radially away from the source, 𝐔and d𝑆′are always perpendicular. Furthermore, since d𝑆and d𝑆′are arbitrarily small, they can be considered flat surfaces, and 𝐔is assumed constant throughout their span. Note that Fig. 7 is not a 2D simplification of the 3D case, it represents the actual 3D flat differential surfaces seen from a convenient point of view: perpendicular to the plane that contains the normal vectors of both d𝑆and d𝑆′. Looking at Fig. 7, it is not difficult to see that since −𝐍⋅𝐔=𝑈cos 𝛾, and cos 𝛾d𝑆=d𝑆′, the differential force can be expressed in terms of a differential surface element of the beam as d𝐅𝐵=𝑚𝐼𝑛𝐼𝐔𝑈d𝑆′.(17) Now, spherical coordinates as in Fig. 8 are defined, where 𝑂𝐶𝑥𝐶𝑦𝐶𝑧𝐶represents the CLVLH frame and 𝑑𝑉is the distance from the Fig. 8. Definition of spherical coordinates for the ion beam computation. Table 1 Ion beam characteristics. 𝑅0[m] 𝑚𝐼[kg] 𝑛0[m−3]𝑈0[m/s] 𝛼0[deg] 0.1894 2.18 ⋅10−25 4.13 ⋅1015 71,580 7 COM of the chaser to the cone vertex of the ion beam. The value of 𝛼 ranges between 0 and 𝛼0, which is the angle of divergence of the beam, and 𝜌∈ [0,∞),𝛽∈ [0,2𝜋). In these coordinates, a differential surface element of the beam, the ion number density, and the ion velocity are given by d𝑆′=𝜌2sin 𝛼d𝛼d𝛽, (18) 𝑛𝐼= 𝑅2 0𝑛0 𝜌2tan2𝛼0cos2𝛼exp (−𝐶tan2𝛼 2 tan2𝛼0),(19) 𝐔=𝑈0 cos 𝛼 𝝆,(20) where  𝝆is the radial unit vector, 𝑅0is the radius at a reference plane of the plasma tube that contains 95% of the ion beam mass flow, 𝑛0is the density at the center of that reference plane, 𝐶≈ 6 is an integration constant related to that 95%, and 𝑈0is the ion velocity magnitude at the cone axis (see [16]). Therefore, the differential force can be computed as d𝐅𝐵= 𝑚𝐼𝑛0𝑅2 0𝑈2 0tan 𝛼 tan2𝛼0cos3𝛼exp (−𝐶tan2𝛼 2 tan2𝛼0) 𝝆d𝛼d𝛽. (21) Note that this expression does not depend on the distance between the ion beam impact point and the source point. However, it is still necessary to compute the distance for two reasons. The first is that, even if the distance itself does not appear, it must be known if the beam in that specific direction (given by 𝛼and 𝛽) intersects the target, because if it does not, then d𝐅𝐵is zero. The second reason is that even if the distance is not needed for the force computation, it is needed for the torque, and that is probably why algorithms in [16,17] do not compute it. For the actual computation, taking into account that the chosen ion thruster is Xenon-based, the same value as in [16,17] is used for the ion mass. The ion density at the center of the reference plane, the ion velocity magnitude at the cone axis, and the divergence angle are also taken from those studies. What is modified, then, is 𝑅0, which is adapted so that the force that the ion beam exerts in an object when it impacts entirely is 𝐹𝐼𝑇 𝑇 = 164.8mN (the value in Table 6). The characteristics of the beam are included in Table 1. 3.3. Ray-marching to compute the distance For a cylinder of radius 𝑅𝑐and height 𝐻𝑐centered at the origin and with its symmetry axis along 𝑧, the distance 𝜌(𝛼, 𝛽)from an arbitrary point 𝐐𝐵in the target body axes to the cylinder can be computed with just one formula as 𝜌=‖‖‖‖‖‖‖⎛⎜⎜⎝ max (√𝑄𝐵 𝑥 2+𝑄𝐵 𝑦 2−𝑅𝑐,0) max (|𝑄𝐵 𝑧|−𝐻𝑐∕2,0)⎞⎟⎟⎠‖‖‖‖‖‖‖ ,(22) Acta Astronautica 225 (2024) 689–706 694 J. Urrios et al. and, with the cylinder center in a certain position with respect to the CLVLH frame 𝐫𝐶 𝑂𝐶𝑦𝑙 and an arbitrary attitude 𝑅𝐵 𝐶,𝐐𝐵can be easily obtained as 𝐐𝐵= −𝑅𝐵 𝐶(𝐫𝐶 𝑂𝐶𝑦𝑙 −𝐝𝐶 𝑉),(23) where 𝐝𝐶 𝑉is the position vector of the ion beam cone vertex in the CLVLH frame. If the ion beam is always pointed towards −𝑦𝐶, then 𝐝𝐶 𝑉= [0,−𝑑𝑉,0]𝑇. 3.4. Discretizing the beam To discretize the beam, simple linear spacings in both 𝛼and 𝛽are taken. In Section 3.5 it is shown that 𝑁𝛼= 100 and 𝑁𝛽= 20 points for 𝛼and 𝛽respectively are enough to approximate the beam force to an accuracy of around 0.15% of relative error, while still being fairly fast. For each 𝛼𝑖,𝛽𝑗, one computes  𝝆𝐶 𝑖𝑗 =⎛⎜⎜⎝ cos 𝛽𝑗sin 𝛼𝑖 − cos 𝛼𝑖 sin 𝛽𝑗sin 𝛼𝑖⎞⎟⎟⎠ ,(24) 𝛥𝐅𝐶 𝐵,𝑖𝑗 = 𝑚𝐼𝑛0𝑅2 0𝑈2 0tan 𝛼𝑖 tan 𝛼2 0cos3𝛼𝑖 exp (−𝐶tan2𝛼𝑖 2 tan2𝛼0) 𝝆𝐶 𝑖𝑗 𝛥𝛼𝛥𝛽𝑏𝑖𝑗 ,(25) where 𝑏𝑖𝑗 is a boolean variable that is 0 if the beam does not intersect the cylinder and 1 if it does. The force in the CLVLH frame, and torque in the TBA frame, which are the outputs of the algorithm, can be retrieved as 𝐅𝐶 𝐵= 𝑁𝛼 ∑ 𝑖=1 𝑁𝛽 ∑ 𝑗=1 𝛥𝐅𝐶 𝐵,𝑖𝑗 ,(26) 𝐌𝐵 𝐵=𝑅𝐵 𝐶 𝑁𝛼 ∑ 𝑖=1 𝑁𝛽 ∑ 𝑗=1 𝐫𝐶 𝑂𝐶𝑦𝑙,𝑃 𝑖𝑗 ×𝛥𝐅𝐶 𝐵,𝑖𝑗 ,(27) where 𝐫𝐶 𝑂𝐶𝑦𝑙,𝑃 𝑖𝑗 is the vector that goes from the cylinder center to the impact point of the beam, and can be computed as 𝐫𝐶 𝑂𝐶𝑦𝑙,𝑃 𝑖𝑗 =𝜌𝑖𝑗  𝝆𝐶 𝑖𝑗 −(𝐫𝐶 𝑂𝐶𝑦𝑙 −𝐝𝐶 𝑉).(28) 3.5. Sensitivity analysis and results validation After applying the discretization to the model, it is left to choose the appropriate values of 𝑁𝛼and 𝑁𝛽. Logically, for 𝑁𝛼→∞, 𝑁𝛽→∞, the discretized solution converges to the real one. However, as the number of points increases, so does the computation time, and, given that this function is called multiple times by the integrator, it is mandatory that it takes very few to execute. Because of this, even if the programming environment used for the development of this paper is MATLAB, a C++ MEX function has been programmed precisely for this computation of force and torque. To choose 𝑁𝛼and 𝑁𝛽, a large value is first set for both 𝑁𝛼= 𝑁𝛽=10,000 and the force and torque are calculated with them, results that can be assumed to be very close to the real solution and thus are denoted 𝐅𝐵,ℎ𝑝 and 𝐌𝐵,ℎ𝑝, with the subscript ℎ𝑝 denoting high precision. Then a series of values lower than 10,000 are tested for both 𝑁𝛼and 𝑁𝛽, and the force 𝐅𝐵and torque 𝐌𝐵are computed for each of them. Lastly, the relative error between each obtained result and the one corresponding to the high precision computation is calculated and shown in Table 2 for the force and Table 3 for the torque, to calculate the quality of the results. For these computations the relative position of the target with respect to the chaser in a chaser-centered LVLH frame is 𝐫𝐶 𝑂𝐶𝑦𝑙 = [1,−𝑑𝐶𝑇 ,0]𝑇(with 𝑑𝐶𝑇 = 12 m), and its attitude, described in Euler angles, is 𝜃1= −45◦, 𝜃2= −30◦, 𝜃3= 0 for the YXZ sequence. They are chosen so that the beam hits the target in a rather arbitrary attitude and just partially so that torque is generated (otherwise, torque relative errors would just be computer numeric noise), but the same trends are found in other relative positions and attitudes. Table 2 Values of 100 ‖𝐅𝐵−𝐅𝐵,ℎ𝑝 ‖ ‖𝐅𝐵,ℎ𝑝 ‖[%], where 𝐅𝐵is computed for each 𝑁𝛼, 𝑁𝛽, and 𝐅𝐵,ℎ𝑝 is the value for 𝑁𝛼=𝑁𝛽= 10000. 𝑁𝛼\𝑁𝛽10 20 50 100 1000 10,000 10 0.8638 0.8637 0.7932 0.7931 0.7614 0.7618 20 0.5330 0.5330 0.5426 0.5426 0.5319 0.5336 50 0.3418 0.2981 0.2621 0.2621 0.2648 0.2643 100 0.1955 0.1513 0.1394 0.1394 0.1358 0.1357 1000 0.0635 0.0213 0.0133 0.0133 0.0122 0.0124 10,000 0.0521 0.0097 0.0014 0.0010 0.0000 0 Table 3 Values of 100 ‖𝐌𝐵−𝐌𝐵,ℎ𝑝 ‖ ‖𝐌𝐵,ℎ𝑝 ‖[%], where 𝐌𝐵is computed for each 𝑁𝛼, 𝑁𝛽, and 𝐌𝐵,ℎ𝑝 is the value for 𝑁𝛼=𝑁𝛽= 10000. 𝑁𝛼\𝑁𝛽10 20 50 100 1000 10,000 10 0.6860 0.7484 0.6410 0.6583 0.6112 0.6132 20 0.3350 0.4071 0.4377 0.4526 0.4368 0.4399 50 0.3450 0.2878 0.2309 0.2323 0.2377 0.2369 100 0.2049 0.1386 0.1244 0.1270 0.1212 0.1212 1000 0.1133 0.0234 0.0105 0.0123 0.0105 0.0109 10,000 0.1112 0.0205 0.0031 0.0022 0.0000 0 Table 4 Computation time [ms] computed for each 𝑁𝛼, 𝑁𝛽. 𝑁𝛼\𝑁𝛽10 20 50 100 1000 10,000 10 0.049 0.058 0.124 0.233 2.005 18.43 20 0.054 0.079 0.190 0.373 3.691 34.53 50 0.110 0.191 0.460 0.882 8.871 85.78 100 0.202 0.375 0.904 1.776 17.01 168.7 1000 1.806 3.606 8.932 17.86 168.7 1659 10,000 17.61 36.04 87.61 173.2 1699 16,653 It can be seen that the model is relatively precise even for a small number of points, achieving a less than 1% error with only 𝑁𝛼= 𝑁𝛽= 10. The computation times for each combination are included in Table 4. The calculations were performed using a Windows 11 laptop equipped with an Intel Core i7-12700H processor and 32 GB of RAM. Noting that the accuracy is much more sensitive to 𝑁𝛼than to 𝑁𝛽, and taking a look at the computational times, it seems reasonable to choose 𝑁𝛼= 100, 𝑁𝛽= 20 as the values with which to perform the simulations. Now, to validate the results of the model with 𝑁𝛼= 100, 𝑁𝛽= 20, one can check the force it gives compared to the results obtained in [16] using both direct integration and central projection. For this task the cylinder is modified to have the same measurements, which are a radius of 1.1 m and a height of 2.6 m; and the beam characteristics are changed to coincide with those in [16]. The forces obtained for the three methods and the 15 cases they present (which are depicted in Fig. 9) are included in Table 5, concluding that the mean relative error is around 0.15%, as expected. However, this beam marching algorithm is also able to compute torque while still being efficient. Note that in Table 5, cases 1, 4, 7, 10 and 13 give the same exact force result. This is no coincidence. With the hypotheses made in Section 3.2, in which the ion beam expands as a cone, and the force magnitude does not depend on the distance, the force can be computed by projecting the debris surface impacted by the ion beam on a certain perpendicular plane (that is the base of how they compute the force in [16]). Therefore, if the ion beam completely impacts the target, the projected shape is a circle regardless of the attitude, and the force is the same. Thus, as long as the ion beam completely impacts the target, the target’s attitude does not change its acceleration. If the impact is partial, attitude will influence the result. 4. Control design To control the relative position of the chaser with respect to the target, a Robust MPC for Tracking controller is designed, based on Acta Astronautica 225 (2024) 689–706 695 J. Urrios et al. Table 5 Comparison of force results obtained for both methods 1 and 2 in [16] vs using beam marching. Case Direct integration Central projection Beam marching 𝐹𝐵,𝑥 [N] 𝐹𝐵,𝑦 [N] 𝐹𝐵,𝑧 [N] 𝐹𝐵,𝑥 [N] 𝐹𝐵,𝑦 [N] 𝐹𝐵,𝑧 [N] 𝐹𝐵,𝑥 [N] 𝐹𝐵,𝑦 [N] 𝐹𝐵,𝑧 [N] 1 0.000 0.000 2.986e−2 0.000 0.000 2.975e−2 0.000 0.000 2.979e−2 2 0.000 3.431e−5 2.943e−2 0.000 3.459e−5 2.943e−2 0.000 3.451e−5 2.947e−2 3 0.000 5.332e−4 1.764e−2 0.000 5.327e−4 1.766e−2 0.000 5.336e−4 1.768e−2 4 0.000 0.000 2.974e−2 0.000 0.000 2.975e−2 0.000 0.000 2.979e−2 5−7.700e−6 8.587e−5 2.888e−2−7.487e−6 8.636e−5 2.887e−2−7.681e−6 8.671e−5 2.891e−2 6−7.496e−6 5.313e−4 1.834e−2−7.120e−6 5.318e−4 1.833e−2−3.458e−6 5.359e−4 1.825e−2 7 0.000 0.000 2.975e−2 0.000 0.000 2.975e−2 0.000 0.000 2.979e−2 8 6.939e−6 4.858e−6 2.967e−2 6.720e−6 5.464e−6 2.968e−2 7.354e−6 6.098e−6 2.971e−2 9 1.314e−4 1.352e−4 2.747e−2 1.303e−4 1.374e−4 2.745e−2 1.297e−4 1.368e−4 2.752e−2 10 0.000 0.000 2.975e−2 0.000 0.000 2.975e−2 0.000 0.000 2.979e−2 11 6.939e−6 4.858e−6 2.967e−2 6.720e−6 5.464e−6 2.968e−2 7.354e−6 6.098e−6 2.971e−2 12 1.314e−4 1.352e−4 2.747e−2 1.303e−4 1.374e−4 2.745e−2 1.297e−4 1.368e−4 2.752e−2 13 0.000 0.000 2.975e−2 0.000 0.000 2.975e−2 0.000 0.000 2.979e−2 14 1.311e−5 1.321e−5 2.959e−2 1.297e−5 1.297e−5 2.959e−2 1.243e−5 1.243e−5 2.964e−2 15 1.575e−4 2.259e−4 2.578e−2 1.566e−4 2.261e−4 2.578e−2 1.529e−4 2.281e−4 2.583e−2 Fig. 9. The 15 different configurations in which the force results are compared to those in [16]. the formulation in [34]. As typical in MPC, the system model for the controller is a linear approximation of the real system. Errors due to these approximations and from other sources are bounded and fed to the controller to make it robust. In addition, one must take into account that the controller output is continuous and is applied to the system via a zero-order-hold (ZOH). However, the cold-gas thrusters operate in a binary on–off mode, regulating the time they are actuating. To convert from one type of actuation to the other, an on–off filter is also included. This section is organized as follows: the on–off filter is explained in Section 4.1, the linear model approximation of the system is presented in Section 4.2, the system constraints are introduced in Section 4.3, the set that bounds the disturbances is computed in Section 4.4, and finally, Section 4.5 includes the controller design. 4.1. On–off filter for the cold-gas thrusters The control output as computed by the MPC lies in the reals and is assumed to be applied as a ZOH. Indeed, this is possible in the in-track direction, where the force of the ITT can be modulated as needed. However, one must remember that in the radial and out-of- plane directions the actuators operate in a binary on–off mode, which means they apply only either zero or maximum thrust, and are able to correctly actuate on the system by conveniently adjusting the times they are switched on and off. To convert the MPC control output into a feasible actuation, a centered equal-area filter is applied, so that the real control commanded to the actuators is 𝑢𝑖(𝑡) = ⎧ ⎪ ⎨ ⎪ ⎩ 0 0 ≤𝑡<𝑇𝑜𝑛 𝑖 𝑢𝑚𝑎𝑥,𝐶𝐺 ⋅sgn(𝑢𝑖)𝑇𝑜𝑛 𝑖≤𝑡≤𝑇𝑜𝑓𝑓 𝑖 0𝑇𝑜𝑓𝑓 𝑖< 𝑡 ≤𝑇𝑐, (29) where 𝑇𝑜𝑛 𝑖=𝑇𝑐 2(1 − |𝑢𝑖| 𝑢𝑚𝑎𝑥,𝐶𝐺 ),𝑇𝑜𝑓𝑓 𝑖=𝑇𝑐 2(1 + |𝑢𝑖| 𝑢𝑚𝑎𝑥,𝐶𝐺 );𝑖= {𝑥, 𝑧}, that is, this filter is applied only in the radial and out-of-plane directions; and 𝑢𝑚𝑎𝑥,𝐶𝐺 =𝐹𝐶𝐺∕𝑚− 𝐶is the maximum acceleration the cold-gas thrusters can produce at the control sampling time. The control is assumed to be recomputed every 𝑇𝑐seconds. Note that this filter centers the actuation in the middle of the control interval and adjusts the control so that ‘‘control areas’’ are preserved (that is, ∫𝑇𝑐 0𝑢𝑖d𝑡=∫𝑇𝑐 0𝑢𝑖), which forces both controls to have a ‘‘similar’’ impact on the system. There is one last feature of this filter. If the control action ||𝑢𝑖|| commanded by the MPC is so small that 𝑇𝑜𝑓𝑓 𝑖−𝑇𝑜𝑛 𝑖< 𝛥𝑇𝑙𝑖𝑚 (where 𝛥𝑇𝑙𝑖𝑚 is a small value), the filter is ordered not to act at all. This behavior is added to deal with the fact that cold-gas thrusters have a minimumwidth pulse that they can perform. More details on this are given in Section 4.4.1. 4.2. Linear model Sensors on-board the chaser estimate its absolute state 𝐫𝐼 𝐶,𝐯𝐼 𝐶and the relative state of the target with respect to the chaser-centered LVLH frame 𝐫𝐶 𝑇,𝐯𝐶 𝑇. The target is to be kept at a distance 𝑑𝐶𝑇 behind the chaser in the negative in-track direction, so 𝐫𝐶 𝑇 ,𝑟𝑒𝑓 = [0,−𝑑𝐶𝑇 ,0]𝑇, 𝐯𝐶 𝑇 ,𝑟𝑒𝑓 = [0,0,0]𝑇. Therefore, for the controller, the virtual state 𝐱= (𝐫𝐶 𝑉,𝐯𝐶 𝑉), is defined, where 𝐫𝐶 𝑉=𝐫𝐶 𝑇 ,𝑟𝑒𝑓 −𝐫𝐶 𝑇,𝐯𝐶 𝑉=𝐯𝐶 𝑇 ,𝑟𝑒𝑓 −𝐯𝐶 𝑇. The task of the controller is to regulate the virtual state 𝐱to the origin. The dynamics of the virtual state can be approximated by a modified version of the Hill–Clohessy–Wiltshire (HCW) equations:  𝐱=𝐴𝑐𝐱+𝐵𝑐𝐮+𝐰𝑐,(30) where 𝐴𝑐=⎡⎢⎢⎢⎢⎢⎢⎣ 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 3𝑛2+𝑐0 0 0 2𝑛0 0 −2𝑐0 −2𝑛0 0 0 0 −𝑛2+𝑐0 0 0 ⎤⎥⎥⎥⎥⎥⎥⎦ , 𝐵𝑐=⎡⎢⎢⎢⎢⎢⎢⎣ 000 000 000 100 010 001 ⎤⎥⎥⎥⎥⎥⎥⎦ , (31) and 𝐰𝑐is an unknown but bounded state disturbance. In 𝐴𝑐,𝑛is the mean motion of the chaser, which can be obtained from 𝐫𝐼 𝐶,𝐯𝐼 𝐶, and 𝑐is a new term that accounts for the gradient effect the ion beam has on the target when the chaser and the target are not properly Acta Astronautica 225 (2024) 689–706 696 J. Urrios et al. aligned (see [20]). This term 𝑐is computed assuming that the target is perfectly spherical (the so-called cannonball model), which is logically not the case, but allows the controller to have some kind of information about the effect of the ion beam on the system, enhancing performance. Coming from a cylinder, setting the radius of the apparent sphere is not straight-forward. In this work, it is chosen to be the mean between the cylinder radius and its semiheight, 𝑅𝑎𝑠 = (𝑅𝑐+𝐻𝑐∕2)∕2. This is a good compromise value for its simplicity and because, for cylindrical targets with height to radius ratios between 3.5 and 16.3, it lies between the equivalent sphere radius obtained by computing the mean frontal area of the cylinder (𝑅𝐴=√2𝑅𝑐∕𝜋(𝑅𝑐∕𝜋+ 2𝐻𝑐)), and the equivalent radius of a sphere with the same volume (𝑅𝑉=3 √3𝑅2 𝑐𝐻𝑐∕4). Another way to compute this equivalent sphere radius may be convenient for other 𝐻𝑐∕𝑅𝑐proportions. The system in (30), with 𝑛𝑥= 6 states and 𝑛𝑢= 3 control inputs, is an LTI continuous system, while the controller acts discretely. The discrete version of the system is 𝐱𝑘+1 =𝐴𝐱𝑘+𝐵𝐮𝑘+𝐰𝑘,(32) where, again, 𝐰𝑘is an unknown but bounded state disturbance. Assuming that the controller acts every 𝑇𝑐seconds, matrices 𝐴and 𝐵can be computed using the State Transition Matrix (STM) 𝛷(𝑡)as  𝛷=𝐴𝑐𝛷, 𝛷(0) = 𝐈6,(33) 𝐴=𝛷(𝑇𝑐),(34) 𝐵=∫𝑇𝑐 0 𝛷(𝜏)d𝜏 𝐵𝑐,(35) where 𝐈6is the 6 ×6 identity matrix. Although there are simple analytic solutions for 𝐴, 𝐵 in the classical HCW model, these solutions become much more cumbersome with the addition of the 𝑐term, and their numerical calculation becomes a better choice. However, (35) is not an efficient way of computing 𝐵, since it performs a quadrature of a function that performs a numerical integration each time it is evaluated. Instead, differentiating (35), applying the Fundamental Theorem of Calculus, and reordering the equations, one obtains the next system of Ordinary Differential Equations (ODEs):  𝛷=𝐴𝑐𝛷, 𝛷(0) = 𝐈6,(36)  𝛱=𝐴𝑐𝛱, 𝛱(0) = 𝟎6×3,(37)  𝛹=𝛱, 𝛹(0) = 𝐵𝑐,(38) in which 𝛷(𝑇𝑐) = 𝐴,𝛹(𝑇𝑐) = 𝐵. Note that, since the system is LTI, these matrices are calculated once and then stored; they do not need to be recomputed at each sampling time, except when the controller is re-tuned (since the chaser and target orbits slowly change, after some time the mean motion 𝑛must be updated and the controller re-tuned). 4.3. System constraints In accordance with the MPC formulation, the constraints that the system must meet at each time instant are expressed as polytopes. In the case of the state, a certain polytope that ensures full impact of the ion beam is defined for position, and maximum values are set for velocity. For the control, there are maximum values that the actuators can input into the system. 4.3.1. State constraints In terms of position, the chaser should not deviate significantly from the virtual reference, since that would imply only partially hitting the target (or potentially not hitting it at all), which in turn means a longer mission time and more fuel consumption, since the ion thrusters would be more time activated. Note, however, that it may not be convenient to impose very strict constraints in position either. The chaser has some clearance to move while still completely hitting the target (or, more Fig. 10. Visualization of the position constraints polytope. precisely, the cone that represents the 95% of the ion beam contribution completely hits the target). So, if the chaser is in that vicinity, it can be left unactuated while not draining mission performance; quite the contrary, it would increase efficiency because less fuel would be used. To compute the clearance region, it is assumed that within that region, the ion beam completely hits the target regardless of the target’s attitude. This is equivalent to saying that the ion beam must hit the sphere inscribed in the cylinder, that is, a sphere with radius 𝑅𝑐. With a divergence angle 𝛼0, recalling that in the nominal relative position the COMs of the target and chaser are separated by 𝑑𝐶𝑇 , and that the distance between the chaser’s COM and the ion beam source is 𝑑𝑉, the maximum displacement in the in-track direction while still entirely hitting the target can be computed to be 𝛥𝑦𝑚𝑎𝑥 =𝑅𝑐 sin 𝛼0 −(𝑑𝐶𝑇 −𝑑𝑉).(39) From there, if one moves towards the negative in-track direction with an angle less than 𝛼0, the target would still be completely hit, so this region is actually a cone. Another limit must be imposed on the other side, since the chaser should not get too close to the target. Assuming that in the worst attitude case the ion beam source should be at 𝑑𝑚𝑖𝑛 = 3 m of the target’s closest point and not less, the maximum movement in the negative in-track direction is 𝛥𝑦𝑚𝑖𝑛 = − ⎛⎜⎜⎝ 𝑑𝐶𝑇 −𝑑𝑉−√𝑅2 𝑐+(𝐻𝑐 2)2 −𝑑𝑚𝑖𝑛⎞⎟⎟⎠ .(40) For the simulation case studied in Section 5, this truncated cone is shown in Fig. 10, along with an inner approximation made with an eight-plane polytope, which is the actual region input into the controller. To build the inner approximation, the next approach was followed: •The cone vertex is denoted as 𝑃0={0, 𝛥𝑦𝑚𝑎𝑥,0}. •Eight points are defined at the minimum 𝑦distance, inscribing an octagon inside the cone intersection circle: 𝑃𝑖={(𝛥𝑦𝑚𝑎𝑥 −𝛥𝑦𝑚𝑖𝑛) tan 𝛼0cos 𝛽𝑖, 𝛥𝑦𝑚𝑖𝑛,(𝛥𝑦𝑚𝑎𝑥 −𝛥𝑦𝑚𝑖𝑛)tan 𝛼0sin 𝛽𝑖}, where 𝛽𝑖= 0,𝜋 4,𝜋 2, …,7𝜋 4for 𝑖= 1,…,8, and 𝑃9∶= 𝑃1. •Eight planes are defined with the three-point definition: 𝑃0𝑃𝑖𝑃𝑖+1, for 𝑖= 1,…,8. •Eight plane inequalities are set using those planes, such that all of them contain the origin. •The last inequality is given by 𝑦 > 𝛥𝑦𝑚𝑖𝑛. Position constraints are therefore given by 𝐴𝑟𝐫≤𝐛𝑟, where 𝐴𝑟and 𝑏𝑟are defined according to the polytope in Fig. 10. In terms of velocity, a maximum velocity 𝑣𝑚𝑎𝑥 = 1 m/min is imposed, to increase safety in operations. Thus, velocity restrictions are expressed as 𝐴𝑣𝐱≤𝐛𝑣, where 𝐴𝑣=[𝟎3𝐈3 𝟎3−𝐈3], 𝑏𝑣=𝑣𝑚𝑎𝑥𝟏6×1.(41) Acta Astronautica 225 (2024) 689–706 697 J. Urrios et al. Table 8 Total fuel consumption by radial CGTs, ion beams and out-of-plane CGTs. Duration [days] Fuel by radial CGTs [kg] Fuel by ion beams [kg] Fuel by out-of-plane CGTs [kg] 165.6 0.67 136.3 0.0003 Fig. 16. Performance of the controller for different initial relative states of the chaser. The chaser starts from different positions of a 1.5 m in radius circle in the CLVLH XY-plane (left) and YZ-plane (right). Red region denotes the position constraint polytope. Only the first orbit revolution is shown. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) control is recomputed (and thus sensor measurements are used), and excluding the departing asterisks, all the rest lie inside the red-shaded region, which also validates the boundings performed for the controller. As expected, symmetric behavior is seen in the YZ-plane, excluding small differences due to sensor and force noise and hitting the target from slightly different positions. 6. Conclusions The proposed Ion Beam Shepherd system with an enhanced actuator configuration, incorporating cold-gas thrusters alongside the typical Impulse Transfer and Impulse Compensation Thrusters, has proven to be a robust and viable solution for the deorbiting of large space debris such as the Zenit-2 second stage rocket body. Simulation results demonstrate that the system can successfully deorbit this statistically concerning 9000 kg object in approximately 165.6 days, with a fuel consumption of around 137 kg, less than 1% of which is attributed to the cold-gas thrusters. The inclusion of cold-gas thrusters together with a novel Model Predictive Controller ensure robustness without significant fuel penalties, making the system practical for near-future missions. The development of a fast and accurate ray-marching-based ion beam force and torque computation algorithm speeds up simulation times. While the results confirm the feasibility of using IBS for large debris deorbiting, further investigation is required on aspects such as higherfidelity ion beam models, ion beam sputtering effects on the chaser, long-term reliability of constant-thrust ion beams, or the use of state observers to loosen the sensoring requirements. Nevertheless, the study indicates that the IBS system is a promising candidate for future space debris removal missions. CRediT authorship contribution statement Javier Urrios: Software, Methodology, Conceptualization. Rafael Vazquez: Software, Methodology, Conceptualization. Francisco Gavilan: Software, Methodology, Conceptualization. Ignacio Alvarado: Software, Methodology, Conceptualization. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments The authors gratefully acknowledge support by grant TED2021- 132099B-C33 funded by MICIU/AEI/10.13039/501100011033 and by ‘‘European Union NextGenerationEU/PRTR’’. Appendix A. Abbreviations and nomenclature Vector quantities are represented as bold symbols. Except for Section 4.4.1, if a vector has a superscript, it indicates the reference frame in which it is measured, while the subscript indicates the object it is related to (𝐶for chaser and 𝑇for target); i.e., 𝐚𝐼 𝐶denotes the acceleration of the chaser measured from the ECI frame. The modulus of a vector is simply denoted as the vector itself unbold. For clarity, in the notation list, superscripts of vector quantities are omitted. ADR Active Debris Removal CGT Cold-Gas Thruster CLVLH Chaser-centered LVLH frame COM Center Of Mass DCM Direction Cosine Matrix ECEF Earth-Centered Earth-Fixed frame ECI Earth-Centered-Inertial frame HCW Hill-Clohessy-Wiltshire IBS Ion Beam Shepherd ICT Impulse Compensation Thruster ISA International Standard Atmosphere ITT Impulse Transfer Thruster LEO Low Earth Orbit LMI Linear Matrix Inequality Acta Astronautica 225 (2024) 689–706 704 J. Urrios et al. LQR Linear Quadratic Regulator LTI Linear Time-Invariant LVLH Local-Vertical Local-Horizontal frame MPC Model Predictive Control ODE Ordinary Differential Equation QP Quadratic Programming STM State Transition Matrix TBA Target Body Axes frame TLVLH Target-centered LVLH frame TRL Technology Readiness Level ZOH Zero-Order-Hold 𝟎𝑖×𝑗Zero matrix of size 𝑖×𝑗 𝟏𝑖×𝑗Ones matrix of size 𝑖×𝑗 𝐴𝑓,𝑖 Frontal area of object 𝑖 𝐚𝑖,𝑗 Acceleration of object 𝑖due to factor 𝑗 𝑐𝐷,𝑖 Drag coefficient of object 𝑖 𝑑𝐶𝑇 Distance between chaser and target COMs 𝑑𝑉Distance between chaser COM and ITT 𝐅𝐵Force exerted by the ICT on the target 𝐹𝐼𝐶𝑇 Force magnitude of ICT 𝐹𝐼𝑇 𝑇 Force magnitude of ITT 𝐻𝑐Target cylinder height ℎ𝑖Height of object 𝑖 𝐈𝑖Identity matrix of 𝑖th order 𝐼𝑠𝑝,𝐶𝐺 Specific impulse of CGTs 𝐼𝑠𝑝,𝐼𝐵 Specific impulse of ion beams 𝐌𝐵Torque exerted by the ICT on the target 𝑚𝑖Mass of object 𝑖 𝑁MPC control horizon 𝑛Instantaneous orbit mean motion 𝑛0Ion number density at reference plane center 𝐪𝐵 𝐴Attitude quaternion from A to B frame 𝑅0Radius of the 95% ion mass flow tube at the reference plane 𝑅𝐵 𝐴DCM from A to B frame 𝑅𝑐Target cylinder radius 𝐫𝑖Position of object 𝑖 𝑇𝑐Control refreshing period 𝐮Control variables vector 𝐯𝑖Velocity of object 𝑖 𝐯𝑟𝑖 Relative velocity between object 𝑖and air 𝐰𝑘State disturbance at instant 𝑘 𝐱Virtual state of the chaser 𝛼0Ion beam divergence angle 𝝎𝑇Target angular velocity Target’s inertia matrix Control constraints polytope State constraints polytope References [1] ESA Space Debris Office, ESA’s Annual Space Environment Report, 2023, https://www.sdo.esoc.esa.int/environment_report/Space_Environment_Report_ latest.pdf, (Accessed 27 September 2024). [2] D.J. Kessler, B.G. Cour-Palais, Collision frequency of artificial satellites: The creation of a debris belt, J. Geophys. Res. Space Phys. 83 (A6) (1978) 2637–2646, http://dx.doi.org/10.1029/JA083iA06p02637. [3] C.P. Mark, S. Kamath, Review of active space debris removal methods, Space Policy 47 (2019) 194–206, http://dx.doi.org/10.1016/j.spacepol.2018.12.005. [4] J. Locke, T.J. Colvin, L. Ratliff, A. Abdul-Hamid, C. Samples, Cost and benefit analysis of mitigating, tracking, and remediating orbital debris, Cost Benefit Anal. Mitig. Track. Remediat. Orbital Debris (2024). [5] S. Kitamura, Large space debris reorbiter using ion beam irradiation, in: 61 St International Astronautical Congress, Prague, Czech Republic, 2010. [6] C. Bonnal, J. Ruault, P. Bultel, M. Desjean, High level requirements for an operational space debris deorbiter, in: Proc. NASA/DARPA International Conference on Orbital Debris Removal, Chantilly, VA, 2009. [7] C. Bombardelli, J. Peláez, Ion beam shepherd for contactless space debris removal, J. Guid. Control Dyn. 34 (3) (2011) 916–920, http://dx.doi.org/10. 2514/1.51832. [8] S. Khoroshylov, Relative control of an ion beam shepherd satellite in eccentric orbits, Acta Astronaut. 176 (2020) 89–98, http://dx.doi.org/10.1016/j.actaastro. 2020.06.027. [9] V. Aslanov, A. Ledkov, Detumbling of axisymmetric space debris during transportation by ion beam shepherd in 3D case, Adv. Space Res. 69 (1) (2022) 570–580, http://dx.doi.org/10.1016/j.asr.2021.10.002. [10] M. Merino Martínez, E. Ahedo Galilea, C. Bombardelli, H. Urrutxua Cereijo, J. Peláez Álvarez, Hypersonic Plasma Plume Expansion in Space, 32nd International Electric Propulsion Conference, 2011. [11] M. Merino, F. Cichocki, E. Ahedo, A collisionless plasma thruster plume expansion model, Plasma Sources Sci. Technol. 24 (3) (2015) 035006, http: //dx.doi.org/10.1088/0963-0252/24/3/035006. [12] S. Kitamura, Y. Hayakawa, S. Kawamoto, A reorbiter for large GEO debris objects using ion beam irradiation, Acta Astronaut. 94 (2) (2014) 725–735, http://dx.doi.org/10.1016/j.actaastro.2013.07.037. [13] M. Merino, E. Ahedo, C. Bombardelli, H. Urrutxua, J. Peláez, Ion beam shepherd satellite for space debris removal, Prog. Propuls. Phys. 4 (2013) 789–802, http://dx.doi.org/10.1051/eucass/201304789. [14] F. Cichocki, M. Merino, E. Ahedo, M. Smirnova, A. Mingo, M. Dobkevicius, Electric propulsion subsystem optimization for ‘‘ion beam shepherd’’ missions, J. Propuls. Power 33 (2) (2017) 370–378, http://dx.doi.org/10.2514/1.B36105. [15] H. Urrutxua, C. Bombardelli, J.M. Hedo, A preliminary design procedure for an ion-beam shepherd mission, Aerosp. Sci. Technol. 88 (2019) 421–435, http: //dx.doi.org/10.1016/j.ast.2019.03.038. [16] A. Alpatov, F. Cichocki, A. Fokov, S. Khoroshylov, M. Merino, A. Zakrzhevskii, Determination of the force transmitted by an ion thruster plasma plume to an orbital object, Acta Astronaut. 119 (2016) 241–251, http://dx.doi.org/10.1016/ j.actaastro.2015.11.020. [17] M. Redka, S. Khoroshylov, Determination of the force impact of an ion thruster plume on an orbital object via deep learning, Space Sci. Technol. 28 (5) (2022) 15–26, http://dx.doi.org/10.15407/knit2022.05.015. [18] F. Cichocki, M. Merino, E. Ahedo, Spacecraft-plasma-debris interaction in an ion beam shepherd mission, Acta Astronaut. 146 (2018) 216–227, http://dx.doi.org/ 10.1016/j.actaastro.2018.02.030. [19] A. Melnikov, V. Abgaryan, A. Mogulkin, O. Peysakhovich, V. Svotina, Experimental study of ion beam interaction with target surface aimed at developing a contactless method for space debris removal by an ion beam, Acta Astronaut. 216 (2024) 120–128, http://dx.doi.org/10.1016/j.actaastro.2024.01.004. [20] C. Bombardelli, H. Urrutxua, M. Merino Martínez, E. Ahedo Galilea, J. Peláez Álvarez, Relative Dynamics and Control of an Ion Beam Shepherd Satellite, vol. 143, Advances in the Astronautical Sciences, 2012. [21] S. Khoroshylov, Relative motion control system of spacecraft for contactless space debris removal, Sci. Innov 14 (4) (2018) 5–16, http://dx.doi.org/10.15407/ scine14.04.005. [22] V.S. Aslanov, A.S. Ledkov, et al., Attitude motion of cylindrical space debris during its removal by ion beam, Math. Probl. Eng. 2017 (2017) http://dx.doi. org/10.1155/2017/1986374. [23] A. Ledkov, V. Aslanov, Attitude motion of space debris during its removal by ion beam taking into account atmospheric disturbance, in: J. Phys. Conf. Ser., 1050, IOP Publishing, 2018, 012041, http://dx.doi.org/10.1088/1742-6596/1050/1/ 012041. [24] V.S. Aslanov, A.S. Ledkov, Space debris attitude control during contactless transportation in planar case, J. Guid. Control Dyn. 43 (3) (2020) 451–461, http://dx.doi.org/10.2514/1.G004686. [25] A. Alpatov, S. Khoroshylov, C. Bombardelli, Relative control of an ion beam shepherd satellite using the impulse compensation thruster, Acta Astronaut. 151 (2018) 543–554, http://dx.doi.org/10.1016/j.actaastro.2018.06.056. [26] S. Khoroshylov, Out-of-plane relative control of an ion beam shepherd satellite using yaw attitude deviations, Acta Astronaut. 164 (2019) 254–261, http://dx. doi.org/10.1016/j.actaastro.2019.08.016. [27] A.S. Ledkov, V.S. Aslanov, Active space debris removal by ion multi-beam shepherd spacecraft, Acta Astronaut. 205 (2023) 247–257, http://dx.doi.org/10. 1016/j.actaastro.2023.02.003. [28] E. Camacho, C. Alba, Model Predictive Control, Springer, 2013. [29] F. Gavilan, R. Vazquez, E.F. Camacho, Chance-constrained model predictive control for spacecraft rendezvous with disturbance estimation, Control Eng. Pract. 20 (2) (2012) 111–122, http://dx.doi.org/10.1016/j.conengprac.2011.09.006. [30] R. Vazquez, F. Gavilan, E.F. Camacho, Model Predictive Control for Spacecraft Rendezvous in Elliptical Orbits with On/Off Thrusters, IFAC-PapersOnLine 48 (9) (2015) 251–256, http://dx.doi.org/10.1016/j.ifacol.2015.08.092. [31] M. Schwenzer, M. Ay, T. Bergs, D. Abel, Review on model predictive control: An engineering perspective, Int. J. Adv. Manuf. Technol. 117 (5) (2021) 1327–1349, http://dx.doi.org/10.1007/s00170-021-07682-3. [32] I. Alvarado, P. Krupa, D. Limon, T. Alamo, Tractable robust MPC design based on nominal predictions, J. Process Control 111 (2022) 75–85, http://dx.doi.org/ 10.1016/j.jprocont.2022.01.006. Acta Astronautica 225 (2024) 689–706 705 J. Urrios et al. [33] D. Limón, I. Alvarado, T. Alamo, E.F. Camacho, MPC for tracking piecewise constant references for constrained linear systems, Automatica 44 (9) (2008) 2382–2387, http://dx.doi.org/10.1016/j.automatica.2008.01.023. [34] D. Limón, I. Alvarado, T. Alamo, E.F. Camacho, Robust tube-based MPC for tracking of constrained linear systems with additive disturbances, J. Process Control 20 (3) (2010) 248–260, http://dx.doi.org/10.1016/j.jprocont.2009.11. 007. [35] G. Borelli, M. Trisolini, M. Massari, C. Colombo, et al., A comprehensive ranking framework for active debris removal missions candidates, in: 8th European Conference on Space Debris, ESA/ESOC, ESA, 2021, pp. 1–10. [36] J. Wong, Ray marching and signed distance functions, 2016, https://jamiewong.com/2016/07/15/ray-marching-signed-distance-functions/, (Accessed 27 September 2024). [37] Space-Track.org, Orbital elements data, 2024, https://www.space-track.org, (Accessed 27 September 2024). [38] K. Anflo, B. Crowe, In-space demonstration of an ADN-based propulsion system, in: 47th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, 2011, p. 5832, http://dx.doi.org/10.2514/6.2011-5832. [39] W. Tighe, K.-R. Chien, E. Solis, P. Rebello, D. Goebel, J. Snyder, Performance evaluation of the XIPS 25-cm thruster for application to NASA discovery missions, in: 42nd AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, 2006, p. 4666, http://dx.doi.org/10.2514/6.2006-4666. Acta Astronautica 225 (2024) 689–706 706