Time-Fixed Glideslope Guidance for Approaching the Proximity of an Asteroid

: The guidance and control problem of spacecraft approaching an asteroid using constant continuous thrust is studied in this work. The range of interest is from hundreds of kilometers to several kilometers, in which relative measurements of much higher accuracy than based on Earth can be used to facilitate further hovering or landing operations. Time-fixed glideslope guidance algorithm is improved by introducing a substitute of an existing control parameter and combined with elliptical relative orbital dynamics to rendezvous the spacecraft with a prescribed location in the proximity of a given asteroid. A vast range of values for the control parameters are explored and suitable combinations are found. To fully validate the robustness and accuracy of the proposed control algorithm, Monte Carlo simulations are done with the navigational error and implementation error considered.


Introduction
With an increasing interest in asteroid exploration, research on the design and control of trajectories towards the close vicinity of asteroids has gained much attention [1][2][3]. One of the major issues is how to guide the spacecraft from hundreds of kilometers to the vicinity of a given asteroid. This phase relies not on Earth-based telemetry, but on relative measurements and autonomous control by the spacecraft itself. It paves the way for even closer observational phases or landing [4]. Since the key problems and conditions are identical to the near-Earth rendezvous problem [5], where constraints are imposed on both the final relative position and velocity, the control issue of the transfer to the vicinity of an asteroid can be regarded as a rendezvous problem in the heliocentric two-body dynamics background, where the Tschauner-Hempel (TH) equations [5] provide good, approximated solutions to the first order for elliptical orbits. It should be noted that the relative range under consideration is from hundreds of kilometers to several kilometers, which belongs to the classical near-range rendezvous process. The aim for this stage is not to land the spacecraft on the asteroid, but to make observations or get preparation at a required distance. The sphere of influence of the asteroid in this work is estimated to be less than 90 m, as will be mentioned later. Therefore, the asteroid's gravity is ignored in this study.
Glideslope guidance is a multi-impulse trajectory transfer algorithm and has been applied to near Earth rendezvous studies. Li et al. [6] applied the glideslope method to asteroid approaching. Hablani et al. [7] proposed algorithms based on traditional glideslope guidance with astronauts in the loop, for spacecraft to approach, to fly around, and to depart from a target vehicle in a near-earth circular orbit. The closed-form solution of the CW equations was utilized, and this was based on impulsive assumption. The flight time therein is free to vary. The major drawback of previous methods is the total time of flight cannot be known before its implementation. Lian et al. [8] improved the method to make the time of flight as one of the control parameters and extended its application to libration point rendezvous problem. However, one of the key parameters, * ρ , which denotes the final distance-to-go, relies largely on specific rendezvous scenario, and may have different best values for different problems. In this work, we relate * ρ to the initial distance and replace it with a newly defined parameter.
The rest of the work is organized as follows. Section 2 describes the relative dynamic model and state space form of the classical TH equation is given. Section 3 improves the time-fixed glideslope guidance algorithm and combines it with the TH equation. Practical issues concerned with constant thrust are discussed. Based on the large number of numerical simulations, the efficiency and robustness of the proposed method are analyzed in Section 4. Section 5 concludes this study.

Relative Equations of Motion
Denoting the heliocentric position vector of the spacecraft by 1 r , and that of the asteroid by 2 r , the relative position vector (denoted by ρ ) of the spacecraft is computed Using the two-body dynamical equations, one has In this work, the development of the guidance algorithm is based on the orbital frame whose origin is at the center of mass of the asteroid, x-axis aligns with the heliocentric position vector, y-axis lies in the orbital plane and points towards the motion direction, and z-axis aligns with the orbital moment of momentum.
The approximated relative equations of motion are given by [9]  is the control acceleration vector, while ω and ξ are the angular velocity and angular acceleration, respectively.
Consider the following change of variable, (4) and substitute them into Equation (3), then we have the so-called T-H equation [5]. For simplicity, denote         , , , , , T x y z x y z X (6)        1  2  3  4  5  6 , , , , , The solution to Equation (3) [5] can be expressed as  A X D (8) where  (9) and where , , e f E are the eccentricity, the true anomaly, and the eccentric anomaly of the target, respectively,

Glideslope Guidance Based on TH Equation
The time-fixed glideslope guidance algorithm was initiatively proposed in [8]. For the completeness of this work, the major steps are reviewed in the following, where improvements are made in the computation of the control parameters.

Guidance Law Parameter Computation
Assume the spacecraft has an initial position, Note that * where  ρ ρ ρ , ( ) ρ t is called the distance-to-go, and should satisfy the following con- Design the following linear relationship connecting ρ and  ρ , Out of safety concern, it is required that Combining Equations(18)-(20)the following equations can be solved: Denoting the total number of segments of the trajectory by N , there will be 1 N thrust firing. Denote the final distance-to-go by * ρ (i.e., the last segment), and define According to [8], if From Equation (26), it is required that In this work, we introduce a ration factor  1 ε , then To sum up, there are four independent design parameters, T, N, 0 ρ , and ε , that fully determine the time-fixed glideslope guidance algorithm.

Guidance Law Design and Continuous Thrust Implementation
Given the flight time, T , the guided trajectory will be divided into N segments or arcs, for which the flight time is . According to Equations (16) and (21), From Equation (14), one has Denoting the current velocity by In real applications, the implementation of the theoretical instantaneous delta-v needs a period of time via a constant continuous thrust. Denote the firing time duration , which can be deducted based on the Tsiolkovsky equation and computed by , c is the exhaust velocity, and i m is the spacecraft mass before firing. The fuel consumption is computed by where F is the constant thrust magnitude.
The direction of the thrust is computed by Since the installation of the thrusters or the determination of the spacecraft attitude contains inevitable errors, both the magnitude and direction of the thrust may not be implemented fully as required. The impact of such errors will be considered in the following simulations.

Baseline Trajectory
In the process of approaching the asteroid, the spacecraft can be assumed as a point mass whose size and shape are negligible. The target asteroid is chosen to be 2000 SG344 based on an accessibility analysis whose details are out of the scope of this work and therefore will not be given. The asteroid has an estimated diameter of 40 m [10]. Unfortunately, the mass of the asteroid has no data available yet. However, if we assume the density of 2000 SG344 is the same as 433 Eros, one asteroid that has been closely visited by the NEAR mission and has an approximate gravitational constant of 4.46 10 m /s and a reference radius of 16 km [11]. Considering the mass is in cubic relationship to the radius, the sphere of influence of the asteroid 2000 SG344 can be estimated to be only 74 m. Since our final required location is kilometers away from the asteroid, its gravity thus can be ignored. The initial heliocentric state of both the asteroid and spacecraft are given in Table  1, where the last column means the relative state of the spacecraft with respect to the asteroid in the heliocentric ecliptic inertial (HEI) frame.
The values of the mass, the magnitude of the constant thrust, and the specific impulse are given in Table 3. Using initial conditions given in Tables 1-3 . Note that for ε , we only sample it at a 1 6 -step, although it is continuous within the range (0,1) . The time-fixed glideslope guidance algorithm will generate accordingly grids of data for the terminal position deviation ( e r ), the terminal velocity deviation ( e v ), the total fuel consumption ( Δm ), and the total delta-v ( ΔV ), whose relationships are shown in Figures 1-3  ε .   . The data is reorganized to show in Figure 4, where one can first notice that the cost indices ( Δm and ΔV ) seem independent with N because all curves overlap with each other, and as TOF grows, Δm and ΔV become smaller. Therefore, we take TOF = 40 h. Then we seek to select the value of N from the upper two subplots. Since e v is smaller than 1.5 × 10 −7 m/s for all values of N, we mainly focus on the e r subplot. It should be noted that we do not hope the number of segments is too large because too many thrust firings bring about larger probability of failure. Seen from Figure 4, the curves gather except for N = 2 which has obvious larger position deviations for all TOF values. Starting from N = 4, the differences among curves are negligible. Therefore, we choose N = 4. To sum up, taking the design parameters as  2 3 ε ,  4 N , and  TOF 40 h , the resultant terminal position deviation ( e r ) is less than 0.01 m, terminal velocity deviation ( e v ) is less than   9 2 10 m/s, the fuel consumption ( Δm ) is less than 1 kg, and the total deltav ( ΔV ) is less than 2 m/s. The approaching trajectory in the orbital frame and the associated relative distance are given in Figures 5 and 6, respectively, where in Figure 5 the red crosses denote the locations where the constant thrust is on, and the green circle displays the required location. Table 4 gives the exact value of each delta-v along the routine, Δ i V   As can be seen from Figure 5, the relative trajectory seems to follow a straight line. This is mainly because in a decade of hours, the orientation of the orbital reference frame changes merely about 0.47 deg (due to a very small orbital angular velocity    1. 3 5 deg/ s f E ), which makes the orbital frame could be viewed as an "inertial" frame, where the spacecraft experiences an almost force-free environment because the 100-km-distance is too small to make a large difference between the Sun's gravitational accelerations on the spacecraft and the asteroid, respectively. Also noting that in Table 1 the initial inertial velocity is zero (the associated relative velocity is close to zero as shown in Table 2), the glideslope guidance algorithm therefore generates a path almost coinciding with a straight line determined by the Newton's first law.
If the initial relative velocity (inertial or non-inertial) takes a non-zero value,   [50.0, 60.0,80.0] T r V m/s, for instance, the trajectory in the orbital frame will first have turning segments to diminish the relative velocity, then approach the target location along the line of sight, as shown in Figure 7. The first thrusting arc (marked by red crosses) takes a much longer time than other arcs to remove the relative non-zero velocity. Note that the trajectory ever since the completion of the second thrusting starts to align with the line of sight.
If we enlarge the time of flight to the order of decades of days, 100 days, for instance, the trajectory will be no longer straight lines. For the purpose of comparison, we still use the zero-initial inertial relative velocity. The result is given in Figure 8, where one can see the trajectory does not follow the line of sight, because the flight time is long enough (25 days for each arc) to make large enough changes of the difference between the orbital frame and the inertial frame. Approximations to the conditions in Newton's first law no longer hold and curvature arcs appear in the orbital frame.

Monte Carlo Simulation
In real practice, various errors exist that will have inevitable impact on the guidance performance. As a matter of fact, guidance problem is proposed to deal with errors. In this work, navigational error and implementation error are taken into account, which are assumed to follow a zero-mean normal distribution. Monte Carlo simulations will be done to obtain statistical results.
Navigational error, including position error and velocity error, comes from limitations of both hardware and software which result in difference between the measurement and the real value. In space rendezvous applications, navigational error is related to the relative distance, i.e., the closer the distance, the higher the accuracy. Implementation error comes from the misalignment of the engine or attitude control error and will be reflected in both the magnitude and direction of the thrust. Since the glideslope algorithm is in essence an impulsive guidance method, the implementation error will be simulated based on delta-v. Assume the above errors follow zero-mean normal distributions with different standard deviations, which are given in Table 5, where rj σ corresponds to the position error, vj σ the velocity error, Δvj σ the implementation error, and  , , j x y z . As shown in Table 5, there are four combinations of different values of the standard deviations, in order to analyze the individual impact of the errors.  Table 5, and obtained results are given in Tables 6-9. In order to check the individual impact of each error, three conditions are taken into account, (1) only the navigational error, (2) only the implementation error, (3) both errors. The maximum, mean, minimum, and standard deviation of each 300-run is given in Tables 6-9. It should be noted that in real practice, it is the maximum value, not the mean value, of each performance index that should be paid most attention to. For the comparison convenience, results of the baseline trajectory are also given in Tables 6-9. For error set I, one can see from Table 6 and Figure 9 that, the terminal position deviation can reach 183 m if all errors are taken, while the maximum of the terminal velocity deviation is 0.006 m/s, and 1.961 m/s and 0.938 kg for ΔV and Δm , respectively. Small standard deviations (0.009 m/s and 0.004 kg) implies that the navigational error and implementation error have little impact on ΔV and Δm , but on terminal position instead. Another fact is that the last corrective control is very effective to diminish the terminal velocity deviation.  For error set II, results are given in Table 7 and Figure 10. Compared with error set I, the only change is the standard deviation of the position navigational error is enlarged from 0.1 m to 1 m. One can see that the variations of the control performance indices when compared with their counterparts in Table 6 are negligible, which implies that the terminal position accuracy is not sensitive to the position navigational error.  For error set III, results are given in Table 8 and Figure 11. Compared with error set II, the standard deviation of the velocity navigational error is increased from 0.001 m/s to 0.01 m/s. Results show that the terminal position deviation grows from a maximum of 183 m to 1248 m, the terminal velocity deviation from 6.453 × 10 −3 m/s to 3.534 × 10 −2 m/s, while the fuel consumption and total delta-v show no obvious change, which indicates that the velocity navigational error has a dominant impact on the terminal position deviation, while the last correction remains effective to the terminal velocity deviation.  For error set IV, results are given in Table 9 and Figure 12. Compared with error set II, the only difference is the implementation error increases from 0.5% to 2%. As a result, the terminal position deviation grows from 183 m to 594 m, the terminal velocity deviation increases from 6.453 × 10 −3 m/s to 2.511 × 10 −2 m/s, and the fuel consumption and total delta-v remain almost the same. This indicates that the implementation error has a very large impact on the terminal position deviation, while the other indices are little affected.  To sum up, the terminal position deviation ( e r ) is the only performance index sensible to the navigational error and implementation error, while the terminal velocity deviation ( e v ), the total fuel consumption ( Δm ), and the total delta-v ( ΔV ) are much less affected by the variations of the errors. In the case of 0.1 m (standard deviation) of position uncertainty, 0.001 m/s of velocity uncertainty, 0.5% of implementation uncertainty, e r can be less than 185 m, e v less than 7 × 10 −3 m/s, Δm less than 1 kg, and ΔV less than 2 m/s. Velocity navigation error is the most powerful factor that destructs the final position accuracy. A standard deviation of 0.01 m/s of the velocity navigation uncertainty will result in almost 1300 m of final position deviation.

Conclusions
This work initiatively combines the elliptical orbital relative dynamics and time-fixed glideslope guidance method to study the control problem of approaching the proximity of an asteroid. A substitutional ratio factor is proposed to improve the design flexibility of time-fixed glideslope guidance law. A large number of Monte Carlo simulations are done with taking into account both the navigational error and the implementation error. Results show that different combinations of the two types of error mainly affect the terminal position accuracy, while the terminal velocity accuracy, the fuel consumption, and the total delta-v vary in narrow ranges. For the spacecraft that maneuvers from 100 km to 1 km of relative distance, the total delta-v is approximately 2 m/s, and fuel consumption is around 1 kg for a spacecraft of 1-ton mass and 300-N constant thrust. The velocity navigational error and the implementation error play key roles in affecting the terminal position accuracy. When the velocity navigational error increases from 1 mm/s to 1 cm/s, the maximum terminal position deviation changes from 184 m to 1248 m. On the other hand, when the implementation error increases from 0.5% to 2%, the maximum terminal position deviation grows from 184 m to 594 m. In all error settings, the terminal velocity deviation is less than 0.035 m/s, indicating the last corrective delta-v is very effective.
In conclusion, the proposed method proves to be promising in the real application of asteroid proximity approaching, because it only requires a small number of thrust burns with low fuel cost, and the terminal state accuracy is good enough to initialize following a hovering or landing phase of the spacecraft.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy.