# Maneuver Strategy for Active Spacecraft to Avoid Space Debris and Return to the Original Orbit

Authors: Qun Fang; Zhen Zhang; Haodong Meng; Xiaolong Wang; Xiuwei Zhang
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

During normal operation of the on-orbit spacecraft, if some satellite in a nearby orbit suddenly breaks apart, its debris will threat the safe operation of the on-orbit spacecraft. Therefore, it is necessary to study the active spacecraft’s avoidance of the space debris group and returning to the original orbit. In this way, the safe operation of on-orbit spacecraft will be guaranteed. However, as the geometric structure of the space debris group is constantly changing, it is hard to accurately demonstrate the changing shape of the debris group, let alone determine the unreachable domain. Traditional obstacle avoidance problems involve low speed of the vehicle; so, the application of artificial potential field and particle swarm algorithms is suitable for such problems. However, these two methods are not applicable to the maneuver strategy of spacecraft with high initial velocity. Therefore, to help spacecraft avoid the space debris group, a new method is required. This paper has established a simplified model to simulate the unreachable domain of the space debris group. It has modified the artificial potential field (APF) method and particle swarm optimization algorithm, with an aim to help spacecraft avoid the space debris group and return to the original orbit. Based on the method, the paper has proposed a three-stage maneuver strategy for the spacecraft to avoid the debris. To show the effectiveness of the method, this paper has simulated an on-orbit spacecraft’s avoidance of the space debris group nearby and returning to its original orbit. Through simulation, the feasibility of the maneuver strategy for spacecraft in the geosynchronous orbit is evaluated. The simulation results show that the method proposed in this paper can effectively accomplish the task.


---

## Body

## 1. Introduction

With the rapid development of the space technology, the number of space orbit debris keeps increasing. According to the United States Space Surveillance Network, by the end of December 2021, there were more than 36,500 space debris with diameters larger than 10 cm littering space. More than 630 dangerous events have happened due to collisions or breakups of space debris [1]. As space debris fragments are travelling at fast speed, once they crash into the on-orbit spacecraft, the spacecraft will be severely damaged or even shattered. In history, many space collisions happened. For example, in 1991, Russian satellite COSMOS1934 collided with debris No. 13475 [2]; in 1996, the French spacecraft “Cerise” and space debris No. 18208 hit each other [3]; in 2005, THORBURNER 2A rocket collided with debris No. 26207 [3], and in 2009, two communications satellites—the commercial Iridium 33 and the Russian military Cosmos 2251—accidentally collided [4]. Every high-speed collision will scatter a larger number of smaller debris and form a vicious cycle. For example, the collision of Iridium 33 and Cosmos 2251 in 2009 generated 2,201 pieces of debris in total [5]. The collision between space objects not only poses a huge threat to the on-orbit spacecraft but also becomes the biggest contributor to the increasing of space debris [6, 7].With the development of space technology, the number of spacecrafts increases gradually. When the on-orbit spacecraft is working, if there is a sudden breakup of satellite near the orbit, a debris group will be produced. The debris group will severely threat the safe operation of on-orbit spacecraft. To reduce the threat caused by the space debris from satellite breakup, we need to study the maneuver strategy for the on-orbit spacecraft to avoid the space debris group and return to its original orbit. This is an important precondition for the safe operation of the spacecraft. Space debris and debris group pose different threats to spacecraft in orbit. The problem of active spacecraft evading space debris is a point-to-point problem. Only a fixed safety threshold needs to be considered when the spacecraft avoids space debris. In contrast, the space debris group spreads over time; so, it is necessary to consider the constantly changing geometry of the space debris group to avoid them. Currently, many scholars have studied the spacecraft avoidance strategy. Russell focused on calculating the collision probability among space objects, and he proposed to replace sphere colliders with cylindrical colliders [8]. After considering the general thrust and the control and maintenance of the spacecraft’s orbit and location, Chan proposed a maneuvering speed expression. His expression was based on the collision probability [9]. Alfano and Mueller proposed other appropriate strategies to solve the collision avoidance problem [10, 11]. Alfano analyzed the instantaneous maneuvering speed [10] while Mueller et al. used the standard form of optimal control to express the avoidance problem [11]. In Mueller’s study, the problem was discretized and transformed from a control problem into a nonlinear planning problem. Mueller et al. proposed a method to control the spacecraft and avoid on-orbit collision [12]. This method is applied to both situations when the two objects had a high relative velocity or a low one. It has a certain robustness. Besides, this method would prevent a second collision for a long time after the first maneuvering. Graziano et al. thoroughly introduced the basic principle for detecting and estimating the risks of satellite collision. He analyzed the operating procedure of collision prevention and summarized the rules for the maneuver strategy of avoidance [13]. Kelly and Piciotto adopted nonlinear optimization technology to study the optimal collision avoidance maneuvering calculation [14]. Based on the state-transition matrix, Gonzalo and Colombo designed an orbit maneuver strategy under the pulse thrust and analyzed its effectiveness [15]. They also analyzed the orbit maneuvering design under constant thrust provided along the velocity tangential direction [16] and added the results into a computer program to calculate the orbit maneuver [17]. Zhang et al. proposed an efficient method to calculate the diffusion law of space debris groups over time [18]. Dharmarajan et al. studied the method of satellite formation to avoid space debris. He proposed that the time of applying correction should be determined based on the capability of the satellite performing the maneuver. Moreover, he offered a method to solve eigenvalues in the optimization of maximum range and minimum velocity [19]. Most scholars mainly concentrated on how spacecraft can avoid space debris or regular space objects, but few of them studied how spacecraft can avoid the space debris group from the breakup of satellites in the neighboring orbit.This paper focuses on the maneuver strategy of spacecraft to avoid the space debris group, which are generated by the breakup of objects in neighboring orbit, and then return to its original orbit. It proposes a simplified method of the unreachable domain for the space debris group. Through the modification of artificial potential field (APF) and particle swarm optimization (PSO) algorithms, the spacecraft will complete the maneuvering process through three stages: the stage of avoidance, the stage of returning from a distant proximity, and the stage of returning from a close proximity. Next, the paper verifies the effectiveness of this method through a simulation experiment.


## 2. Dynamic Modeling

### 2.1. Selection of the Coordinate System

To accurately demonstrate the orbital motion of the active spacecraft and the space debris group, the earth-centered inertial (ECI) coordinate systemOXYZ is defined as can be shown in Figure 1: the coordinate system has its origin O at the center of the Earth, the OX axis is aligned with the mean equinox at epoch J2000, the OZ axis is aligned with the Earth’s North Pole, the OY axis is established by the right-hand rule, and P is an on-orbit spacecraft.Figure 1
Diagram of the Earth-Centered Inertial (ECI) coordinate system.To describe the orbital motion of an active spacecraft during the close return stage (Figure1), the CW coordinate system O′X′Y′Z′ is defined (Figure 2): the coordinate system has its origin O′ at the center of the mass of the target point, the O′X′ axis is aligned with the center of mass of the target along the center of the earth, the O′Y′ axis is aligned with the direction of motion of the target point, which is in the orbital plane of the target point and perpendicular to the O′X′ axis, and the O′Z′ axis is established by the right-hand rule.Figure 2
Diagram of the CW coordinate system.


### 2.2. The Dynamic Equation of the Active Spacecraft and the Space Debris Group

If the spacecraft can complete its task within a short period, it is reasonable to suppose that the influence of the perturbation term can be ignored. As the spacecraft needs to urgently avoid a close debris group, it is assumed that the spacecraft and space debris group are not affected by perturbation forces during their motion.The ECI system is selected to illustrate the orbital motion of the active spacecraft and space debris. The dynamic equation is(1)r¨=−μr3r,where r¨ is the acceleration vector of the space goal, μ denotes the geocentric gravitational constant, r means the geocentric distance of the space goal, and r represents the distance vector from the space goal to the center of Earth.


### 2.3. Discretization of the Relative Dynamic Equation

When the spacecraft is relatively close to the goal, the CW equation of the goal relative to the spacecraft is usually established under the CW coordinate system.(2)x¨−2nẏ−3n2x=fx,y¨+2nẋ=fy,z¨+n2z=fz,where x¨, y¨, and z¨ are the relative accelerations of the spacecraft to the goal in three directions of the CW coordinate system; ẋ and ẏ are the relative velocities of the spacecraft to the goal; x and z are the relative location of the spacecraft to the goal in the directions of O′X′ and O′Z′ under the CW coordinate system; fx, fy, and fz are the thrust accelerations of the spacecraft relative to the goal in the three directions of the CW coordinate system; and n is the mean angular velocity.
(3)n=μr03,where μ denotes the geocentric gravitational constant, and r0 is the goal’s geocentric distance.Based on linear systems theory, formula (2) is discretized [20], and the state-transition matrix of the spacecraft is obtained.
(4)Xt=Γt−t0Xt0+∫t0tΓt−τBuτdτ,where(5)Γt−t0=ΓΔt=ΓrrΔtΓrvΔtΓvrΔtΓvvΔt,Δt=t−t0,(6)ΓrrΔt=4−3cosnΔt006sinnt−nΔt1000cosnΔt,(7)ΓrvΔt=1nsinnΔt2n1−cosnΔt02ncosnΔt−11n4sinnΔt−3nΔt0001nsinnΔt,(8)ΓvrΔt=3nsinnΔt006ncosnΔt−10000−nsinnΔt,(9)ΓvvΔt=cosnΔt2sinnΔt0−2nsinnΔt4cosnΔt−3000cosnΔt.Considering that when the spacecraft returns to the goal, it will control the velocity impulse for many times, and this paper discretized the system’s continuous time intoT time intervals. The formula (4) can be shown as
(10)Xk+1=ΓkXk+Gkuk,where k indicates the k-th interval, uk=ukx,uky,ukzT denotes the accelerations generated by the spacecraft in three directions under the CW coordinate system, and Γk and Gk are the discretized state-space matrix. Their definitions are as follows:
(11)Γk=ΓT≜const,Gk=GT=∫0TΓτBdτ≜const.As the time intervalT is fixed, when the spacecraft is maneuvered for the N-th times, the state quantity of the spacecraft is
(12)XN=ΓNX0+Γ~u~,where
(13)Γ~=ΓN−1G,ΓN−2G,⋯,G,u~=u0,u1,⋯,uN−1T.Formula (12) is the dynamic model of the relative motion after the discretization of the average time. Using formula (12), we can obtain the motion state of the spacecraft relative to the goal at the time of t=t0+NT under the influence of the velocity impulse. Likewise, if the spacecraft wants to reach a certain state Xt+T at the time of t+T, formula (13) can also generate the amplitude and direction of the velocity impulse ut that is exerted by the spacecraft at the time of t.


## 2.1. Selection of the Coordinate System

To accurately demonstrate the orbital motion of the active spacecraft and the space debris group, the earth-centered inertial (ECI) coordinate systemOXYZ is defined as can be shown in Figure 1: the coordinate system has its origin O at the center of the Earth, the OX axis is aligned with the mean equinox at epoch J2000, the OZ axis is aligned with the Earth’s North Pole, the OY axis is established by the right-hand rule, and P is an on-orbit spacecraft.Figure 1
Diagram of the Earth-Centered Inertial (ECI) coordinate system.To describe the orbital motion of an active spacecraft during the close return stage (Figure1), the CW coordinate system O′X′Y′Z′ is defined (Figure 2): the coordinate system has its origin O′ at the center of the mass of the target point, the O′X′ axis is aligned with the center of mass of the target along the center of the earth, the O′Y′ axis is aligned with the direction of motion of the target point, which is in the orbital plane of the target point and perpendicular to the O′X′ axis, and the O′Z′ axis is established by the right-hand rule.Figure 2
Diagram of the CW coordinate system.


## 2.2. The Dynamic Equation of the Active Spacecraft and the Space Debris Group

If the spacecraft can complete its task within a short period, it is reasonable to suppose that the influence of the perturbation term can be ignored. As the spacecraft needs to urgently avoid a close debris group, it is assumed that the spacecraft and space debris group are not affected by perturbation forces during their motion.The ECI system is selected to illustrate the orbital motion of the active spacecraft and space debris. The dynamic equation is(1)r¨=−μr3r,where r¨ is the acceleration vector of the space goal, μ denotes the geocentric gravitational constant, r means the geocentric distance of the space goal, and r represents the distance vector from the space goal to the center of Earth.


## 2.3. Discretization of the Relative Dynamic Equation

When the spacecraft is relatively close to the goal, the CW equation of the goal relative to the spacecraft is usually established under the CW coordinate system.(2)x¨−2nẏ−3n2x=fx,y¨+2nẋ=fy,z¨+n2z=fz,where x¨, y¨, and z¨ are the relative accelerations of the spacecraft to the goal in three directions of the CW coordinate system; ẋ and ẏ are the relative velocities of the spacecraft to the goal; x and z are the relative location of the spacecraft to the goal in the directions of O′X′ and O′Z′ under the CW coordinate system; fx, fy, and fz are the thrust accelerations of the spacecraft relative to the goal in the three directions of the CW coordinate system; and n is the mean angular velocity.
(3)n=μr03,where μ denotes the geocentric gravitational constant, and r0 is the goal’s geocentric distance.Based on linear systems theory, formula (2) is discretized [20], and the state-transition matrix of the spacecraft is obtained.
(4)Xt=Γt−t0Xt0+∫t0tΓt−τBuτdτ,where(5)Γt−t0=ΓΔt=ΓrrΔtΓrvΔtΓvrΔtΓvvΔt,Δt=t−t0,(6)ΓrrΔt=4−3cosnΔt006sinnt−nΔt1000cosnΔt,(7)ΓrvΔt=1nsinnΔt2n1−cosnΔt02ncosnΔt−11n4sinnΔt−3nΔt0001nsinnΔt,(8)ΓvrΔt=3nsinnΔt006ncosnΔt−10000−nsinnΔt,(9)ΓvvΔt=cosnΔt2sinnΔt0−2nsinnΔt4cosnΔt−3000cosnΔt.Considering that when the spacecraft returns to the goal, it will control the velocity impulse for many times, and this paper discretized the system’s continuous time intoT time intervals. The formula (4) can be shown as
(10)Xk+1=ΓkXk+Gkuk,where k indicates the k-th interval, uk=ukx,uky,ukzT denotes the accelerations generated by the spacecraft in three directions under the CW coordinate system, and Γk and Gk are the discretized state-space matrix. Their definitions are as follows:
(11)Γk=ΓT≜const,Gk=GT=∫0TΓτBdτ≜const.As the time intervalT is fixed, when the spacecraft is maneuvered for the N-th times, the state quantity of the spacecraft is
(12)XN=ΓNX0+Γ~u~,where
(13)Γ~=ΓN−1G,ΓN−2G,⋯,G,u~=u0,u1,⋯,uN−1T.Formula (12) is the dynamic model of the relative motion after the discretization of the average time. Using formula (12), we can obtain the motion state of the spacecraft relative to the goal at the time of t=t0+NT under the influence of the velocity impulse. Likewise, if the spacecraft wants to reach a certain state Xt+T at the time of t+T, formula (13) can also generate the amplitude and direction of the velocity impulse ut that is exerted by the spacecraft at the time of t.


## 3. Generative Model of the Space Debris Group

To simulate the space debris group in the orbit near the active spacecraft, this paper adopts the model of NASA EVOLVE 4.0 [21] as an analog simulator to get information about the space debris fragments, including lengths, areas, mass, and speed increments.The breakup model of NASA EVOLVE 4.0 can be divided into two parts: the explosion breakup and the collision breakup. As they will not affect the application of the proposed algorithm of rapid orbital evolution, the two parts will be considered as one natural breakup, i.e., explosion breakup, without a loss of generality. Formula (14) calculates the size of every debris fragment:
(14)NLc≥Li=τ⋅6⋅Li−1.6,i=1,2⋯n,NLc=Li=NLc≥Li−NLc≥Li+1,where τ is the constant coefficient, 4.5 is the usually adopted value, Lc is the characteristic length, which is derived by taking the average number of the debris lengths in three directions, namely, Lc=lx+ly+lz/3, and NLc≥Li indicates the number of debris fragments with characteristic length longer than Li. The shorter the interval of Li is, the closer it is to the reality.Using the lengths of space debris, we can further determine the area-to-mass ratio (AMR) of each debris with binormal distribution (15):
(15)χ=αNμ1,σ1+1−αNμ2,σ2,Nμi,σi=12πσie−χ−μi22σi2,i=1,2,where α is the weight coefficient, and μi and σi are the mean and variance of the normal distribution, respectively. The characteristic lengths of the three parameters are obtained through the logarithm of the characteristic lengths log10Lc [19], and χ denotes the logarithm of the AMR with χ≙log10A/M.The mass of space debris is obtained by formula (16).
(16)M=Axη,where Ax is the average cross-sectional area. The average cross-sectional area of each piece of the space debris is determined by its characteristic length [19] and acquired through formula (17): η=10χ.
(17)Ax=0.540424Lc2,Lc<0.00167m,0.556945Lc2.0047077,Lc≥0.00167m.After the explosion and breakup, space debris will pick up a speed increment, which is obtained by formula (9):
(18)ϒ=12πσe−χ−μ22σ2,where ϒ is the logarithm of the speed increment, i.e., ϒ=log10Δv; μ is determined by χ, the logarithm of the AMR, i.e., μ=0.2χ+1.85; and σ is the constant of 0.4. The direction of the speed increment satisfies three dimensional random uniform distribution [22].


## 4. Maneuver Strategy of the Active Spacecraft’s Avoidance and Return

This section will establish a simplified model to simulate the unreachable domain of the space debris group. It will introduce the spacecraft maneuver strategy at the stage of avoidance, the stage of returning from a distant proximity, and the stage of returning from a proximity.


### 4.1. Modeling of the Unreachable Domain

For the space debris group, the state of every space debrisΩit can be simulated with a four-dimensional vector:
(19)Ωit=fxi,yi,zi,t⊂ℝ,where Ωit can be obtained through numerical integration of formula (1). Ωit represents the state of the i-th space debris at the time of t; xi, yi, and zi indicate the components of location vectors of the i-th space debris during its entire motion at the ECI coordinate system, respectively, and ℝ shows the state space consisting of all the space debris states during the entire motion.The orbital plane of the on-orbit spacecraft is defined as spaceℤ. Thus, the intersection of each space debris Ωi and the orbital plane is
(20)Qi=Ωi∩ℤ,whereQi=xi,yi,zi indicates the intersection of the i-th space debris and the spacecraft’s orbital plane, and xi, yi, and zi are the three coordinates of the intersection.The problem of active spacecraft’s avoidance of the space debris group is a four-dimensional problem; so, the active spacecraft needs to stay away from the space debrisΩi at the time of t. However, in reality, the calculation costs are high, as there are a large number of space debris in a debris group, which are hard to be analyzed one by one. Besides, the geometry of the space debris group will considerably change with the altering initial value. It is difficult to simulate using a unified mathematic rule. Therefore, dimension reduction is used in the model to present the threat of the space debris group. In other words, the two-dimensional space avoidance problem is transformed into an avoidance problem on a four-dimensional plane.According to formula (20), if Qi=Ωi∩ℤ=∅, there is no threat from the space debris to collide with the spacecraft; if Qi=Ωi∩ℤ≠∅, then Qi serves as the dangerous point for further discussion of modeling the unreachable domain.The set of coordinates for all dangerous points can be expressed as(21)ℤs_d=Q1,Q2,⋯,Qi,⋯,QN=ℝ∩ℤ.It should be noted thatℤs_d shows the set space of all the dangerous points. All the points inside the space are on the orbital plane ℤ:
(22)ℤs_d⊂ℤ..A circular envelope is used to encompass all the dangerous pointsℤs_d to form an unreachable domain ℤud. Thus, this unreachable domain is a circular no-fly zone on the orbital plane ℤ. The centre of the circle Qud is as follows:
(23)Qud=∑i=1NQiN.As for thei-th dangerous point Qi, which is the distance from Qi to the center of the unreachable circle, Qud can be derived as follows:
(24)ri=xi−xud2+yi−yud2+zi−zud2.Thus, the radius of the unreachable domain is(25)rud=maxr1,r2,⋯,ri,⋯,rN.Thus, the no-fly zone of the space debris group on the spacecraft’s orbital plane is obtained. If the active spacecraft can avoid the unreachable domain, of which the center of the circle isQud and the radius is rud, the danger from the space debris group can be avoided.


### 4.2. The Stages of the Spacecraft’s Task

Based on the requirement of the problem, the task of the spacecraft can be divided into the stage of avoidance, the stage of returning from a distant proximity, and the stage of returning from a close distance. To enable the spacecraft’s avoidance of the space debris group, an unreachable domain is established as the no-fly zone in Section3.1; to enable the spacecraft’s return to the original orbit, the goal is that the spacecraft should move back to the initial orbit. In other words, reaching back to the location and state of the spacecraft when it has not maneuvered is the goal.Suppose that the distance from the spacecraft to the no-fly zone at the time oftf is drepf, the distance at the time of tf+1 is drepf+1. If drepf>drepf+1, then the spacecraft should be considered as moving closer to the no-fly zone. At this time, the task of the spacecraft is to avoid the no-fly zone, and this period can be regarded as the stage of avoidance; conversely, if drepf≤drepf+1, the spacecraft is considered to be flying away from the no-fly zone, and this period can be regarded as the stage of returning from a distant proximity.Suppose that the distance from the spacecraft to the goal at the time oftk is dattk, and at the time of tk+1, its distance to the goal becomes dattk+1. If dattk>dattk+1, the spacecraft is considered to be closer to the goal. This is when the spacecraft is at the stage of returning from a distant proximity; conversely, if dattk≤dattk+1, the adopted method fails, and the spacecraft enters the stage of returning from a close proximity.


### 4.3. The Maneuver Strategy of the Spacecraft at the Stage of Avoidance

At the stage of avoidance, an artificial repulsive potential function is adopted to analyze the maneuver strategy for active spacecraft to avoid the space debris group. However, traditional repulsive potential functions are designed to solve problems such as the motion of robots with relatively small initial velocities or the rendezvous and docking of spacecrafts with a small relative velocity. For an active spacecraft to avoid the space debris group, there are some difficulties that have never occurred in traditional scenarios.The spacecraft itself is running at a fast speed and has a relatively high inertia, and the artificial repulsive potential field is only effective when the two objects are close to each other. Thus, if the spacecraft only receives notable impact from the repulsive potential field when it is close to the no-fly zone, it needs a thrust higher than its upper limit to overcome its inertia and avoid the no-fly zone.To solve these problems, the artificial potential field function needs to be modified. First, the distance values in the traditional repulsive potential are mapped to minimize the impact of the distance on the repulsive force; second, in contrast to the traditional repulsive potential field, the new method calculates the total force (thrust) instead of the repulsive potential field; third, the thrusts of the spacecraft at very moment are divided into the spacecraft velocity in the normal line direction [16].To help the spacecraft stay away the no-fly zone, an artificial repulsive force field is established. Its repulsive potential functionUrepq to the spacecraft is
(26)Urepq=25⋅krep⋅1qdrep2.5,drep≤d0,0,drep>d0,where krep is the gain coefficient, and drep indicates the Euclidean distance from the active spacecraft to the no-fly zone. The specific definition is shown in formula (27), and qdrep denotes the mapping of the distance drep within the interval of 1,2. The detailed definition is in formula (28), and d0 is the influence distance from the no-fly zone’s repulsive force field to the active spacecraft. If the distance between the spacecraft and the no-fly zone is smaller than this value, namely, drep≤d0, the spacecraft will be affected by the repulsive force. Conversely, the spacecraft will not be affected by the repulsive force.
(27)drep=x−xud2+y−yud2+z−zud2−rud,where x, y, and z are spacecraft’s location components under the three axes of the ECI system; xud, yud, and zud indicate the coordinates of the center of the no-fly zone under the three axes of the ECI system; and rud represents the radius of the no-fly zone.
(28)qdrep=drepd0+1.Due to the special condition of the spacecraft’s motion, when it is close to the no-fly zone, it is not able to maneuver and avoid the zone. Therefore, exerting distance from the repulsive force field needs to be projected to increase the influence of the field on the spacecraft when it enters the range of the repulsive force field.The negative gradient of the repulsive potential field is the repulsive force received by the active spacecraft, namely,(29)Frepq=−∇Urepq=−∂Urepq∂q=krep⋅1qdrep3.5,drep≤d00,drep>d0.Particularly, in formula (29), as qdrep⊂1,2, we get 1/qdrep3.5≤1; so, the gain coefficient of repulsive potential field krep can obtain the upper limit of the thrust from the spacecraft. In this way, the thrust of the spacecraft will be less than its upper limit.For further description, we define two concepts. At timet0, if the spacecraft is not maneuvering, the geocentric distance of spacecraft at the next moment t1 is Dn. If the spacecraft maneuvers, the geocentric distance of spacecraft at the next moment t1 is Dm. If Dm is longer than Dn, e the altered orbit of spacecraft at time t0 is defined as the ascending orbit. If Dm is shorter than Dn, the altered orbit of spacecraft at time t0 is defined as descending orbit.After the thrust is derived by formula (29), the direction of the thrust’s distribution is to be determined. As can be shown in Figure 3, when the center P of the unreachable domain is located at the original orbit (the initial orbit), where the spacecraft has not maneuvered to change orbit, the altered orbit of spacecraft should be the descending orbit (orbit 1); when the center P of the unreachable domain is inside the initial orbit (orbit 2), the altered orbit of spacecraft should be the ascending orbit. Next, the distribution direction of the thrust is to be determined based on whether the altered orbit of spacecraft is the ascending or descending orbit.Figure 3
Illustration of thrust distribution.Furthermore, the distribution of the thrust is illustrated in detail.Three nonlinear points on the orbital plane of the spacecraft are obtained. They arePsc1xsc1,ysc1,zsc1, Psc2xsc2,ysc2,zsc2, and Psc3xsc3,ysc3,zsc3, respectively. The normal vector of the spacecraft is
(30)nsc⟶=αn,βn,γn=Psc1Psc2⟶×Psc1Psc3⟶.The velocity vector of the active spacecraft at a certain time isvsc⟶=αv,βv,γv. It is assumed that the direction of the exerted thrust is Frep′⟶=αrep′,βrep′,γrep′. The thrust direction at this moment and the velocity vector vsc⟶ are perpendicular to the normal vector nsc⟶ of the orbital plane. Suppose that γrep′=1, then the dot product of vsc⟶ and nsc⟶ with Frep′⟶ is
(31)αrep′=βv⋅γv−αv⋅γnαv⋅αv⋅βn−βv−γvαv,βrep′=γv−αv⋅γnαv⋅βn−βv.Formula (31) is further processed by normalizing the thrust direction vector Frep′⟶=αrep′,βrep′,γrep′. The unit vector of the thrust direction Frep⟶=kFrep⋅αrep,βrep,γrep is
(32)αrep=αrep′αrep′2+βrep′2+γrep′2=αrep′αrep′2+βrep′2+1,βrep=βrep′αrep′2+βrep′2+γrep′2=βrep′αrep′2+βrep′2+1,γrep=γrep′αrep′2+βrep′2+γrep′2=1αrep′2+βrep′2+1.In Equation (32), kFrep is 1 or -1, representing two-unit vectors with different directions.As shown in Figure4, when the spacecraft moves to a random point S, the thrust calculated based on the repulsive potential function is either F1⟶ or F2⟶, determined by the value of kFrep.Figure 4
Diagram of the selection of the unit vector direction of the thrust direction.To further specify the thrust direction, takekv2>kv1>0, and the coordinates of A is kv1αv,kv1βv,kv1γv and B is kv2αv,kv2βv,kv2γv, while kFrep=1, the coordinates of point C is αrep,βrep,γrep.To figure out whether the pointC’s relative position to the vector AB on the orbital plane, project the points A, B, and C onto the xoy plane of the ECI coordinate system. Upon the projection, the coordinates of points A, B, and C are
(33)A′=kv1αv,kv1βv,B′=kv2αv,kv2βv,C′=αrep,βrep.ConnectingA′,B′,andC′ on the plane in the sequence to form an area which is defined as
(34)SA′,B′,C′=kv1αv−αrep·kv2βv−βrep−kv1βv−βrep·kv2αv−αrep.WhenA′,-B′,and -C′ are connected in a counterclockwise order, SA′,B′,C′ is positive. When A′,-B′,and -C′ are connected in a clockwise order, SA′,B′,C′ is negative. In other words, if the vector starts from A′ and ends at B′, then the C′ is on the right side of the vector A′B′ when Equation (33) is greater than 0, and the thrust in this direction will push the altered orbit of spacecraft ascend. If the C′ is on the left side of the vector A′B′ when Equation (34) is less than 0, the thrust in this direction will make the altered orbit of spacecraft descend.In summary, when the spacecraft needs to avoid collision by ascending it orbit, ifkFrep=1, the thrust makes the altered orbit of spacecraft ascend. Based on Equations (29)and (32), the thrust in the three directions of the spacecraft is
(35)Fx=αrep·Frepq,Fy=βrep·Frepq,Fz=γrep·Frepq,and vice versa.


### 4.4. Analysis on the Maneuver Strategy of Spacecraft at the Stage of Returning from a Distant Proximity

To design the orbit of the spacecraft returning from a distant proximity with artificial gravitational potential function, the center of the gravitational potential field will be set as the goal of the spacecraft when returning from a distant proximity. In order to keep a stable motion state for the spacecraft, the gravitational potential field center is designed to move along the initial orbit of the spacecraft. In other words, the position state of the nonmaneuvering spacecraft is considered as the center of the gravitational potential field.If the state vector at the initial moment of the spacecraft isx0=x0,y0,z0,vx0,vy0,vz0, the initial state vector of the gravitational potential field center will be xa0. By substituting the initial state vector x0 into Equation (1), the state vector of the center of the gravitational potential field at the time of tf will be xaf=xaf,yaf,zaf,vxaf,vyaf,vzaf.The gravitational potential functionUattq is
(36)Uattdatt=25·katt·datt2.5where katt is the gain coefficient of the gravitational potential field, datt is the Euclidean distance between the spacecraft and the center of the gravitational potential field, as defined in Equation (37).
(37)datt=xf′−xaf2+yf′−yaf2+zf′−zaf2,where xf′, yf′, and zf′ represent the location of the spacecraft at the moment of tf.Therefore, the gravitationFattdatt on the spacecraft is the negative gradient of the repulsive potential field function Uattdatt, denoted as
(38)Fattdatt=−∇Uattdatt=−katt·datt1.5.The result of Equation (38) represents the combined gravitational force on the spacecraft, whose components in the three directions in the ECI coordinate system are distributed according to the unit vector from the spacecraft to the center of the gravitational potential field.Accordingly, the vector from the spacecraft to the gravitational potential field center is the gravitational direction vectorFatt′⟶, which is denoted as
(39)Fatt′⟶=αatt′,βatt′,γatt′=xaf−xf′,yaf−yf′,zaf−zf′.After normalizing the gravitational direction vectorFatt′⟶, the unit vector along the gravitational direction is Fatt⟶=αatt,βatt,γatt, where
(40)αatt=αatt′αatt′2+βatt′2+γatt′2,βatt=βatt′αatt′2+βatt′2+γatt′2,γatt=γatt′αatt′2+βatt′2+γatt′2.According to Equations (38) and (40), the gravitational forces in the three directions of the spacecraft are
(41)Fx=αatt·Fattd,Fy=βatt·Fattd,Fz=γatt·Fattd.


### 4.5. Analysis on the Maneuvering Strategy of the Spacecraft at the Stage of Returning from a Close Proximity

Limited by the artificial gravitational potential function, the gravitational force of the spacecraft calculated when approaching the goal may be small, not accurate enough to guarantee the spacecraft return exactly to goal. Therefore, the study on the maneuvering strategy of spacecraft in the close proximity collision avoidance phase uses the PSO algorithm.For a better description, this section discusses the relative motion model under the CW coordinate system, while the selection of the goal remains the same as those in the previous section. Assuming that the spacecraft enters the stage of returning from a close proximity, the state quantity of the goal and the spacecraft in the ECI coordinate system isXA=rA,vA and XA=rB,vB, respectively. The CW coordinate system takes the coordinates of the target point rA as the origin; so, the relative state of the spacecraft relevant to the goal in the CW coordinate system is X0=r0,v0, where
(42)r0=rB−rA,v0=vB−vA−rA×vArA2×r0.To obtain the current desired outputU0 based on the current state quantity X0 and the state quantity X1 at the next moment, deform Equation (10) into
(43)U0=G−1X1−ΓX0Similarly, the spacecraft’s pulse maneuver time interval isT. For any state of the spacecraft at any t=t0+kT moment calculated through Equation (12), the output set at each time interval is U~,; so,
(44)U~=Γ~−1Xk−ΓkX0,where u~=u0,u1,⋯,uN−1T.Assuming that the spacecraft needs a total ofNpulses to complete the mission, to reach the goal step-by-step, an optimization model for the spacecraft reaching the goal will be built as
(45)s.t.0≤uk≤umax.In the model,DS,Pk denotes the distance to the goal at the k-th pulse of the spacecraft, uk denotes the 2-norm of the pulse vector at the k-th maneuver of the spacecraft, umax denotes the maximum input control that can be generated by the spacecraft, and α1 and α2 are weighting factors.To study the maneuver strategy of the spacecraft in respect of returning to the goal, the particles in the PSO algorithm are given actual physical meanings. Each particle position represents a possible location where the spacecraft might reach. Based on that, the solution of the particle after PSO optimization will be a node along the maneuvering path of the spacecraft. Through Equation (45), the optimal value searched by each particle in the particle swarm exchanges information, and the optimal value produced in every search will be the end position of the motion path for spacecraft at this moment. The search process will not stop repeating until it reaches the goal. The Equation (43) calculates the pulse velocity increment at each moment while the spacecraft is moving to the goal.To update the velocity and position of each particle in the PSO algorithm, the equation is(46)vi,dk+1=ωvi,dk+c1r1pi,dk−xi,dk+c2r2pg,dk−xi,dk,xi,dk+1=xi,dk+vi,dk+1,where the physical meaning of symbols in the equation is as shown in Table 1.Table 1
Physical meaning of symbols in Equation (46).
SymbolMeaningiParticle numberkThe times of current iterationsdThe dimension of the solutionvi,dkVelocity information of thei-th particle on the d-th dimensional solution at the k-th iterationxi,dkThe position of thei-th particle on the d dimensional solution at the k-th iterationpi,dkThe individual optimal position of thei-th particle for the previous k iterationspg,dkThe global optimal positions of all particles for the previousk iterationsωInertia weightsc1,c2Nonnegative constantsr1,r2Any constants between [0,1]The inertia weightω is important for particle convergence and mitigating the contradiction between particle swarm search and exploitation. A larger ω is helpful to increase the diversity of the swarm, while a smaller ω is beneficial for the exploitation of local optima. When planning the return path of the spacecraft to the goal, the ω requires specifically settings: at the initial stage, a larger search range is needed for spacecraft to quickly approach the goal; however, when the spacecraft is approximate to the goal, local optimum solution can be produced by searching in a smaller range, thus reducing the repeated oscillations. To enhance the search rate and accuracy, the study adopts the inertia weights dynamically changing:
(47)ω=ωkmaxkmaxωminmaxmin,where ωmin represents the minimum inertia weight, ωmax is the maximum inertia weight, k is the number of current iterations, and kmax is the maximum times of iterations.Supposef· denotes the fitness function and npnp is the number of particles in the particle swarm, to find the minimal value, the individual optimal position pi,dkpi,dk for the first k iterations of the i-th particle is calculated in the equation:
(48)pi,dk+1=pi,dk,iffxi,dk+1≥fpi,dk,xi,dk+1,iffxi,dk+1<fpi,dk.For the global optimal positionpg,dk in the first k iterations of all particles, the equation is
(49)pg,dk∈p1,dk,p2,dk,⋯,pnp,dkfpg,dk=minfp1,dk,fp2,dk,⋯,fpnp,dk.When the velocity of the particle swarm is updating, it might be very large, even causing the particles to rush out of the solution range or even diverge. In order to control the search action of each particle, the particle flight speed has to be controlled within a specific range. Therefore, the particle velocity equation will be updated as(50)vi,dk+1=vi,dk+1,ifvi,dk+1≤Vmax,d,Vmax,d,ifvi,dk+1>Vmax,d,−Vmax,d,ifvi,dk+1<−Vmax,d,where Vmax,d is the maximum velocity that a particle can have in the d-th dimension.Since the search distance of particle swarm is usually restricted to certain constraints, for example, at the initial search moment of the PSO algorithm, it is expected that a larger search radius can motivate the random walk of particles. As the spacecraft approaches the goal, a larger search radius may lead to a large computational cost caused by the PSO algorithm when converging to the goal.As a result, this paper constrains the initial values of each particle in the swarm to a sphere region with radius ofRPSO. The position of the sphere center is determined by the position information of the spacecraft at the moment, while the radius of the sphere is depending on the distance between the spacecraft and the corresponding goal. In this paper, the radius is proportional to the distance between the spacecraft and the goal. The radius of the particle swarm algorithm search sphere is RPSO; so,
(51)RPSO=β11+e−β2·DS,Pk+β3+β4,where β1，β2，β3，β4 are coefficients.It should be noted that the particles in the PSO algorithm might move out of the search radius of the particle swarm when performing random walk; so, boundary treatment is also needed.As shown in Figure5, if the spacecraft is located at A, at the time t, the sphere represents the search range of PSO algorithm. If a location B is searched by particle i in the swarm, it indicates that the particle has reached beyond the search boundary of the particle swarm, requiring boundary treatment. In order to reduce the influence of the treatment on the iterations and evolutionary direction of the particle swarm, this paper contracts the position of particle i along the vector AB⟶ and corrects the location of particle i from point B to point B′.Figure 5
Diagram of modifying particles that are beyond the search range.Suppose the coordinates of points arex1,y1,z1, x2,y2,z2, and B′=a,b,cx2,y2,z2, then the unit vector η=ηx,ηy,ηz of vector AB⟶ is
(52)ηx=x2−x1x2−x12+y2−y12+z2−z12,ηy=y2−y1x2−x12+y2−y12+z2−z12,ηz=z2−z1x2−x12+y2−y12+z2−z12.AsOB′⟶=OA⟶+AB′⟶,, the modified coordinates of the particles are
(53)a=x1+ηx·RPSO,b=x2+ηy·RPSO,c=x3+ηz·RPSO.Figure6 shows the PSO search path used by the spacecraft. In conclusion, at the moment of t0t0, the spacecraft is at the initial position, which is taken as the sphere center. In this sphere with the radius of R0, np particles are randomly generated. The global optimal particle position will be selected as the spacecraft position at t0+T after kmax iterations. Based on this and repeat the computation, we can plan the tracking and return path of the spacecraft to the goal.Figure 6
Diagram of the Spacecraft Using PSO Algorithm to Search the Path.


## 4.1. Modeling of the Unreachable Domain

For the space debris group, the state of every space debrisΩit can be simulated with a four-dimensional vector:
(19)Ωit=fxi,yi,zi,t⊂ℝ,where Ωit can be obtained through numerical integration of formula (1). Ωit represents the state of the i-th space debris at the time of t; xi, yi, and zi indicate the components of location vectors of the i-th space debris during its entire motion at the ECI coordinate system, respectively, and ℝ shows the state space consisting of all the space debris states during the entire motion.The orbital plane of the on-orbit spacecraft is defined as spaceℤ. Thus, the intersection of each space debris Ωi and the orbital plane is
(20)Qi=Ωi∩ℤ,whereQi=xi,yi,zi indicates the intersection of the i-th space debris and the spacecraft’s orbital plane, and xi, yi, and zi are the three coordinates of the intersection.The problem of active spacecraft’s avoidance of the space debris group is a four-dimensional problem; so, the active spacecraft needs to stay away from the space debrisΩi at the time of t. However, in reality, the calculation costs are high, as there are a large number of space debris in a debris group, which are hard to be analyzed one by one. Besides, the geometry of the space debris group will considerably change with the altering initial value. It is difficult to simulate using a unified mathematic rule. Therefore, dimension reduction is used in the model to present the threat of the space debris group. In other words, the two-dimensional space avoidance problem is transformed into an avoidance problem on a four-dimensional plane.According to formula (20), if Qi=Ωi∩ℤ=∅, there is no threat from the space debris to collide with the spacecraft; if Qi=Ωi∩ℤ≠∅, then Qi serves as the dangerous point for further discussion of modeling the unreachable domain.The set of coordinates for all dangerous points can be expressed as(21)ℤs_d=Q1,Q2,⋯,Qi,⋯,QN=ℝ∩ℤ.It should be noted thatℤs_d shows the set space of all the dangerous points. All the points inside the space are on the orbital plane ℤ:
(22)ℤs_d⊂ℤ..A circular envelope is used to encompass all the dangerous pointsℤs_d to form an unreachable domain ℤud. Thus, this unreachable domain is a circular no-fly zone on the orbital plane ℤ. The centre of the circle Qud is as follows:
(23)Qud=∑i=1NQiN.As for thei-th dangerous point Qi, which is the distance from Qi to the center of the unreachable circle, Qud can be derived as follows:
(24)ri=xi−xud2+yi−yud2+zi−zud2.Thus, the radius of the unreachable domain is(25)rud=maxr1,r2,⋯,ri,⋯,rN.Thus, the no-fly zone of the space debris group on the spacecraft’s orbital plane is obtained. If the active spacecraft can avoid the unreachable domain, of which the center of the circle isQud and the radius is rud, the danger from the space debris group can be avoided.


## 4.2. The Stages of the Spacecraft’s Task

Based on the requirement of the problem, the task of the spacecraft can be divided into the stage of avoidance, the stage of returning from a distant proximity, and the stage of returning from a close distance. To enable the spacecraft’s avoidance of the space debris group, an unreachable domain is established as the no-fly zone in Section3.1; to enable the spacecraft’s return to the original orbit, the goal is that the spacecraft should move back to the initial orbit. In other words, reaching back to the location and state of the spacecraft when it has not maneuvered is the goal.Suppose that the distance from the spacecraft to the no-fly zone at the time oftf is drepf, the distance at the time of tf+1 is drepf+1. If drepf>drepf+1, then the spacecraft should be considered as moving closer to the no-fly zone. At this time, the task of the spacecraft is to avoid the no-fly zone, and this period can be regarded as the stage of avoidance; conversely, if drepf≤drepf+1, the spacecraft is considered to be flying away from the no-fly zone, and this period can be regarded as the stage of returning from a distant proximity.Suppose that the distance from the spacecraft to the goal at the time oftk is dattk, and at the time of tk+1, its distance to the goal becomes dattk+1. If dattk>dattk+1, the spacecraft is considered to be closer to the goal. This is when the spacecraft is at the stage of returning from a distant proximity; conversely, if dattk≤dattk+1, the adopted method fails, and the spacecraft enters the stage of returning from a close proximity.


## 4.3. The Maneuver Strategy of the Spacecraft at the Stage of Avoidance

At the stage of avoidance, an artificial repulsive potential function is adopted to analyze the maneuver strategy for active spacecraft to avoid the space debris group. However, traditional repulsive potential functions are designed to solve problems such as the motion of robots with relatively small initial velocities or the rendezvous and docking of spacecrafts with a small relative velocity. For an active spacecraft to avoid the space debris group, there are some difficulties that have never occurred in traditional scenarios.The spacecraft itself is running at a fast speed and has a relatively high inertia, and the artificial repulsive potential field is only effective when the two objects are close to each other. Thus, if the spacecraft only receives notable impact from the repulsive potential field when it is close to the no-fly zone, it needs a thrust higher than its upper limit to overcome its inertia and avoid the no-fly zone.To solve these problems, the artificial potential field function needs to be modified. First, the distance values in the traditional repulsive potential are mapped to minimize the impact of the distance on the repulsive force; second, in contrast to the traditional repulsive potential field, the new method calculates the total force (thrust) instead of the repulsive potential field; third, the thrusts of the spacecraft at very moment are divided into the spacecraft velocity in the normal line direction [16].To help the spacecraft stay away the no-fly zone, an artificial repulsive force field is established. Its repulsive potential functionUrepq to the spacecraft is
(26)Urepq=25⋅krep⋅1qdrep2.5,drep≤d0,0,drep>d0,where krep is the gain coefficient, and drep indicates the Euclidean distance from the active spacecraft to the no-fly zone. The specific definition is shown in formula (27), and qdrep denotes the mapping of the distance drep within the interval of 1,2. The detailed definition is in formula (28), and d0 is the influence distance from the no-fly zone’s repulsive force field to the active spacecraft. If the distance between the spacecraft and the no-fly zone is smaller than this value, namely, drep≤d0, the spacecraft will be affected by the repulsive force. Conversely, the spacecraft will not be affected by the repulsive force.
(27)drep=x−xud2+y−yud2+z−zud2−rud,where x, y, and z are spacecraft’s location components under the three axes of the ECI system; xud, yud, and zud indicate the coordinates of the center of the no-fly zone under the three axes of the ECI system; and rud represents the radius of the no-fly zone.
(28)qdrep=drepd0+1.Due to the special condition of the spacecraft’s motion, when it is close to the no-fly zone, it is not able to maneuver and avoid the zone. Therefore, exerting distance from the repulsive force field needs to be projected to increase the influence of the field on the spacecraft when it enters the range of the repulsive force field.The negative gradient of the repulsive potential field is the repulsive force received by the active spacecraft, namely,(29)Frepq=−∇Urepq=−∂Urepq∂q=krep⋅1qdrep3.5,drep≤d00,drep>d0.Particularly, in formula (29), as qdrep⊂1,2, we get 1/qdrep3.5≤1; so, the gain coefficient of repulsive potential field krep can obtain the upper limit of the thrust from the spacecraft. In this way, the thrust of the spacecraft will be less than its upper limit.For further description, we define two concepts. At timet0, if the spacecraft is not maneuvering, the geocentric distance of spacecraft at the next moment t1 is Dn. If the spacecraft maneuvers, the geocentric distance of spacecraft at the next moment t1 is Dm. If Dm is longer than Dn, e the altered orbit of spacecraft at time t0 is defined as the ascending orbit. If Dm is shorter than Dn, the altered orbit of spacecraft at time t0 is defined as descending orbit.After the thrust is derived by formula (29), the direction of the thrust’s distribution is to be determined. As can be shown in Figure 3, when the center P of the unreachable domain is located at the original orbit (the initial orbit), where the spacecraft has not maneuvered to change orbit, the altered orbit of spacecraft should be the descending orbit (orbit 1); when the center P of the unreachable domain is inside the initial orbit (orbit 2), the altered orbit of spacecraft should be the ascending orbit. Next, the distribution direction of the thrust is to be determined based on whether the altered orbit of spacecraft is the ascending or descending orbit.Figure 3
Illustration of thrust distribution.Furthermore, the distribution of the thrust is illustrated in detail.Three nonlinear points on the orbital plane of the spacecraft are obtained. They arePsc1xsc1,ysc1,zsc1, Psc2xsc2,ysc2,zsc2, and Psc3xsc3,ysc3,zsc3, respectively. The normal vector of the spacecraft is
(30)nsc⟶=αn,βn,γn=Psc1Psc2⟶×Psc1Psc3⟶.The velocity vector of the active spacecraft at a certain time isvsc⟶=αv,βv,γv. It is assumed that the direction of the exerted thrust is Frep′⟶=αrep′,βrep′,γrep′. The thrust direction at this moment and the velocity vector vsc⟶ are perpendicular to the normal vector nsc⟶ of the orbital plane. Suppose that γrep′=1, then the dot product of vsc⟶ and nsc⟶ with Frep′⟶ is
(31)αrep′=βv⋅γv−αv⋅γnαv⋅αv⋅βn−βv−γvαv,βrep′=γv−αv⋅γnαv⋅βn−βv.Formula (31) is further processed by normalizing the thrust direction vector Frep′⟶=αrep′,βrep′,γrep′. The unit vector of the thrust direction Frep⟶=kFrep⋅αrep,βrep,γrep is
(32)αrep=αrep′αrep′2+βrep′2+γrep′2=αrep′αrep′2+βrep′2+1,βrep=βrep′αrep′2+βrep′2+γrep′2=βrep′αrep′2+βrep′2+1,γrep=γrep′αrep′2+βrep′2+γrep′2=1αrep′2+βrep′2+1.In Equation (32), kFrep is 1 or -1, representing two-unit vectors with different directions.As shown in Figure4, when the spacecraft moves to a random point S, the thrust calculated based on the repulsive potential function is either F1⟶ or F2⟶, determined by the value of kFrep.Figure 4
Diagram of the selection of the unit vector direction of the thrust direction.To further specify the thrust direction, takekv2>kv1>0, and the coordinates of A is kv1αv,kv1βv,kv1γv and B is kv2αv,kv2βv,kv2γv, while kFrep=1, the coordinates of point C is αrep,βrep,γrep.To figure out whether the pointC’s relative position to the vector AB on the orbital plane, project the points A, B, and C onto the xoy plane of the ECI coordinate system. Upon the projection, the coordinates of points A, B, and C are
(33)A′=kv1αv,kv1βv,B′=kv2αv,kv2βv,C′=αrep,βrep.ConnectingA′,B′,andC′ on the plane in the sequence to form an area which is defined as
(34)SA′,B′,C′=kv1αv−αrep·kv2βv−βrep−kv1βv−βrep·kv2αv−αrep.WhenA′,-B′,and -C′ are connected in a counterclockwise order, SA′,B′,C′ is positive. When A′,-B′,and -C′ are connected in a clockwise order, SA′,B′,C′ is negative. In other words, if the vector starts from A′ and ends at B′, then the C′ is on the right side of the vector A′B′ when Equation (33) is greater than 0, and the thrust in this direction will push the altered orbit of spacecraft ascend. If the C′ is on the left side of the vector A′B′ when Equation (34) is less than 0, the thrust in this direction will make the altered orbit of spacecraft descend.In summary, when the spacecraft needs to avoid collision by ascending it orbit, ifkFrep=1, the thrust makes the altered orbit of spacecraft ascend. Based on Equations (29)and (32), the thrust in the three directions of the spacecraft is
(35)Fx=αrep·Frepq,Fy=βrep·Frepq,Fz=γrep·Frepq,and vice versa.


## 4.4. Analysis on the Maneuver Strategy of Spacecraft at the Stage of Returning from a Distant Proximity

To design the orbit of the spacecraft returning from a distant proximity with artificial gravitational potential function, the center of the gravitational potential field will be set as the goal of the spacecraft when returning from a distant proximity. In order to keep a stable motion state for the spacecraft, the gravitational potential field center is designed to move along the initial orbit of the spacecraft. In other words, the position state of the nonmaneuvering spacecraft is considered as the center of the gravitational potential field.If the state vector at the initial moment of the spacecraft isx0=x0,y0,z0,vx0,vy0,vz0, the initial state vector of the gravitational potential field center will be xa0. By substituting the initial state vector x0 into Equation (1), the state vector of the center of the gravitational potential field at the time of tf will be xaf=xaf,yaf,zaf,vxaf,vyaf,vzaf.The gravitational potential functionUattq is
(36)Uattdatt=25·katt·datt2.5where katt is the gain coefficient of the gravitational potential field, datt is the Euclidean distance between the spacecraft and the center of the gravitational potential field, as defined in Equation (37).
(37)datt=xf′−xaf2+yf′−yaf2+zf′−zaf2,where xf′, yf′, and zf′ represent the location of the spacecraft at the moment of tf.Therefore, the gravitationFattdatt on the spacecraft is the negative gradient of the repulsive potential field function Uattdatt, denoted as
(38)Fattdatt=−∇Uattdatt=−katt·datt1.5.The result of Equation (38) represents the combined gravitational force on the spacecraft, whose components in the three directions in the ECI coordinate system are distributed according to the unit vector from the spacecraft to the center of the gravitational potential field.Accordingly, the vector from the spacecraft to the gravitational potential field center is the gravitational direction vectorFatt′⟶, which is denoted as
(39)Fatt′⟶=αatt′,βatt′,γatt′=xaf−xf′,yaf−yf′,zaf−zf′.After normalizing the gravitational direction vectorFatt′⟶, the unit vector along the gravitational direction is Fatt⟶=αatt,βatt,γatt, where
(40)αatt=αatt′αatt′2+βatt′2+γatt′2,βatt=βatt′αatt′2+βatt′2+γatt′2,γatt=γatt′αatt′2+βatt′2+γatt′2.According to Equations (38) and (40), the gravitational forces in the three directions of the spacecraft are
(41)Fx=αatt·Fattd,Fy=βatt·Fattd,Fz=γatt·Fattd.


## 4.5. Analysis on the Maneuvering Strategy of the Spacecraft at the Stage of Returning from a Close Proximity

Limited by the artificial gravitational potential function, the gravitational force of the spacecraft calculated when approaching the goal may be small, not accurate enough to guarantee the spacecraft return exactly to goal. Therefore, the study on the maneuvering strategy of spacecraft in the close proximity collision avoidance phase uses the PSO algorithm.For a better description, this section discusses the relative motion model under the CW coordinate system, while the selection of the goal remains the same as those in the previous section. Assuming that the spacecraft enters the stage of returning from a close proximity, the state quantity of the goal and the spacecraft in the ECI coordinate system isXA=rA,vA and XA=rB,vB, respectively. The CW coordinate system takes the coordinates of the target point rA as the origin; so, the relative state of the spacecraft relevant to the goal in the CW coordinate system is X0=r0,v0, where
(42)r0=rB−rA,v0=vB−vA−rA×vArA2×r0.To obtain the current desired outputU0 based on the current state quantity X0 and the state quantity X1 at the next moment, deform Equation (10) into
(43)U0=G−1X1−ΓX0Similarly, the spacecraft’s pulse maneuver time interval isT. For any state of the spacecraft at any t=t0+kT moment calculated through Equation (12), the output set at each time interval is U~,; so,
(44)U~=Γ~−1Xk−ΓkX0,where u~=u0,u1,⋯,uN−1T.Assuming that the spacecraft needs a total ofNpulses to complete the mission, to reach the goal step-by-step, an optimization model for the spacecraft reaching the goal will be built as
(45)s.t.0≤uk≤umax.In the model,DS,Pk denotes the distance to the goal at the k-th pulse of the spacecraft, uk denotes the 2-norm of the pulse vector at the k-th maneuver of the spacecraft, umax denotes the maximum input control that can be generated by the spacecraft, and α1 and α2 are weighting factors.To study the maneuver strategy of the spacecraft in respect of returning to the goal, the particles in the PSO algorithm are given actual physical meanings. Each particle position represents a possible location where the spacecraft might reach. Based on that, the solution of the particle after PSO optimization will be a node along the maneuvering path of the spacecraft. Through Equation (45), the optimal value searched by each particle in the particle swarm exchanges information, and the optimal value produced in every search will be the end position of the motion path for spacecraft at this moment. The search process will not stop repeating until it reaches the goal. The Equation (43) calculates the pulse velocity increment at each moment while the spacecraft is moving to the goal.To update the velocity and position of each particle in the PSO algorithm, the equation is(46)vi,dk+1=ωvi,dk+c1r1pi,dk−xi,dk+c2r2pg,dk−xi,dk,xi,dk+1=xi,dk+vi,dk+1,where the physical meaning of symbols in the equation is as shown in Table 1.Table 1
Physical meaning of symbols in Equation (46).
SymbolMeaningiParticle numberkThe times of current iterationsdThe dimension of the solutionvi,dkVelocity information of thei-th particle on the d-th dimensional solution at the k-th iterationxi,dkThe position of thei-th particle on the d dimensional solution at the k-th iterationpi,dkThe individual optimal position of thei-th particle for the previous k iterationspg,dkThe global optimal positions of all particles for the previousk iterationsωInertia weightsc1,c2Nonnegative constantsr1,r2Any constants between [0,1]The inertia weightω is important for particle convergence and mitigating the contradiction between particle swarm search and exploitation. A larger ω is helpful to increase the diversity of the swarm, while a smaller ω is beneficial for the exploitation of local optima. When planning the return path of the spacecraft to the goal, the ω requires specifically settings: at the initial stage, a larger search range is needed for spacecraft to quickly approach the goal; however, when the spacecraft is approximate to the goal, local optimum solution can be produced by searching in a smaller range, thus reducing the repeated oscillations. To enhance the search rate and accuracy, the study adopts the inertia weights dynamically changing:
(47)ω=ωkmaxkmaxωminmaxmin,where ωmin represents the minimum inertia weight, ωmax is the maximum inertia weight, k is the number of current iterations, and kmax is the maximum times of iterations.Supposef· denotes the fitness function and npnp is the number of particles in the particle swarm, to find the minimal value, the individual optimal position pi,dkpi,dk for the first k iterations of the i-th particle is calculated in the equation:
(48)pi,dk+1=pi,dk,iffxi,dk+1≥fpi,dk,xi,dk+1,iffxi,dk+1<fpi,dk.For the global optimal positionpg,dk in the first k iterations of all particles, the equation is
(49)pg,dk∈p1,dk,p2,dk,⋯,pnp,dkfpg,dk=minfp1,dk,fp2,dk,⋯,fpnp,dk.When the velocity of the particle swarm is updating, it might be very large, even causing the particles to rush out of the solution range or even diverge. In order to control the search action of each particle, the particle flight speed has to be controlled within a specific range. Therefore, the particle velocity equation will be updated as(50)vi,dk+1=vi,dk+1,ifvi,dk+1≤Vmax,d,Vmax,d,ifvi,dk+1>Vmax,d,−Vmax,d,ifvi,dk+1<−Vmax,d,where Vmax,d is the maximum velocity that a particle can have in the d-th dimension.Since the search distance of particle swarm is usually restricted to certain constraints, for example, at the initial search moment of the PSO algorithm, it is expected that a larger search radius can motivate the random walk of particles. As the spacecraft approaches the goal, a larger search radius may lead to a large computational cost caused by the PSO algorithm when converging to the goal.As a result, this paper constrains the initial values of each particle in the swarm to a sphere region with radius ofRPSO. The position of the sphere center is determined by the position information of the spacecraft at the moment, while the radius of the sphere is depending on the distance between the spacecraft and the corresponding goal. In this paper, the radius is proportional to the distance between the spacecraft and the goal. The radius of the particle swarm algorithm search sphere is RPSO; so,
(51)RPSO=β11+e−β2·DS,Pk+β3+β4,where β1，β2，β3，β4 are coefficients.It should be noted that the particles in the PSO algorithm might move out of the search radius of the particle swarm when performing random walk; so, boundary treatment is also needed.As shown in Figure5, if the spacecraft is located at A, at the time t, the sphere represents the search range of PSO algorithm. If a location B is searched by particle i in the swarm, it indicates that the particle has reached beyond the search boundary of the particle swarm, requiring boundary treatment. In order to reduce the influence of the treatment on the iterations and evolutionary direction of the particle swarm, this paper contracts the position of particle i along the vector AB⟶ and corrects the location of particle i from point B to point B′.Figure 5
Diagram of modifying particles that are beyond the search range.Suppose the coordinates of points arex1,y1,z1, x2,y2,z2, and B′=a,b,cx2,y2,z2, then the unit vector η=ηx,ηy,ηz of vector AB⟶ is
(52)ηx=x2−x1x2−x12+y2−y12+z2−z12,ηy=y2−y1x2−x12+y2−y12+z2−z12,ηz=z2−z1x2−x12+y2−y12+z2−z12.AsOB′⟶=OA⟶+AB′⟶,, the modified coordinates of the particles are
(53)a=x1+ηx·RPSO,b=x2+ηy·RPSO,c=x3+ηz·RPSO.Figure6 shows the PSO search path used by the spacecraft. In conclusion, at the moment of t0t0, the spacecraft is at the initial position, which is taken as the sphere center. In this sphere with the radius of R0, np particles are randomly generated. The global optimal particle position will be selected as the spacecraft position at t0+T after kmax iterations. Based on this and repeat the computation, we can plan the tracking and return path of the spacecraft to the goal.Figure 6
Diagram of the Spacecraft Using PSO Algorithm to Search the Path.


## 5. Simulation Analysis

The analysis first used the NASA EVOLVE 4.0 model to simulate the generation of space debris group. Table2 shows the six orbital elements of the active spacecraft at the initial stage and of the spacecraft that are disintegrating.Table 2
The six elements of initial moment active spacecraft and disintegrating spacecraft disintegrating moment orbit.
Active spacecraftDisintegrating spacecraftSemimajor axis/km4216442014Eccentricity00Inclination/°1560Right ascension of the ascending node/°00Argument of perigee/°38245True anomaly/°316142The disintegrating spacecraft spawns space debris group disintegration point, creating a total of 15,573 pieces of debris, while debris less than 0.5 g are not counted in the analysis. The circle center of the unreachable domain is (42014, 0, 0) km with a radius of 158 km.The mass of the active spacecraft is 1000 kg, the spacecraft thrust is no more than 100 N, the repulsive potential field gain coefficientkrep is 100, taken from the upper limit of spacecraft thrust, and the repulsive potential field action range is d0=3000km.At the stage of avoidance and the stage of returning from the distant proximity, the thrust magnitude in the three directions and total thrust of the spacecraft are as shown in Figure7. It should be noted that the thrust and the acceleration are in different directions, indicating that the spacecraft is under the repulsive potential field from 173 to 1200 seconds and under the gravitational potential field from 1200 to 2446 seconds.Figure 7
The thrust of the active spacecraft at the stage of avoidance and the stage of returning from the distant proximity: (a) thrust in theX direction, (b) thrust in the Y direction, (c) thrust in the Z direction, and (d) total thrust (scalar).
(a)(b)(c)(d)In Figure7, the spacecraft is less than 3000 km away from the unreachable domain at about 173 seconds when entering the stage of avoidance. At this timestamp, the artificial repulsion function is used to make the spacecraft gradually move away from the no-fly zone, during which the spacecraft enters the stage of approaching from the distant proximity from 1200 seconds to 2446 seconds. By then, the artificial gravitational function enables the spacecraft to obtain a gravity pointing at the target point, thus gradually approaching the goal. Limited by the artificial potential function, the distance between the spacecraft and the goal is the shortest at 2446 seconds when the gravity provided by the gravitational potential field is insufficient to pull the spacecraft back to the goal in such a short distance. After that, the spacecraft enters the stage of returning from the stage of close proximity.It is worth mentioning that when the active spacecraft enters the stage of returning from the distant proximity, the spacecraft’s initial velocity is too big for the gravity calculated with the gravitational potential field function to change the spacecraft’s motion trend in a short period of time. As a result, the distance between the spacecraft and the target point will gradually extend at the beginning stage of avoiding from a distant proximity. Until about 1679 seconds, the gravity is sufficient to change the moving direction of the spacecraft, leading the spacecraft to approach the goal. Therefore, between 1200 and 1679 seconds, although the spacecraft is provided with a thrust pointing towards the goal (which is the gravity generated by the artificial gravitational potential field), it moves away from the goal at a decreasing speed until the 1679 seconds when the spacecraft stops to approach the goal.Figure8 shows the schematic diagram of the active spacecraft transferring orbit with modified artificial potential field method. The diagram covers the stage of avoidance and the stage of returning from a distant proximity.Figure 8
Active spacecraft transfers orbit based on the modified artificial potential field method.Figure9 shows how the distance between the active spacecraft, and the unreachable domain varies during the maneuver, while Figure 10 illustrates how the distance between the active spacecraft and the goal changes.Figure 9
The variation of distance between active spacecraft and unreachable domain.Figure 10
The variation of distance between the active spacecraft and the goal.In Figure9, the closest distance between the active spacecraft and the unreachable domain is about 3.3 km; so, there is no risk of collision. As shown in Figure 10, the active spacecraft enters the stage of avoidance from 173 seconds to 1200 seconds, but it gradually moves away from the goal due to the impact of repulsive potential field. Following the stage from 1200 to 2446 seconds, the active spacecraft is under the gravitational potential field. However, during 1200 to 1679 seconds, the large inertia of the spacecraft prevents its motion state from an instant change caused by the gravitational potential field. As the spacecraft continues to move towards the goal, the distance between the spacecraft and the goal continues increasing. In spite of this, the increasing speed gradually slows down. At the 1679 seconds, the gravitational potential field has offset the inertia of the original motion when the spacecraft stops moving away from the goal but starts to approach the target. Afterwards, the distance between the spacecraft and the target narrows to the minimum at 2446 seconds when the spacecraft ends its stage of returning from the distant proximity.At 2446 seconds, the spacecraft will enter the stage of returning from a close proximity. The states of the spacecraft and the goal under the ECI coordinate system at the time are shown in Table3.Table 3
States of the spacecraft and the goal under the ECI coordinate system at the initial moment of the stage of returning from a close proximity.
x/my/mz/mvx/m⋅s−1vy/m⋅s−1vz/m⋅s−1Spacecraft41990424.0423692364.155989081.845-319.571-319.571-319.571Goal41990353.8153692445.699989387.824-278.757-278.757-278.757The CW coordinate system is established with the goal as the center of the circle. Substituting the initial values of the states in Table3 into Equation (42), the initial states of the spacecraft relative to the goal under the CW coordinate system is shown in Table 4.Table 4
The state of the spacecraft relative to the goal in the CW coordinate system at the initial moment of the stage of returning from the close proximity.
x/my/mz/mvx/m⋅s−1vy/m⋅s−1vz/m⋅s−1Relative state70.227-81.544-305.978-40.826-0.680-2.963It should be noted that since the CW coordinate system is established with the goal’s position as the circle center, the state of the goal’s position is0,0,0.Subsequently, the PSO algorithm carries out the simulation analysis on the path planning of the active spacecraft, and parameter values used in the PSO algorithm are(54)np=100,c1=2,c2=2,Vmax=10,d=3,kmax=100,ωmin=0.3,ωmax=0.75.For the search radius of the particle swarm, givenβ1=80,β2=0.02, β2=0.02, β3=0, and β4=−39, the relationship between spacecraft-goal distance and the particle swarm search radius is as shown in Figure 11.Figure 11
Relationship between the spacecraft-goal distance and the particle swarm search radius.In Equation (45), for the function to be optimized, the parameters are α1=0.8 and α2=0.2. The maximum pulse maneuvering volume for each maneuver of the spacecraft is defined as Δu<20m/s. In the simulation, the spacecraft performs pulse maneuvering at a time interval of 1 second. It takes the calculation of pulse maneuvering size and direction for about 0.48 seconds per time on average, while the maximum time is about 0.52 seconds, meeting the real-time requirements of the mission.The PSO algorithm also solves the maneuver mode and movement trajectory of the spacecraft in each iteration, and the results are shown in Figures12–14. Figure 12 shows the trajectory of the spacecraft to the goal, Figure 13 demonstrates the variation curve of distance between the spacecraft and the goal, and Figure 14 illustrates the spacecraft’s total pulse of the orbital maneuver at the stage of approaching from the close proximity.Figure 12
The trajectory of the spacecraft to the goal.Figure 13
The distance between the spacecraft and the goal: (a) distance in theX direction, (b) distance in the Y direction, (c) distance in the Z direction, and (d) total distance.
(a)(b)(c)(d)Figure 14
The time variation curve of the spacecraft’s total orbital maneuver pulse at the stage of approaching from a close proximity.As shown in Figures12 and 13, the PSO algorithm searches radius with dynamic variation helping the spacecraft return to the goal in about 13 seconds, thus effectively transferring the original orbit of returning from a close proximity. It should be noted that as the initial velocity of the spacecraft relative to the goal is too large, the motion curve of the spacecraft might not be smooth at the initial stage. The relative distance between the spacecraft and the goal in the O′X′ direction of the CW coordinate system increases before shrinking. However, in Figure 14, maneuver pulse amount is always less than 20 m/s during the process of the spacecraft maneuvering.As shown in Figure13, both the spacecraft-goal distance and velocity are close to 0 at about 13 seconds. Similarly in Figure 14, the pulse amount required for orbital maneuvers is around 0 after 13 seconds. Based on the two conditions above, the improved PSO algorithm can successfully push the spacecraft to back to its original orbit at about 13 seconds; so, no further maneuvers will be needed afterwards. As a result, the spacecraft is managed to continue its flight along its original orbit while being successfully prevented from the collision with the debris group.To sum up, this paper proposes an effective maneuver strategy for spacecraft to avoid space debris group and return to the original orbit. In existing studies, studies on how spacecraft evade space debris groups are not sufficient, as most research focuses on space debris. The simulation analysis effectively describes the unreachable domain of the space debris group by reducing the dimensionality of the changing geometry configuration of the space debris group. To calculate the avoidance orbit and return orbit of spacecraft, the APF method and PSO algorithm are adopted in the analysis. During the process, the APF method and PSO algorithm are improved. The APF method is suitable for spacecraft with high velocity by distance projection, while the calculation efficiency of the PSO algorithm is improved by dynamic searching radius.


## 6. Conclusions

This paper proposes a method of orbital transfer for spacecraft to urgently avoid the debris group in the near space and return to its original orbit. The method innovatively highlights the following: (1) the changing geometry of the space debris group is simplified by dimensionality reduction while the mathematical method describes the space debris group as an unreachable domain, (2) the method figures out a maneuver strategy covering three stages from avoidance to return, (3) the method suggests to use distance projection approach to break the limitations of the artificial repulsive potential function based on the special features of the spacecraft motion at the stage of avoidance, and (4) the method proposes a modified function to search with dynamic changing search radius and particle correction approach for a spacecraft at the stage of returning from a close proximity, which ensures the random walk of the PSO algorithm in the initial phase and the computational efficiency in the ending phase, to avoid local oscillation.



---

*Source: 4650730-2022-04-13.xml*