Simulating Space Debris Removal Maneuvers: A Multi-Orbit Approach for Active Debris Clearance in LEO and GEO

1. Introduction

The ever-growing accumulation of space debris—comprising defunct satellites, discarded rocket stages, and fragments from past collisions or explosions—poses an escalating threat to the safety and sustainability of current and future space operations. As of May 2022, there were over 130 million pieces of debris larger than 1 mm in Earth’s orbit. Research indicates that the debris population in low Earth orbit (LEO) is unstable and will continue to grow, even without new launches, due to collisions between existing debris. To maintain a sustainable LEO environment, the removal of at least five large debris objects annually is deemed essential. Active debris removal (ADR) has thus become a critical focus for space agencies, with several methods being explored, including robotic arms, flying nets, tethers, harpoons, lasers, ion beam shepherds, and electrodynamic tethers.

A major challenge in ADR is the detumbling of space debris, which often rotates or tumbles due to residual angular momentum and space disturbances. Tumbling debris complicates capture and deorbit operations, leading to potential risks such as entanglement of tether systems or damage to robotic arms. To address this, a variety of detumbling techniques have been proposed, ranging from biologically inspired mechanisms like gecko adhesive grippers and soft tentacles to methods utilizing magnetic eddy currents and thruster plume impingement. Despite their promise, these methods face significant limitations. For example, contact-based methods often result in debris moving away from the detumbling equipment due to reactive forces, necessitating frequent positional adjustments by the servicing satellite, which increases fuel consumption and reduces operational life. Furthermore, most existing methods require large and complex systems, making them dependent on major space agencies and expensive to implement.

To overcome these challenges, a novel nanosatellite-based detumbling approach is proposed, utilizing adhesion and harpoon-based mechanisms. Unlike conventional methods, this strategy employs a nanosatellite equipped with multiple harpoons and thrusters. Deployed from a servicing satellite, the nanosatellite attaches itself to the debris by embedding harpoons into the target’s surface. This robust attachment allows the nanosatellite’s thrusters to detumble the debris effectively without the risk of drifting away. The harpoon-based attachment method also avoids generating new debris, as tests have shown that any fragments created remain contained within the target.

In this study, we build on previous work by developing a high-fidelity dynamic model for the combined system of the nanosatellite and the debris object. This model accounts for coupled orbital and attitude dynamics during the detumbling process, incorporating factors such as atmospheric drag, gravity-gradient torques, and variations in inertial parameters due to fuel consumption. An adaptive sliding mode control scheme, integrated with pulse-width pulse-frequency (PWPF) modulation, is proposed to handle uncertainties in the model and ensure robust control of the debris object. Unlike previous approaches, this method does not require precise knowledge of the debris’ inertial parameters, making it highly effective for detumbling operations without extensive prior assessment.

The remainder of this paper is structured as follows: Section 2 provides an overview of the dynamic model and preliminary considerations; Section 3 details the collision modeling between the detumbling device and debris; Section 4 derives the coupled orbit-attitude dynamics of the system; Section 5 introduces the adaptive sliding mode control approach; Section 6 presents simulation results demonstrating the effectiveness of the proposed strategy; and Section 7 concludes with a discussion of the findings and future work.

2. Preliminaries

Figure 1. An example of the adhesion detumbling approach based on harpoons.

Figure 1. An example of the adhesion detumbling approach based on harpoons.

2.1. System Overview and Reference Frames

This study presents a detumbling strategy for stabilizing tumbling space debris, such as a rocket upper stage in low Earth orbit (LEO). The proposed method involves a detumbling device that secures itself to the debris by embedding multiple harpoons into its surface, creating a single, rigid system. The dynamics of this system are described using several reference frames:

• Inertial Frame (${\Sigma }_I$): A fixed frame centered at Earth’s center, used for describing orbital motion.
• Target Frame (${\Sigma }_t$): Attached to the center of mass (CoM) of the debris, describing its rotation.
• Device Frame (${\Sigma }_d$): Attached to the detumbling device’s CoM, used for thruster control.
• Combined System Frame (${\Sigma }_s$): Represents the newly formed system after attachment, aligned with the target frame.

The detumbling device is equipped with thrusters oriented along specific axes of the device frame, enabling controlled maneuvers to reduce the debris’s angular velocity. For simplicity, both the device and debris are modeled as rigid bodies, and thruster operations are considered instantaneous and precise.

2.2. Notations

Figure 2. Position vectors and reference frames for the combined system.

Figure 2. Position vectors and reference frames for the combined system.

In this paper, various symbolic representations are used to describe vectors, matrices, and transformations between different reference frames. To ensure clarity, the key notations are defined as follows:

For any vector $\mathbf{a}$ and matrix $\mathbf{A}$ expressed in the inertial reference frame ${\Sigma }_I$, their representations in another reference frame ${\Sigma }_i$ are denoted by ⁱ$\mathbf{a}$ and ⁱ$\mathbf{A}$, respectively. The relationship between the frames is described by a rotation matrix. Specifically, the rotation matrix that transforms coordinates from frame ${\Sigma }_i$ to frame ${\Sigma }_j$ is represented as ^j${\mathbf{R}}_i$. For simplicity, if ${\Sigma }_j$ is the inertial frame ${\Sigma }_I$, the rotation matrix is simply written as ${\mathbf{R}}_i$.

The transformation rules for converting vectors and matrices between frames are defined as:

$${}^j\mathbf{a}{= }^j{{\mathbf{R}}_i}^i\mathbf{a}{, }^j\mathbf{A}{= }^j{{\mathbf{R}}_i}^i\mathbf{A}{{(}^j{\mathbf{R}}_i)}^T.$$

For three-dimensional vectors, the skew-symmetric matrix representation is an important tool for expressing cross products. Given a vector $\mathbf{a}={[a_1,a_2,a_3]}^T$, the skew-symmetric matrix ${\mathbf{a}}^{\times }$ is defined as:

$${\mathbf{a}}^{\times }=\left[\begin{array}{ccc}0& -a_3& a_2\\ a_3& 0& -a_1\\ -a_2& a_1& 0\end{array}\right].$$

This matrix is used to compute the cross product between two vectors. For any vectors $\mathbf{b},\mathbf{c}\in {\mathbb{R}}^3$, the cross product can be represented using the skew-symmetric matrix as $\mathbf{b}\times \mathbf{c}={\mathbf{b}}^{\times }\mathbf{c}$. Note that the cross product is anti-commutative, meaning $\mathbf{b}\times \mathbf{c}=-\mathbf{c}\times \mathbf{b}$.

3. Collision Model

To ensure a successful attachment of the detumbling device to the debris, the harpoons must collide with the target at a high relative velocity. This high-impact collision can induce rapid changes in the motion of the target. However, given the extremely short penetration time of the harpoons (approximately 2 ms), any change in the target’s position and orientation during the collision is negligible. Therefore, we assume that the positions (${\mathbf{r}}_t^{+}$ and ${\mathbf{r}}_t^{-}$) and orientations (${\varphi }_t^{+}$ and ${\varphi }_t^{-}$) of the target before and after the collision are nearly identical. Additionally, external forces acting on the target, such as gravitational and space disturbance forces, are considered insignificant during this short interaction period.

Since the collision force between the detumbling device and the target is internal to the system, the total momentum of the system is conserved during the collision. This principle of momentum conservation allows us to compute the velocity of the combined system after the collision based on the initial velocities of the detumbling device and the target.

The linear momentum of the system before the collision can be expressed as:

$${\mathbf{P}}^{-}=m_t{\mathbf{v}}_t^{-}+m_d{\mathbf{v}}_d^{-}$$

where $m_t$ and $m_d$ are the masses of the target and the detumbling device, respectively, and ${\mathbf{v}}_t^{-}$ and ${\mathbf{v}}_d^{-}$ are their velocities before the collision.

The angular momentum of the system with respect to the target’s center of mass (CoM) before the collision is given by:

$${\mathbf{L}}^{-}={\mathbf{I}}_t{\mathit{\omega }}_t^{-}+{\mathbf{I}}_d{\mathit{\omega }}_d^{-}+m_d\ae \times {\mathbf{v}}_d^{-}$$

where ${\mathbf{I}}_t$ and ${\mathbf{I}}_d$ are the moments of inertia of the target and the detumbling device, respectively, ${\mathit{\omega }}_t^{-}$ and ${\mathit{\omega }}_d^{-}$ are their angular velocities, and $\ae$ is the vector from the target’s CoM to the detumbling device’s CoM.

After the collision, the detumbling device becomes rigidly attached to the target, implying that their angular velocities are the same (${\mathit{\omega }}_d^{+}={\mathit{\omega }}_t^{+}$). Consequently, the linear velocity of the detumbling device can be described as:

$${\mathbf{v}}_d^{+}={\mathbf{v}}_t^{+}+{\mathit{\omega }}_t^{+}\times \ae ={\mathbf{v}}_t^{+}-{\ae }^{\times }{\mathit{\omega }}_t^{+}$$

where ${\ae }^{\times }$ denotes the skew-symmetric matrix of $\ae$.

The linear and angular momentum of the combined system after the collision are then given by:

$${\mathbf{P}}^{+}=m_s{\mathbf{v}}_t^{+}-m_d{\ae }^{\times }{\mathit{\omega }}_t^{+}$$

$${\mathbf{L}}^{+}=({\mathbf{I}}_t+{\mathbf{I}}_d-m_d{\ae }^{\times }{\ae }^{\times }){\mathit{\omega }}_t^{+}+m_d{\ae }^{\times }{\mathbf{v}}_t^{+}$$

where $m_s=m_t+m_d$ is the total mass of the combined system.

Using the conservation of linear and angular momentum (i.e., ${\mathbf{P}}^{+}={\mathbf{P}}^{-}$ and ${\mathbf{L}}^{+}={\mathbf{L}}^{-}$), we can derive the velocities of the combined system after the collision:

$$\left[\begin{array}{c}{\mathbf{v}}_t^{+}\\ {\mathit{\omega }}_t^{+}\end{array}\right]={\mathbf{M}}^{-1}\left[\begin{array}{c}{\mathbf{P}}^{-}\\ {\mathbf{L}}^{-}\end{array}\right]$$

where $\mathbf{M}\in {\mathbb{R}}^{6\times 6}$ is defined as:

$$\mathbf{M}=\left[\begin{array}{cc}{m}_s{\mathbf{E}}_{3\times 3}& -m_d{\ae }^{\times }\\ m_d{\ae }^{\times }& {\mathbf{I}}_t+{\mathbf{I}}_d-m_d{\ae }^{\times }{\ae }^{\times }\end{array}\right].$$

It can be shown that the matrix $\mathbf{M}$ is positive definite, ensuring that its inverse always exists. This formulation provides a more stable and reliable approach for calculating the post-collision velocities compared to previous methods, as it reduces the risk of numerical instability caused by an ill-conditioned matrix.

4. Dynamic Model

4.1. Dynamic Equations of the Combined System

In this section, we develop a high-precision coupled dynamics model to describe the orbital and rotational motion of the combined system formed after the detumbling device attaches to the target debris. To simplify the notation, we omit the superscript “+” for state vectors describing the system after the collision. Given the assumption that the change in the center of mass (CoM) due to fuel consumption of the thrusters is negligible, the position vector of the combined system’s CoM can be defined as:

$${\mathbf{r}}_s={\displaystyle \frac{m_t{\mathbf{r}}_t+m_d({\mathbf{r}}_t+\mathit{\rho })}{m_t+m_d}}={\mathbf{r}}_t+{\mathbf{b}}_t$$

where ${\mathbf{b}}_t={\displaystyle \frac{m_d}{m_s}}\mathit{\rho }$ and $m_s=m_t+m_d$ is the total mass of the combined system. Because the frames ${\Sigma }_s$ and ${\Sigma }_t$ are parallel, the angular velocities are identical (${\mathit{\omega }}_s={\mathit{\omega }}_t$), and the linear velocity of the system can be expressed as:

$${\mathbf{v}}_s={\mathbf{v}}_t+{\mathit{\omega }}_s\times {\mathbf{b}}_t={\mathbf{v}}_t-{\mathbf{b}}_t^{\times }{\mathit{\omega }}_s$$

To avoid singularities in attitude representation, we use the unit quaternion ${\mathbf{q}}_s={\left[{\mathit{ϵ}}_s^T{\eta }_s\right]}^T$ to represent the orientation of the combined system. The rotation matrix ${\mathbf{R}}_s$ corresponding to this quaternion is given by:

$${\mathbf{R}}_s=(2{\eta }_s^2-1){\mathbf{E}}_{3\times 3}+2{\mathit{ϵ}}_s{\mathit{ϵ}}_s^T-2{\eta }_s{\mathit{ϵ}}_s^{\times }$$

The attitude kinematics equation for the combined system is:

$${\dot{\mathbf{q}}}_s={\displaystyle \frac{1}{2}}\left[\begin{array}{c}{\eta }_s{\mathbf{E}}_{3\times 3}+{\mathit{ϵ}}_s^{\times }\\ -{\mathit{ϵ}}_s^T\end{array}\right]{\mathit{\omega }}_s$$

Using the Newton-Euler approach, the coupled translational and rotational dynamics of the combined system in the inertial frame ${\Sigma }_I$ are given by:

$$m_s{\dot{\mathbf{v}}}_s={\mathbf{f}}_{\mathrm{ext}}+{\mathbf{f}}_{\mathrm{con}}$$

$${\mathbf{I}}_s{\dot{\mathit{\omega }}}_s=-{\mathit{\omega }}_s^{\times }{\mathbf{I}}_s{\mathit{\omega }}_s+{\mathit{\tau }}_{\mathrm{ext}}+{\mathit{\tau }}_{\mathrm{con}}$$

where ${\mathbf{I}}_s$ is the inertia matrix of the combined system, ${\mathbf{f}}_{\mathrm{ext}}$ and ${\mathit{\tau }}_{\mathrm{ext}}$ are external forces and torques, and ${\mathbf{f}}_{\mathrm{con}}$ and ${\mathit{\tau }}_{\mathrm{con}}$ are control forces and torques generated by the thrusters. The inertia matrix ${\mathbf{I}}_s$ is defined as:

$${\mathbf{I}}_s={\mathbf{I}}_t+{\mathbf{I}}_d-m_t{\mathbf{b}}_t^{\times }{\mathbf{b}}_t^{\times }-m_d{\mathbf{b}}_d^{\times }{\mathbf{b}}_d^{\times }$$

where ${\mathbf{b}}_d=\mathit{\rho }-{\mathbf{b}}_t$. To simplify further analysis, the attitude dynamics in the combined frame ${\Sigma }_s$ are expressed as:

$${\mathbf{I}}_s^s{\dot{\mathit{\omega }}}_s^s=-{\left({\mathit{\omega }}_s^s\right)}^{\times }{\mathbf{I}}_s^s{\mathit{\omega }}_s^s+{\mathit{\tau }}_{\mathrm{ext}}^s+{\mathit{\tau }}_{\mathrm{con}}^s$$

where ${\mathbf{I}}_s^s={\mathbf{R}}_s{\mathbf{I}}_s{\mathbf{R}}_s^T$, ${\dot{\mathit{\omega }}}_s^s={\mathbf{R}}_s{\dot{\mathit{\omega }}}_s$, and ${\mathit{\tau }}_{\mathrm{ext}}^s={\mathbf{R}}_s{\mathit{\tau }}_{\mathrm{ext}}$, ${\mathit{\tau }}_{\mathrm{con}}^s={\mathbf{R}}_s{\mathit{\tau }}_{\mathrm{con}}$.

Figure 3. An atmospheric drag schematic diagram.

Figure 3. An atmospheric drag schematic diagram.

4.2. External Forces and Torques

The primary external forces affecting space objects in low Earth orbit (LEO) are atmospheric drag and the J2 gravitational perturbation due to Earth’s oblateness. Thus, the external force ${\mathbf{f}}_{\mathrm{ext}}$ is given by:

$${\mathbf{f}}_{\mathrm{ext}}={\mathbf{f}}_{\mathrm{gra}}+{\mathbf{f}}_{J2}+{\mathbf{f}}_{\mathrm{atm}}$$

where ${\mathbf{f}}_{\mathrm{gra}}$ represents the gravitational force, ${\mathbf{f}}_{J2}$ is the J2 perturbative force, and ${\mathbf{f}}_{\mathrm{atm}}$ is the atmospheric drag force. The gravitational force is computed as:

$$\begin{array}{c}\hfill {\mathbf{f}}_{\mathrm{gra}}=-{\displaystyle \frac{m_s{\mu }_E}{r_s^3}}{\mathbf{r}}_s+{\displaystyle \frac{3{\mu }_E}{2r_s^5}}\left({\mathbf{I}}_s+\mathrm{tr}\left({\mathbf{I}}_s\right){\mathbf{E}}_{3\times 3}\right){\mathbf{r}}_s\\ \hfill -{\displaystyle \frac{15{\mu }_E\left({\mathbf{r}}_s^T{\mathbf{I}}_s{\mathbf{r}}_s\right)}{2r_s^7}}{\mathbf{r}}_s\end{array}$$

(1)

The J2 perturbative force is expressed as:

$$\begin{array}{c}\hfill {\mathbf{f}}_{J2}={\displaystyle \frac{3m_s{\mu }_E{J}_E{R}_E^2}{2r_s^5}}\left(\left(5{\displaystyle \frac{{\left(z_I^T{\mathbf{r}}_s\right)}^2}{r_s^2}}-1\right){\mathbf{r}}_s\right.\\ \hfill \left.-2\left(z_I^T{\mathbf{r}}_s\right)z_I\right)\end{array}$$

(2)

where $J_E$ is the J2 perturbation coefficient, $R_E$ is Earth’s equatorial radius, and $z_I={[0,0,1]}^T$ is the unit vector in the inertial frame.

Atmospheric drag depends on the shape, altitude, and velocity relative to the atmosphere. If the debris is approximated as a cylinder, and the detumbling device is considered negligible in size, the atmospheric drag is calculated as:

$${\mathbf{f}}_{\mathrm{atm}}=-{\displaystyle \frac{1}{2}}{\rho }_{\mathrm{atm}}{C}_{\mathrm{atm}}S\parallel {\mathbf{v}}_{\mathrm{rel}}\parallel {\mathbf{v}}_{\mathrm{rel}}$$

where ${\rho }_{\mathrm{atm}}$ is the atmospheric density, $C_{\mathrm{atm}}$ is the drag coefficient, S is the cross-sectional area, and ${\mathbf{v}}_{\mathrm{rel}}$ is the relative velocity between the combined system and the atmosphere.

The total external torque ${\mathit{\tau }}_{\mathrm{ext}}$ acting on the system includes contributions from gravity-gradient torque and atmospheric drag torque:

$${\mathit{\tau }}_{\mathrm{ext}}={\mathit{\tau }}_{\mathrm{gra}}+{\mathit{\tau }}_{\mathrm{atm}}$$

where the gravity-gradient torque is:

$${\mathit{\tau }}_{\mathrm{gra}}={\displaystyle \frac{3{\mu }_E}{r_s^5}}{\mathbf{r}}_s^{\times }{\mathbf{I}}_s{\mathbf{r}}_s$$

and the atmospheric drag torque is:

$${\mathit{\tau }}_{\mathrm{atm}}=\sum _{i=1}^3{\mathbf{p}}_i^{\times }{\mathbf{f}}_{\mathrm{atm},i}$$

with ${\mathbf{p}}_i$ denoting the vector from the CoM to the center of pressure of each face.

5. Controller Design

5.1. Adaptive Sliding Mode Control

Space debris typically exhibits non-cooperative behavior with uncertain inertial parameters, making it challenging to precisely control its detumbling. A robust control strategy is required to handle these parametric uncertainties and external disturbances. Sliding Mode Control (SMC) is effective for such systems with high nonlinearity and modeling uncertainties and has been widely applied in fields like robotics, aerospace, and industrial control. However, conventional SMC approaches require a large control gain based on prior knowledge of disturbance boundaries to maintain stability. If the gain is not properly tuned, this can lead to severe system chattering.

To overcome these limitations, we propose an adaptive sliding mode control scheme that improves detumbling performance while minimizing system chattering. The detumbling device is equipped with five thrusters capable of providing constant thrust $f_d$ along the $\pm x_d$, $\pm y_d$, and $+z_d$ directions. We define a continuous control input vector $\mathbf{u}={[u_x,u_y,u_z]}^T$, where $-f_d\le u_x,u_y\le f_d$ and $0\le u_z\le f_d$. The control force and torque vectors can then be expressed as:

$${\mathbf{f}}_{\mathrm{con}}={\mathbf{R}}_d\mathbf{u}={\mathbf{R}}_I^{-1}{\mathbf{R}}_d\mathbf{u},$$

$${\mathit{\tau }}_{\mathrm{con}}={\mathbf{R}}_I{\mathbf{b}}_d^{\times }\mathbf{u}={\mathbf{b}}_d^{\times }{\mathbf{R}}_d\mathbf{u}.$$

Given the limited thrust capacity, input saturation must be considered. To facilitate the controller design, the dynamic equation from Equation (15) is reformulated as:

$${\mathbf{I}}_s{\dot{\mathit{\omega }}}_s=-{\mathit{\omega }}_s^{\times }{\mathbf{I}}_s{\mathit{\omega }}_s+\mathrm{sat}\left({\mathit{\tau }}_{\mathrm{con}}\right)+\mathbf{d},$$

where $\mathrm{sat}\left({\mathit{\tau }}_{\mathrm{con}}\right)$ denotes the saturated control torque, and $\mathbf{d}$ represents the lumped uncertainty caused by external disturbances and model inaccuracies.

To handle uncertainties in the control torque due to the position uncertainty of the system’s CoM, we express ${\mathit{\tau }}_{\mathrm{con}}$ as:

$${\mathit{\tau }}_{\mathrm{con}}={\mathbf{b}}_d\times {\mathbf{R}}_d\mathbf{u}=\widehat{{\mathbf{b}}_d}\times {\mathbf{R}}_d\mathbf{u}+\widetilde{{\mathbf{b}}_d}\times {\mathbf{R}}_d\mathbf{u},$$

where the hat (^) denotes estimation, and tilde (˜) denotes uncertainty. The bounds on $\mathrm{sat}\left({\mathit{\tau }}_{\mathrm{con}}\right)$ can be determined using $\widehat{{\mathbf{b}}_d}$ and $\mathbf{u}$:

$$\mathrm{sat}\left({\mathit{\tau }}_{\mathrm{con}}\right)\in [{\mathit{\tau }}_{\mathrm{con}}^{\min },{\mathit{\tau }}_{\mathrm{con}}^{\max }].$$

To manage control constraints, we introduce a new control input $\mathbf{w}={[w_1,w_2,w_3]}^T$ with $\mathrm{sat}\left(\mathbf{w}\right)\in [-1,1]$:

$$\mathrm{sat}\left({\mathit{\tau }}_{\mathrm{con}}\right)=\mathbf{A}\mathrm{sat}\left(\mathbf{w}\right)+\mathbf{h},$$

where $\mathbf{A}$ is a diagonal matrix and $\mathbf{h}$ is a constant vector. The adaptive sliding mode control law is given by:

$$\mathbf{w}=-{\mathbf{K}}_1\sigma -{\mathbf{K}}_2\mathrm{sgn}\left(\sigma \right)-\mathit{\xi },$$

where $\sigma =k{\mathit{\omega }}_s$ is the sliding surface, ${\mathbf{K}}_1$ and ${\mathbf{K}}_2$ are gain matrices, and $\mathit{\xi }$ is the adaptive term. The updating laws for the adaptive gains are defined to ensure convergence and stability based on Lyapunov stability theory.

5.2. PWPF Modulation

Figure 4. A block diagram showing the suggested control plan.

Figure 4. A block diagram showing the suggested control plan.

Figure 5. PWPF modulator.

Figure 5. PWPF modulator.

To convert the continuous control signal $\mathbf{w}$ into a format suitable for on-off thruster operation, Pulse Width Pulse Frequency (PWPF) modulation is employed. PWPF modulation is chosen for its near-linear operation, high efficiency, noise immunity, and tunable pulse width and frequency. The reference control signal ${\mathbf{u}}_{\mathrm{ref}}$ is derived using the pseudo-inverse of $\widehat{{\mathbf{b}}_d}\times$:

$${\mathbf{u}}_{\mathrm{ref}}={\mathbf{R}}_d^{-1}\mathrm{pinv}(\widehat{{\mathbf{b}}_d}\times )\mathrm{sat}\left({\mathit{\tau }}_{\mathrm{con}}\right).$$

The PWPF modulator, as shown in the control block diagram, includes a first-order lag filter, a Schmitt trigger, and a feedback loop. The filter output is:

$${\dot{\mathbf{u}}}_{\mathrm{fil}}=-{\displaystyle \frac{1}{T_m}}{\mathbf{u}}_{\mathrm{fil}}+{\displaystyle \frac{K_m}{T_m}}({\mathbf{u}}_{\mathrm{pwpf}}-{\mathbf{u}}_{\mathrm{ref}}),$$

where $T_m$ is the time constant, and $K_m$ is the gain. The Schmitt trigger switches the thrusters on or off based on the filtered signal, converting the continuous control forces and torques into discrete on-off commands for precise thruster control.

The proposed adaptive sliding mode controller combined with PWPF modulation provides robust, efficient control of the detumbling process for non-cooperative space debris, enhancing stability while minimizing fuel consumption and chattering.

6. Simulation Results

This section presents the numerical simulations performed to analyze the coupled orbit-attitude dynamics of the system and verify the effectiveness of the proposed adaptive sliding mode control. Table 1 summarizes the initial conditions and inertial parameters of the target and the detumbling device.

The orbital elements of the target related to the vectors ${\mathbf{r}}_t$ and ${\mathbf{v}}_t^{-}$ are: semi-major axis of 6878.1363 km, eccentricity of 0.05, inclination of ${10}^{\circ }$, right ascension of ascending node $0^{\circ }$, argument of perigee $0^{\circ }$, and true anomaly ${15}^{\circ }$. The attitude of the detumbling device is defined by Z-Y-X Euler angles, and the associated rotation matrix is:

$${\mathbf{R}}_d\left({\mathit{\varphi }}_d\right)=\left[\begin{array}{ccc}0& 0.9962& -0.0872\\ 0.0872& 0.0868& 0.9924\\ 0.9962& -0.0076& -0.0868\end{array}\right].$$

This results in the relative position vector $\mathit{\rho }={\mathbf{a}}_t+{\mathbf{R}}_d{\mathbf{a}}_d=\left[\begin{array}{c}-0.1980\\ -0.7585\\ 0.3174\end{array}\right]$ m.

The relative velocity between the detumbling device and the target in the detumbling device’s frame ${\Sigma }_d$ is given by $\Delta {\mathbf{v}}_d^{-}$, and the linear velocity of the detumbling device before collision is:

$${\mathbf{v}}_d^{-}={\mathbf{v}}_t^{-}+{\mathbf{R}}_t{\mathbf{R}}_d\Delta {\mathbf{v}}_d^{-}=\left[\begin{array}{c}-1.9956\\ 7.6440\\ 1.3376\end{array}\right]\mathrm{km}/\mathrm{s}.$$

Using these values, the linear and angular velocities of the combined system are computed as:

$${\mathbf{v}}_s\left(0\right)=\left[\begin{array}{c}-1.9733\\ 7.6263\\ 1.3445\end{array}\right]\mathrm{km}/\mathrm{s},{\mathit{\omega }}_s\left(0\right)=\left[\begin{array}{c}20.3921\\ 21.4282\\ 9.2391\end{array}\right]\mathrm{deg}/\mathrm{s}.$$

The initial inertia matrix for the combined system is:

$${\mathbf{I}}_s\left(0\right)=\left[\begin{array}{ccc}526.8332& -1.5972& 0.6305\\ -1.5972& 521.3707& 2.2244\\ 0.6305& 2.2244& 186.4078\end{array}\right]\mathrm{kg}·{\mathrm{m}}^2.$$

The collision primarily caused a significant change in angular velocity, necessitating careful consideration in the control scheme. The robustness of the control strategy is tested by assuming $\mathbf{s}{\widehat{\mathbf{b}}}_d=\left[\begin{array}{c}-0.2732\\ -0.8156\\ 0.2341\end{array}\right]$ m and an adjusted inertia matrix:

$$\mathbf{I}\widehat{\mathbf{s}}=\left[\begin{array}{ccc}790.2498& -4.7915& 1.2609\\ -4.7915& 938.4673& 5.5610\\ 1.2609& 5.5610& 149.1262\end{array}\right]\mathrm{kg}·{\mathrm{m}}^2.$$

For control, the parameters are: $k=1$, $K_1=100\mathbf{I}$, $K_2=0.1\mathbf{I}$, $\beta =10$, $\epsilon =0.01$, $\widehat{\mu }\left(0\right)=100$, $\widehat{\gamma }\left(0\right)=1.0\times {10}^{-5}$, and PWPF modulator settings are $K_m=4.0$, $T_m=0.8$, $U_m=0.10$, $U_{\mathrm{off}}=0.05$, $U_{\mathrm{on}}=0.20$.

The simulation shows that the angular velocities of the system decrease rapidly. After 4000 seconds, they drop to approximately $\left[\begin{array}{c}0.0113\\ 0.1889\\ 0.0148\end{array}\right]$ deg/s, and continue to decrease, indicating the control scheme’s effectiveness against uncertainties and disturbances.

The detumbling mission is completed within one orbital period, with minor changes to the orbit’s perigee and apogee. The fuel consumed for the entire detumbling process is about 0.20 kg, demonstrating that the proposed control method is efficient for handling large debris objects in space.

7. Results and Discussion

The computational aspect of the simulation was performed using numerical integration techniques to model the dynamics of the CubeSat-debris system. Specifically, the simulation was carried out using the Runge-Kutta 4th order (RK4) method with a time step of 0.1 seconds. This method was chosen for its balance between accuracy and computational efficiency, allowing for the precise tracking of the CubeSat’s attitude and motion in response to various external torques, such as gravity-gradient and magnetic disturbances. The integration accounted for both control and disturbance torques to evaluate the robustness of the proposed adaptive controller under realistic space conditions.

Figure 6. Trajectories of debris (blue), attachment point (dashed), and CubeSat (red) in inertial frame $F_i$, highlighting key events during the rendezvous and attachment phases.

Figure 6. Trajectories of debris (blue), attachment point (dashed), and CubeSat (red) in inertial frame $F_i$, highlighting key events during the rendezvous and attachment phases.

Figure 7. CubeSat’s trajectory relative to tumbling debris in body frame $F_2$, showing its maneuver around the debris during rendezvous.

Figure 7. CubeSat’s trajectory relative to tumbling debris in body frame $F_2$, showing its maneuver around the debris during rendezvous.

The Deorbiter CubeSat starts with zero angular velocity ($\mathit{\omega }\left(t_0\right)=0$) after stabilizing post-release. It initiates a rendezvous maneuver with tumbling debris at a starting distance of about 34.6 meters. The debris follows a nearly circular equatorial orbit, but its tumbling motion complicates the CubeSat’s path. The CubeSat completes its synchronization and rendezvous around the debris in approximately 176.7 minutes. In the debris body frame, the CubeSat’s path appears complex due to the debris’s rotation. After attachment, the CubeSat stabilizes and maintains a fixed relative position, crucial for deorbiting maneuvers.

Figure 8. An exemplar reaction wheel momentum unloading process employing magnetorquers.

Figure 8. An exemplar reaction wheel momentum unloading process employing magnetorquers.

The time histories of command magnetic torques and reaction wheel momentum for both proportional and bang-bang controllers reveal key differences in desaturation performance. With wheel saturation, the proportional controller takes around three hours to fully desaturate, while the bang-bang controller achieves this in just about one hour, though it draws maximum current from onboard batteries. Given an orbital period of 1.62 hours, the bang-bang controller can complete desaturation in just over half an orbit, offering a faster but more power-intensive solution.

Figure 9. Using the time-optimal detumbling controller, the debris body rates profile during the detumbling movement.

Figure 9. Using the time-optimal detumbling controller, the debris body rates profile during the detumbling movement.

The body rate profiles of KOMPSAT-1 and Orion 50 XL show a reduction in debris body rate magnitude, $\parallel \omega \parallel$, during detumbling. KOMPSAT-1 achieves detumbling in about 15 days, while Orion 50 XL takes 9 days. However, the full detumbling maneuver includes momentum unloading from reaction wheels, which takes about 25 days for KOMPSAT-1 (3,298 events) and 38 days for Orion 50 XL (1,116 events). The lengthy unloading time is due to the debris’ high moment of inertia relative to the reaction wheels, causing frequent wheel saturation and requiring numerous unloading cycles.

Figure 10. Results of a Monte Carlo simulation for mean 3 bounds (red dashed lines) and errors in angular velocities and quaternion parameters (black solid lines) during a single circle.

Figure 10. Results of a Monte Carlo simulation for mean 3 bounds (red dashed lines) and errors in angular velocities and quaternion parameters (black solid lines) during a single circle.

When the spacecraft is in sunlight, sun sensor measurements enable the estimator to quickly align the predicted and estimated states, minimizing transient errors caused by aggressive initial conditions. However, during an eclipse, the estimator must rely on magnetometer and rate sensor data, which are less accurate than sun sensors. As a result, the accuracy of the estimated quaternion deteriorates.

To evaluate the robustness of the attitude estimator against various initial conditions, a Monte Carlo simulation with 200 runs was performed using a norm-constrained Extended Kalman Filter (EKF) for a nominal orbit. Plotting the angular velocity and quaternion errors from these runs. The quaternion errors, $q^T$, are computed using a multiplicative quaternion error formulation to maintain the unit length constraint.

The study focuses on two debris objects: the defunct KOMPSAT-1 satellite and the Orion 50 XL motor. As of June 1, 2017, KOMPSAT-1 was in a near-circular, highly inclined orbit with an average altitude of 650 km, while Orion 50 XL had a lower inclination but a more eccentric orbit, with a perigee altitude of about 420 km and an apogee altitude of 1400 km. Both orbits are within the operational region for the Deorbiter CubeSat mission.

Figure 11. The system’s low-thrust spiral orbital trajectory over the previous ninety days, along with the semi-major axes of debris being reduced to $a_f$ = 6478:1 km.

Figure 11. The system’s low-thrust spiral orbital trajectory over the previous ninety days, along with the semi-major axes of debris being reduced to $a_f$ = 6478:1 km.

During deorbiting maneuvers, the CubeSat propulsion system, aided by atmospheric drag, reduced KOMPSAT-1’s semi-major axis from 7,045 km to 6,478.1 km over 596 days (1.6 years). For Orion 50 XL, the semi-major axis decreased from 7,290 km to 6,478.1 km within 417 days. Without intervention, natural orbital decay for objects above 500 km altitude typically takes decades. Specifically, under optimal conditions, KOMPSAT-1 would remain in orbit for approximately 95 years and Orion 50 XL for about 22 years, assuming they maintain maximum cross-sectional areas perpendicular to their velocity.

The Deorbiter CubeSat significantly shortens the decay times of both debris objects. The simulation accounts for the J2 perturbation effect, which causes nodal regression as the debris undergoes deorbiting. While J2 causes secular changes in the right ascension of the ascending node (RAAN) and argument of perigee, it only produces periodic variations in the semi-major axis, inclination, and eccentricity. The nodal regression is more pronounced for Orion 50 XL due to its smaller orbital inclination compared to KOMPSAT-1.

8. Conclusions

This paper investigates the dynamics and control mechanisms of a novel adhesion detumbling strategy, specifically using a harpoon-nanosatellite system to stabilize space debris. A used rocket upper stage is selected as the target for this study. Initially, the collision dynamics between the detumbling device and the debris are modeled. Given the extremely brief duration of the collision, external disturbances are neglected. Consequently, the post-collision motion state of the combined system is determined by the principles of conservation of linear and angular momentum.

To precisely simulate the motion of the debris and nanosatellite system during the detumbling process, an orbit-attitude coupled dynamics model is constructed. This model incorporates time-varying external disturbances and variable inertial parameters. Since debris objects are typically non-cooperative and their inertial properties are challenging to determine accurately, an adaptive sliding mode control (SMC) scheme is developed. This control approach is robust to parameter uncertainties and incorporates input saturation to manage control efforts effectively. To convert the continuous control signal into a form suitable for thruster implementation, a Pulse Width Pulse Frequency (PWPF) modulator is employed.

Key findings from the simulations are summarized as follows:

• The collision significantly reduces the angular velocity of the debris from an initial state of ${\left[\begin{array}{ccc}28.6479& 22.9183& 34.3775\end{array}\right]}^T$ deg/s to ${\left[\begin{array}{ccc}20.3921& 21.4282& 9.2391\end{array}\right]}^T$ deg/s. This substantial change in angular velocity must be considered when planning detumbling operations.
• The proposed adaptive sliding mode control (SMC) scheme demonstrates high efficiency, low fuel consumption, and strong robustness. Under this control scheme, the angular velocity of the non-cooperative debris decreases rapidly within 4000 seconds (approximately one orbital period) and reaches near-zero within 8000 seconds.
• The total impulse generated by the thrusters increases sharply to 429.54 Ns at 3760 seconds and gradually rises to 438.35 Ns by 8000 seconds. The fuel consumed in the entire detumbling process is minimal, approximately 0.20 kg.
• The detumbling maneuvers do not significantly alter the target’s orbit, thereby minimizing the risk of potential collisions with other operational spacecraft in different orbits.

These results underscore the effectiveness of the harpoon-nanosatellite-based adhesion detumbling strategy for stabilizing non-cooperative space debris, providing a viable approach for future active debris removal missions.

9. Acknowledgements

I would like to extend their sincere gratitude to ClearSpace and the European Space Agency (ESA) for their invaluable assistance in verifying the simulations and providing essential computational resources. Their guidance and support have been pivotal in ensuring the accuracy and robustness of our results.

I also wish to thank the Indian Space Research Organisation (ISRO) for facilitating connections with these esteemed institutions, fostering an environment conducive to collaboration and knowledge exchange.

Furthermore, I gratefully acknowledge the support from the German Aerospace Center (DLR) and the Jugend Forscht program for providing grants that enabled me to pursue this research. Their funding and encouragement have been crucial in allowing us to explore innovative approaches to space debris removal.

Finally, I appreciate all the feedback and insights provided by various experts and collaborators during the research and manuscript preparation phase.