Spacecraft relative guidance via spatio-temporal resolution in atmospheric density forecasting

Spacecraft equipped with the capability to vary their ballistic coef ﬁ cient can use differential drag as the control force to perform propellant-less relative maneuvers. Because atmospheric drag is proportional to atmospheric density, uncertainty in atmospheric density makes the generation and tracking of drag-based guidances dif ﬁ cult. Spatio-temporal resolution, or the mapping of density information to altitude and time, is shown in this work to improve atmospheric density estimation from forecasted density for spacecraft in LEO. This is achieved by propagating simulated orbits for two spacecraft using forecasted density. Additionally, a receding-horizon control algorithm is introduced, with the goal of improving the tracking of guidances. Using a simulated perfect forecast of the atmospheric density for propagation of the orbits, relative guidance trajectories are generated and tracked, establishing the bene ﬁ t of adding spatio-temporal resolution. Next, imperfect density forecasting is added, indicating that the bene ﬁ t of spatio-temporal resolution is retained in the presence of imperfect forecasting. Finally, a receding-hor-izon control algorithm is used with imperfect forecasting, demonstrating that receding-horizon control improves the tracking of guidances compared to single-horizon control.


Introduction
Formations of small satellites hold the potential for replacing large complex spacecraft, as explained in Refs.[1][2][3][4].On-orbit inspection and maintenance missions and other complex space tasks can be performed by spacecraft flying in formation, providing redundancy in the case of a loss of a satellite.Additionally, smaller satellites are lighter and can be launched as a secondary payload for existing missions, which reduces the cost of orbit injection [5][6][7].Consequently, there is a growing interest in the aerospace community in the development of methods for small spacecraft autonomous formation flying.Formations of spacecraft can cover more ground tracks than any single spacecraft, which may comprise identical or different orbits.Since the atmospheric density varies with both location and time, this implies that spacecraft in a formation will experience different atmospheric density.
Any formation of spacecraft requires the ability for the spacecraft involved to control their relative motion, typically performed using thrusters, requiring propellant to be carried aboard.Hence, alternative means to maneuver spacecraft are of great interest.Leonard et al. [8] proposed varying the cross-wind area of spacecraft to alter the drag force acting on them, as a method for controlling their relative motion at LEO. Differential drag can allow for propellant-less planar relative maneuvering, which can reduce fuel usage and costs for formation flying missions.Sensors mounted onboard spacecraft can also benefit from a cleaner environment due to the lack of thruster plumes.However, using differential drag to maneuver imposes the constraint of operating where the atmosphere is sufficiently dense to generate significant drag forces, and limits the maneuvers to the orbital plane.Moreover, using the drag for maneuvering increases the orbital decay rate of the spacecraft.Despite the downsides of differential drag, future impacts of the ideas here proposed can be foreseen for higher orbits.For example, the concept of exploiting differentials in environmental forces can be also imagined for geosynchronous satellites using solar radiation pressure [9].
In the last few years there have been quite a few papers inspired by Leonard et al [8].Bevilacqua et al. [10,11] used the linear Schweighart and Sedwick model [12] to create a differential drag based rendezvous guidance assuming constant density.Ben-Yaacov and Gurfil [13] studied the use of differential drag for clusterkeeping purposes.Pérez and Bevilacqua [14] developed a Lyapunov-based controller for relative maneuvering of spacecraft using differential drag.Dell'Elce and Kerschen proposed the use of model predictive control [15] and a three-step optimal control approach [16] for drag based rendezvous maneuvers.There have also been a few efforts for exploiting the differential drag concept in real missions.The ORBCOMM [17] constellation used differential drag for constellation keeping.Also, the JC2Sat [18][19][20][21] project developed by the Canadian Space Agency (CSA) and the Japan Aerospace Exploration Agency (JAXA) proposed the use of differential drag for relative maneuvering of spacecraft within close proximity of each other, extending the methodology presented by Leonard and studying implementation issues such as navigation errors.Finally, work by Mazal et al. has focused on using differential drag for long-range maneuvers in the presence of uncertainties in the control forces used [22].
Difficulty in estimating the drag force results in lack of realism in any drag-based guidance trajectory, making tracking more difficult.In the literature on drag-based maneuvering, it is usually assumed that the density is constant for guidance and control purposes (see Refs. [11,17,[23][24][25]). Any guidance trajectory created under the assumption of constant density will be inaccurate due to unrealistic control forces.
In Ref. [26], density from an existing atmospheric model was used for creating a guidance trajectory for a drag based rendezvous.Forecasting was not used; the density was assumed to be known.Pérez at al. subsequently developed an Artificial Neural Network (ANN) forecasting method for atmospheric density along the orbit of a non-maneuvering spacecraft (i.e., with a constant ballistic coefficient and no thrusters) [27].This forecasting method was later combined with the control methods by Pérez and Bevilacqua in [28].
Because the drag-based trajectory followed by a spacecraft is dependent on the atmospheric density it encounters, over or under-forecasting atmospheric density will result in uncertainty in forecasting the trajectory.When a spacecraft leaves the forecasted trajectory due to inaccurate forecasting, even a perfect forecast becomes inaccurate, since it corresponds to a different location.Adding spatio-temporal resolution to atmospheric density forecasting compensates for spacecraft leaving the forecasted trajectory.
In prior work, a single forecasted trajectory was used to create a rendezvous, where the density was forecasted and assumed to be along this trajectory [28].Since any deviation from the forecasted trajectory results in an inaccurate density forecast, and vice versa, using a single trajectory does not have sufficient information to provide a complete density and trajectory forecast.Adding multiple trajectories bounds the motion of the spacecraft.At any given timestep, there exists a set x containing the altitudes of each forecasted trajectory.The deviations between the members of x and the actual trajectory are contained in the set y.For the case of a single forecasted trajectory, both x and y have only a single member, which is by definition the minimum.Increasing the number of forecasted trajectories increases the number of members in both sets, and since the minimum of a set cannot be increased by adding more members, increasing the number of forecasted trajectories can only decrease the minimum deviation between a forecasted trajectory and the actual trajectory.
Previous work has considered atmospheric density as only time-dependent, and independent of location.Since real atmospheric density depends on both time and position, knowledge of a spacecraft's deviation from the expected trajectory can be used to improve the density forecast.This is denoted as spatio-temporal resolution, which reflects both the dependency of the density on both spatial and temporal differences.Using spatio-temporal resolution with an existing differential drag-based relative maneuvering algorithm [29], a rendezvous maneuver is created by modifying the aerodynamic drag on two spacecraft.
Using traditional control methods, or fixed-horizon control, a control algorithm will develop a set of control inputs, which are most effective when the system dynamics nearly match those of the real world.However, in this work, the state at a given timestep is dependent on all previous timesteps, and so small uncertainties in the system dynamics can rapidly result in large uncertainties in the state of the system, which can only be partially compensated with spatio-temporal resolution.Receding-horizon control allows the control algorithm to periodically reset the error state by restarting the control algorithm, which results in control inputs that are more applicable to the real world dynamics.
Receding-horizon control methods have been used previously for autonomous operation and decentralized control.Bellingham et al. [30] used receding-horizon control with fixed wing aircraft to account for disturbances along the flight, resulting in paths that very nearly minimized the chosen cost function, approximating closed-loop control.Additionally, Dunbar and Murray [31] presented an algorithm for multi-vehicle autonomous stable flight.
Receding-horizon control methods will be used when creating the maneuver to show the benefit of periodically updating the control to more accurately reflect the real world dynamics of the maneuver.Using receding-horizon control, a guidance is initially created for the entire maneuver, and updated as the maneuver progresses [32].After each predetermined guidance interval, the guidance is re-computed to account for density estimate inaccuracy.Decreasing the guidance interval approximates closedloop maneuvering as the guidance interval approaches the control update time.Solely closed-loop maneuvering was not used to create the guidance since measured atmospheric density is required for closed-loop control, and this is unavailable at the time of guidance creation.
This process is repeated for the duration of the maneuver.Using receding horizon control minimizes the error between the guidance and execution trajectories by periodically updating the beginning of the execution trajectory.The details of potentially implementing the algorithm onboard a spacecraft are also presented.The density field is uploaded to the spacecraft prior to the maneuver, and only the guidance and control calculations are made onboard.

Assumptions and reference frames used in this research
The guidance methodology presented in this work assumes that two spacecraft (chaser and target) have the capability to vary their ballistic coefficient by deploying or retracting a set of drag surfaces, thus changing the magnitude of the aerodynamic drag.
The reference frame used in this work for representing spacecraft relative motion is the Local Vertical Local Horizontal (LVLH) reference frame, which is assumed to be attached to the target spacecraft.Fig. 1 shows the right-handed Earth-Centered Inertial frame and LVLH frame.The LVLH directions are found from the state of the spacecraft in the ECI frame as follows: Quasi-circular orbits are assumed in this research.Given that the out of plane component of the aerodynamic force is assumed to be zero, and all aerodynamic drag acts antiparallel to the ŷ direction, atmospheric differential drag can provide effective control only in the orbital plane ( x and ŷ).Hence, the discussion presented in this work will be limited to in-plane motion, assuming that out-of-plane ( ẑ ) motion is controlled by different means.The attitude of the spacecraft is also assumed to be stabilized.Additionally, it is assumed that the control input can have three different configurations, either positive maximum (maximum drag on the target and minimum drag on the chaser), negative maximum (maximum drag on the chaser and minimum drag on the target), or zero (minimum drag on both spacecraft), as was previously done in Refs.[10,[19][20][21], neglecting the time required by the surfaces to be deployed or retracted.

Advances in the state of the art
Upcoming and ongoing spacecraft missions, such as the QB50 mission [33], benefit greatly from the ability to accurately anticipate on-orbit atmospheric density to maintain desired formations.Density information can be uploaded for subsequent maneuvers from a ground-based station, reducing the computational load on the spacecraft.
With improvements in density estimation and forecasting, more accurate density information can be uploaded to the spacecraft for subsequent maneuvers, improving estimation as the underlying data continues to improve.In this research, the density is assumed to be uploaded as three density time-series, pairing density and location information.Interpolating these timeseries allows the creation of a density field, which is used to create a relative guidance trajectory.The control inputs required to actuate this guidance trajectory are then applied to the spacecraft during the maneuver, the goal being to match the relative guidance trajectory.
Since density estimation will always have errors, open loop control will lead to errors in the estimation of the final relative position of the spacecraft, and by extension, the real trajectory will not perfectly follow the guidance, hence receding-horizon control is proposed.
This paper builds on work previously presented by the authors [34] and presents the following advances in the state of the art: The addition of spatio-temporal resolution to atmospheric density forecasting is shown to improve the accuracy of density prediction during a rendezvous maneuver, using a simulated perfect forecast.
Using a sample uploaded density field, spatio-temporal resolu- tion is shown to still be beneficial overall when used with imperfect density fields.
A method to implement the addition of spatio-temporal re- solution onboard spacecraft in orbit is developed, with the goals of minimizing the computational load on the spacecraft and improving the accuracy of the density forecasting for differential drag-based maneuvers.This paper is organized as follows.Section 2 details the relationship between atmospheric density and drag, and the effect of solar and geomagnetic effects on atmospheric density.Section 3 introduces the dynamics used in this research, and the algorithm for the creation of differential drag-based guidance trajectories.Section 4 demonstrates the algorithm performance under varying conditions.Finally, the conclusions are presented in Section 5.

Atmospheric density and drag relationship
The drag acceleration affecting spacecraft in LEO is a function of the atmospheric density and winds, orbital velocity, geometry, attitude, drag coefficient and mass of the spacecraft.Many parameters must be estimated, producing errors in the modeling of the drag force.This causes large uncertainties in the control forces available for drag-actuated maneuvers.The differential drag acceleration for target and chaser in similar flow conditions is expressed with the following equation: According to Vallado [35], this approximation of the drag acceleration holds true for a body in viscous flow, in which the time scale of collisions between the body and gas particles is much smaller than the time scale of particle-particle collisions.It is also assumed that the two spacecraft are in similar flow conditions.These equations do not address the effect of atmospheric winds, which affect the velocity of spacecraft relative to a medium.Additionally, uncertainties in the ballistic coefficient and density can affect the drag acceleration.

Solar and geomagnetic activity effects on the atmospheric density
Doornbos et al. [36] indicate that density can vary by more than one order of magnitude at a given time and location as a result of different solar and geomagnetic conditions.Specifically, increased solar and geomagnetic activity tends to increase atmospheric density.Walterscheid [37] provides an explanation of the increase in atmospheric density due to geomagnetic activity.Solar EUV (Extreme UltraViolet) radiation excites air molecules, with effects primarily in the region of high activity.In the regions of high activity, lower molecular weight air molecules absorb energy, causing them to rise to higher altitudes.These are replaced by particles of higher molecular weight, increasing the average local mass density in the area of high activity.An example of the effect of solar and geomagnetic activity can be observed in the evolution of the atmospheric density during the 2003 Halloween solar storm, as explained by Bruinsma et al. [38].

Estimating atmospheric density with density models
The atmospheric density can be estimated using atmospheric models, which can be classified into physics-based and empirical models.The physics-based models can be more accurate, but biases can result from poor representation of the underlying physics (see Ref. [39]), and are computationally expensive due to the high number of variables used and the required size of the modeled spatio-temporal volume.In contrast, empirical models, which rely on observed data, are less computationally expensive.However, since empirical data are not available for all locations and times, biases in the estimated density can result due to behavior not captured in the observed data (as shown in Ref. [40]).In this work three empirical atmospheric models, DTM-2013 [41] (Drag Temperature Model 2013), JB2008 [42] , and NRLMSISE-00 [43] (Naval Research Laboratory Mass Spectrometer and Incoherent Scatter, Exosphere 2000) are used for modeling the density.NRLMSISE-00 is available as a function within ® MATLAB , while the authors have made ® MATLAB wrapper functions avail- able for JB2008 1 and DTM-2013. 2he atmospheric density encountered by a spacecraft determines its drag acceleration.Accurately forecasting atmospheric drag allows a more accurate prediction of the expected trajectory of a spacecraft.The drag force also determines the de-orbit time of the spacecraft, which allows a better estimation of the number of possible maneuvers.Accurate forecasting of the expected atmospheric density onboard the spacecraft is required for autonomous operation.Because atmospheric density depends on location in the atmosphere, any uncertainty in the location results in uncertainties in the atmospheric density and therefore the atmospheric drag acceleration.These uncertainties result in uncertainties in the position of the spacecraft, resulting in further uncertainties in the atmospheric density.
Aside from solar and geomagnetic activity, altitude is the primary driver of density, and both the U.S. Standard Atmosphere [44] and International Standard Atmosphere [45] use altitude as the sole driver of density.At a given altitude, the atmosphere can be represented as a column of air, with the pressure dependent on the weight of the gas particles at higher altitudes.Eq. (7) shows the dependence of altitude on pressure: Using the form of the ideal gas law below, it can be seen that the atmospheric pressure is proportional to the density.Assuming that the composition does not significantly change with altitude, this indicates that the mass density will decrease with increasing altitude.ρ = ( ) P RT 8

Adding spatio-temporal resolution to atmospheric density forecasting
Because of the large effect of altitude on atmospheric density, the density at different altitudes is interpolated to produce a more accurate density estimate.This process is known as adding spatiotemporal resolution, shown in Fig. 2. Beginning from the current location of a spacecraft, the position and the velocity are propagated forward in time for several constant cross-wind areas, which results in several trajectories, shown as dashed lines.Subscripts A, B, and C refer to trajectories propagated with different drag configurations, where the subscript A represents that of the minimum drag case.Each trajectory has a different altitude profile, providing density information for each of these altitudes.When creating the guidance (solid line), the density from each trajectory is interpolated based on its altitude to obtain a single density estimate for each timestep.
The inset in Fig. 2 shows the interpolation data for a single timestep.The altitude for a given trajectory is denoted h timestep level , , and the density is denoted ρ timestep level , .Using the equation below, the density estimate ρ ^can be found as a function of the altitude of the guidance at the same timestep: In prior work by Pérez and Bevilacqua [28], a single forecasted trajectory was used, which corresponds to the orbit previously denoted B, in that the cross-wind area used was halfway between the maximum and the minimum.Since both spacecraft used in a rendezvous maneuver will modulate their cross-wind area throughout the maneuver, this implies that the time-averaged trajectory of each should be near the orbit B, since it represents the trajectory resulting from the average cross-wind area.However, this is not necessarily the case.Because spacecraft in a larger higher-altitude orbit will have a lower average angular velocity, a spacecraft that passes though point h A 1, in Fig. 2 at minimum drag and then increases the drag to descend to point h B 2, will reach this point at a later time than a spacecraft that passes through point h A 3, and then decreases its drag, due to the difference in angular velocity.A single-trajectory forecast does not take this time factor into account.
Increasing the number of trajectories captures more information about the state of the atmosphere.With a single trajectory, the density forecasting error is driven by the altitude error between the actual and forecasted trajectory.As mentioned in Section 1, by using multiple trajectories, the maximum altitude error between the actual trajectory and the set of forecasted trajectories is reduced, which reduces the density error.
Additionally, although the trajectories are closer together in the beginning of the maneuver, spatio-temporal resolution presents more of a benefit after several hours have elapsed because the maximum and minimum drag trajectories are further apart.In the beginning of a maneuver, when the trajectories have not yet diverged much, the maximum possible error in interpolation will be the lowest.All trajectories are very similar, so any interpolated density will be close to all the others.Conversely, at the end of the maneuver, after the trajectories have diverged, it is possible to interpolate the density at a point further from any individual trajectory than in the beginning.Hence, spatial resolution presents a larger benefit later in the maneuver.

Nonlinear orbital dynamics
The nonlinear orbital dynamics of the spacecraft, including two body, J 2 perturbation, and drag acceleration, are represented using the following expression: where B is defined in Eq. ( 6).These nonlinear dynamics (taken from Alfriend et al. [46]) are used to propagate both the target and chaser when tracking and executing the guidance.

Lyapunov controller
In Ref. [47] a Lyapunov approach was used to determine the control input.With Lyapunov control, the derivative of the Lyapunov function V is driven negative, with the goal of eventually driving the function to zero.The Lyapunov function is defined as follows: V is unitless.The control law presented in Ref. [47], necessary to ensure the stability (by driving V negative) is as follows: where x represents the current state of the system.The gain matrix P is a result of the LQR equation.A negative control signal û is interpreted as a signal to increase the area of the chaser relative to the target, and a positive control signal is interpreted as the opposite.A zero control signal can only result from a zero tracking error in relative position and velocity, and so is interpreted as a signal to reduce the area of both spacecraft to preserve the orbit.The error state is defined differently when creating and tracking the guidance.For creation of the guidance, since the objective is to create a rendezvous, the error state is defined as the following: For tracking of the guidance, the error state is defined as follows, where the subscripts g and tr indicate the guidance and tracking states, respectively: Using spatio-temporal resolution and interpolation to improve density estimation.

Creation and execution of differential drag-based relative guidances
The goal of guidance creation is to create the relative path the chaser spacecraft will take to rendezvous with the target.Fig. 3 shows an example guidance (thin line), beginning at an arbitrary location with the goal of completing a rendezvous.This guidance represents the relative position and the velocity of the chaser spacecraft with respect to the target, so reaching a zero-position, zero-velocity state represents completion of a rendezvous.Execution of the guidance (thick line) shows an example of deviations due to uncertainties in the atmospheric density forecast.
The density field is created prior to the maneuver epoch, and provides density information for each potential location-timestep pair.The density timeseries are created using the forecasting techniques described in Section 4. The locations of these timeseries points are determined by propagating the state of the spacecraft using the appropriate cross-wind area.
While creating the guidance, density is interpolated at each maneuver timestep, based on the altitude.After each control update, the error state is used to determine the control inputs to each spacecraft, which determines the cross-wind area used based on the previously described controller.This process is repeated for the length of the horizon.Execution of the guidance is accomplished using the same set of control inputs that were used when calculating the guidance.
After each horizon, the guidance calculation and execution is repeated, starting from the end position of the previous guidance execution.As the maneuver progresses, the guidance for the remainder of the maneuver becomes shorter, since part of the maneuver has already elapsed.Each guidance is executed for the length of the horizon.

Implementation of the algorithm onboard a spacecraft
Existing density models require significant computational capacity.In addition to the three empirical models used in this research [41][42][43], physics-based models require even more computation.For instance, Global Ionosphere-Thermosphere Model (GITM) is intended to be run on many processors by splitting the atmosphere into a grid [48].Additionally, all three empirical models require atmospheric indices to modify the density estimates produced.
Since electrical power onboard a spacecraft is generally limited to what the solar panels can produce, all components must use minimal electrical power, or limit the run-time of calculations.An average of 1367 W/m 2 is generally used in calculations [49], which incorporates the varying distance from a spacecraft to the Sun, as well as the varying solar irradiance.Multiplying by an estimated 15% efficiency gives the maximum power for a representative solar panel area, approximately 205 W/m 2 .As an example, with ap- proximately 0.01 m 2 area, the Clyde Space 1U panels closely track the expected power production at 2 W [50]. Since the solar panels can only produce electrical power when oriented toward the Sun, and when not eclipsed by the Earth, 100 W/m 2 is a more con- servative estimate.Using this estimate of power density, small spacecraft can be assumed to produce only a few watts.Additionally, since convection in orbit can be assumed to be negligible, thermal radiation dominates the heat transfer of spacecraft in LEO.All heat must be dissipated through radiation; limiting the power dissipated as heat minimizes the temperature variation of spacecraft and prevents premature failure of onboard components.
Small spacecraft generally use low-power, legacy hardware for processing.Reliability is paramount on spacecraft since there is no way to manually reset the hardware.Radiation absorbed while in orbit has the potential to affect processing and flip bits in memory [51], causing errors.Some techniques are available to mitigate the effects of radiation, including error-correcting code (ECC) memory and redundant processors, but it is preferred to avoid the problem by using proven hardware with flight heritage.Hence, the processors used on spacecraft often have a low transistor count and are several years out of date, not including possible launch delays.
In order to reduce the computational burden of the spacecraft, the density fields for future orbits of the spacecraft must be uploaded to the spacecraft prior to the maneuver.Equipping the spacecraft with a processor that can calculate the density onboard still requires regular uploads of the atmospheric indices, and it is simpler to upload the density directly.Both spacecraft will therefore relay their positions to the ground station rather than  Fig. 4 shows the steps involved.The target spacecraft is shaded more lightly than the chaser spacecraft.Solid lines indicate the actual path of each spacecraft, and the dashed line indicates the guidance trajectory.
The procedure can be summarized as follows: 1.Both spacecraft send their location to the ground station.
2. The density timeseries are computed on the ground.
3. The appropriate density timeseries is sent to each spacecraft.4. Each spacecraft computes the guidance necessary to create the rendezvous, interpolating the density at each maneuver timestep.5.The guidance may be re-computed after some time has elapsed if receding-horizon control is used.
After receiving the location information from each spacecraft, the density timeseries for each is calculated at the ground station.Once computed, the density timeseries are then relayed back to the target and chaser spacecraft.
Once each spacecraft has received the density timeseries, the guidance for each can be calculated.This guidance will then be followed, either for the duration of the maneuver or until the first control horizon has been reached, depending on the control method chosen.Calculating the guidance in this manner leverages the ground station computing capability while still retaining autonomous operation.
An alternate criterion for recomputing the guidance would be to recalculate once a predetermined error has been reached.However, deviations from the guidance can rapidly accumulate, resulting in control horizons of varying durations.This may result in an increased processing requirement if the duration is allowed to decrease without a lower bound.Since a lower bound is necessary, it is simpler to use a fixed duration control horizon.

Numerical simulations
Before demonstration of the algorithm, the first step is to demonstrate that adding spatio-temporal resolution will in fact improve the density forecasting.To that end, effects due to imperfect forecasting of density will be neglected, but interpolation between the reference time-series will be added to provide the spatiotemporal resolution.By propagating the state of the spacecraft forward, using DTM-2013 directly for density information, perfect forecasting is simulated.DTM-2013 is used due to its improved accuracy over NRLMSISE-00 and JB2008 [41].
After verification of the benefit of adding spatio-temporal resolution, the next step is to consider density forecasting.The density forecast is represented as three density and location timeseries, which can be assumed to be uploaded from the ground station as the spacecraft pass over it.This density timeseries is generated using a similar density model, JB2008, which requires the atmospheric indices F10.7 and Kp.These are assumed to be forecasted using existing sources, such as those forecasted by NOAA [52].Using the uploaded density timeseries, the density is interpolated at each timestep to produce the estimated density field.
Finally, after verification that use of the uploaded density field with spatio-temporal resolution improves the density forecast, the next step is to add in receding horizon control, which is expected to reduce the tracking error.

Simulation parameters
The initial conditions for the target and chaser spacecraft are shown in Table 1.Parameters that differ are bolded.Both spacecraft begin in the same orbital plane to eliminate out of plane effects.These are based on initial conditions used in similar work by Pérez and Bevilacqua [29].Converting to the target LVLH frame, these represent a relative location of approximately À 1 km in the LVLH x direction, and À2 km in the LVLH y direction.
A maneuver length of two days is used.With a maneuver timestep of 60 s, this represents 2880 timesteps.Based on the maneuver initial conditions, using a two day maneuver allows sufficient time to create a rendezvous.A control update (time between control changes) of 10 min is used.This allows sufficient time to propagate the change in orbit, while still allowing sufficiently rapid control updates.Two maneuver epochs are used in this research, due to the effect of geomagnetic activity on the atmospheric density.The high activity case has a maneuver epoch of January 18, 2005, 00:00:00.000UTC.This is chosen because it represents an interval of sustained exceptionally high K p and F10.7 activity, with means of 59.25 and 120.30, respectively (both are unitless).Additionally, the low activity maneuver epoch is January 27, 2005, 00:00:00.000UTC.Similarly, this represents an interval of sustained low activity, with means of 10.06 and 83.35, respectively.This is done to bound the cases that are likely to be observed.
Fig. 5 shows the activity during January 2005, with the high and low activity cases boxed in green (dashed) and red (dot-dashed), respectively, indicating the range of activity.Data are from NASA OMNIWeb [53].
By definition, spatio-temporal resolution is possible with two or more trajectories.However, adding more trajectories involves a tradeoff between improved forecasting precision and increased computation time when generating the density.Adding more trajectories increases the precision but does not necessarily increase the accuracy; this depends on the accuracy of the density model used.For these reasons, three trajectories were used in each case.
Only the J 2 perturbation was included in the orbit propagation.Although both orbits will have similar perturbations due to their relative proximity, the J 2 perturbation is necessary for accurate simulation of the system dynamics.J 3 and higher order perturbations generally are disregarded, as these are several orders of magnitude smaller than J 2 [49].
Both spacecraft used are assumed to have repeatedly retractable drag surfaces, with an area range of 0.174-5.68m 2 , for a ratio of approximately 33:1, which is in line with existing space sail missions such as NanoSail-D [54] and IKAROS [55].The different drag accelerations separate the possible trajectories, maximizing the benefit of spatio-temporal resolution.Since there are three trajectories used to create the density field, these represent the minimum drag case, maximum drag case, and a third case with the average of the minimum and maximum drag (which is not realizable in practice). 3Additional trajectories would fill in the gaps between these three.

Simulated perfect forecasting guidances
Using DTM-2013 directly as the source for density information, which represents using the actual in-orbit density, the state of each spacecraft is propagated for the duration of the maneuver.This is repeated for each of the three trajectories, creating a density timeseries for each.
After generating the density timeseries, the maneuver can be simulated.By interpolating the density at each timestep based on Eq. ( 9), field can be produced, providing a density estimate for any point in the atmosphere.By creating a density estimate of positions along the possible orbits, a guidance can be created for a rendezvous maneuver.The previously described Lyapunov controller (Eq.( 13)) is used to update the cross-wind area of each spacecraft after each control interval.This is performed with the goal of driving the spacecraft to a rendezvous.After the rendezvous is reached, the orbits are propagated with the minimum area to prevent any further relative maneuvers.
The guidance control profile is executed in open loop for the remainder of horizon.With a perfect density forecast and interpolation, the density encountered along each trajectory will be the same between the guidance and execution, although the interpolation will still be an approximation.Assuming the real world dynamics are the same as in the guidance, the relative trajectory reached when executing the maneuver will be the same as that of the guidance.With imperfect forecasting, the relative trajectory will be slightly different, but improvements in the forecasting accuracy will bring the execution relative trajectory closer to that of the guidance.The density information used to execute the maneuver is taken directly from DTM-2013, as an analog to using in-orbit density.The Lyapunov function represents a weighted error state, accounting for both the relative position and velocity error.
Four guidances were created, three using a single forecasted trajectory (representing the previous state of the art [28]), and one using the interpolated density from all three trajectories, for the density estimate.The control time-series from each of these is then used to execute the maneuver.
Table 2 denotes the errors obtained when comparing the guidance and execution relative trajectories.The time averaged Lyapunov function is used to indicate the performance, with a lower mean function indicating better performance.In both the high and low activity cases, the time-averaged Lyapunov function is lowest when using the interpolated density, indicating these guidances were the closest to the execution relative trajectory.
Extending the maneuver by increasing the time interval would represent maintaining an already completed rendezvous.Because the rendezvous is completed at this point, the relative position and velocity will be nearly zero, and so averaging in more timesteps after the rendezvous reduces the overall average.Therefore, it can be concluded that including many additional timesteps after the rendezvous is completed diminishes the importance of the early timesteps.The two day interval was therefore chosen to allow time to complete the maneuver without diminishing the differences between the cases.

Forecasting with spatio-temporal resolution guidances
Once it is established that the addition of spatio-temporal resolution improves the accuracy of guidances as reflected in the mean Lyapunov function, the next step is to repeat the maneuver with forecasted density.An uploaded density field is assumed to be used by the spacecraft.A similar empirical density model, JB2008, is used as a simulated forecast, representing an imperfect but reasonable forecast.The locations of the JB2008 density points are found by propagating the state of the spacecraft in the same manner as before.By interpolating the JB2008-derived density in the same manner as was previously performed, the density field can be found, and a density estimate can be generated for each timestep.
The guidances are then created in the same manner as before, and are then compared to the executed guidance.The mean Lyapunov function is again used to determine the accuracy of each guidance.Again, when using spatio-temporal resolution, the mean Lyapunov function is reduced, demonstrating the benefit of adding spatio-temporal resolution, shown in Table 2.

Adding receding-horizon control
Using the initial density field to create an initial guidance has the advantage of using a single computation, but has the disadvantage of reduced accuracy as the maneuver progresses.This is due to the positive-feedback loop of the density affecting the location, and vice versa.A small error in estimation of the density can easily result in a large error in position.Because the guidance depends on position as much as density, an error in location results in the guidance being poorly suited for the current position.It is therefore beneficial to periodically update the correct relative position of the spacecraft.
To incorporate these periodic corrections into the guidance, receding-horizon control is used.By updating the guidance after the control horizon has elapsed, the guidance is continually recomputed from a relative position that more accurately reflects the relative trajectory of the spacecraft.Density information is taken from JB2008 in the same manner as before; no changes to the interpolation algorithm are required.
Control horizons of 4 and 8 h were used.Shorter control horizons improve the tracking of the guidance and push it closer to closed-loop maneuvering, but reduce the improvement in error from spatio-temporal resolution since there is insufficient time for much deviation from the guidance if the control horizon is too short.This represents a tradeoff between tracking error and computational requirements; using spatio-temporal resolution with a longer control horizon is a way to minimize the tracking error without introducing excessive computational requirements.
For each horizon, the full guidance is calculated for the remainder of the maneuver; the first guidance is the full maneuver duration, and each subsequent guidance is shortened by the length of the horizon.The maneuver is then executed for the duration of the horizon.Since the density forecast is imperfect, the executed trajectory will not perfectly track the guidance.By updating the starting location and restarting the guidance at the correct location, the new guidance will have a reduced error when compared to the single-computation guidance.The mean Lyapunov function is shown for each case in Table 2, with the time evolution shown in Fig. 6.The periodic jumps are due to the correction of the relative position, which resets the error to zero.Only the high activity case is shown for this and all subsequent figures.
These effects can also be seen in Fig. 7; the increased duration of the horizon allows the execution relative trajectory to drift further from the guidance, increasing the mean Lyapunov function.Also observe that the addition of spatio-temporal resolution has a larger benefit with the longer horizon, as well as the cases without receding-horizon control; this is due to the increased proportion of time not near the beginning of the maneuver, when spatio-temporal resolution has less of a benefit.
It is important to note that the addition of spatio-temporal resolution will not always show an improvement.For the low activity case, the atmospheric density is lower than in the high activity case, and so the addition of spatio-temporal resolution has a reduced benefit, since the trajectories do not separate as much as in the high activity case.Additionally, the rationale of spatiotemporal resolution is based on the assumption that the atmospheric density is linear between the interpolation points.This is slightly incorrect, but the improvement from spatio-temporal resolution generally outweighs this inaccuracy.However, for low activity cases, when the horizon is short, this can be reversed, as seen in the low activity case with the 4 h horizon.However, numerically the error is still far smaller than any other case with forecasting, and the mean error with the addition of spatio-temporal resolution is still within ten percent of the minimum error.
The guidance and execution of the relative trajectory is shown in Fig. 8, using the high activity case and an 8 h horizon (this was chosen for clarity).All guidances extend from the beginning of the horizon through the end of the maneuver.The first execution relative trajectory extends from the maneuver epoch to the end of the first horizon, at which point the guidance is recalculated.Note that later guidances are smaller, in that the chaser does not move much relative to the target.Because the spacecraft start closer together, the error order of magnitude is likely to be smaller for the same density forecast error.However, at this point, the forecasted trajectories have begun to separate, increasing the benefit of adding spatio-temporal resolution.Hence, from this point onward, the main contribution to the mean Lyapunov function has already been determined.

General control performance
For almost every case, tracking of the guidance created with interpolated density resulted in the lowest mean Lyapunov function, indicating that the addition of spatio-temporal resolution generally improves the tracking of guidances.Additionally, the mean error is reduced with the addition of receding-horizon control, indicating its benefit for in-orbit use.Both can be seen in Table 2.The minimum mean Lyapunov function in each case is bolded.
When assuming perfect forecasting, the mean Lyapunov function is reduced when adding spatio-temporal resolution.This indicates that the interpolated density estimate had the closest density estimate to the DTM-2013 density, indicating the benefit of spatio-temporal resolution.Note that the mean Lyapunov function is lower for the low activity case than for the high activity case.This is due to the fact that the trajectories in the low activity case do not separate as much as those in the high activity case, meaning that any trajectory that is contained between the maximum and minimum drag cases will not be very far from any trajectory, so the error cannot be as large as in the high activity case.
The results change after the addition of forecasting using JB2008 (see Table 2).The case using spatio-temporal resolution again has the lowest mean Lyapunov function, but the high activity case with spatio-temporal resolution had a lower mean Lyapunov function than that of the low activity case.Because the forecasting is no longer perfect, there is now uncertainty in the atmospheric density, which allows the actual trajectory to deviate from the guidances.This is seen in the fact that the mean Lyapunov functions are both larger than in the case of simulated forecasting.
The addition of receding horizon control sharply reduces the mean Lyapunov function whether or not spatio-temporal resolution is used, and so it must be considered separately from the addition of spatio-temporal resolution.Since receding-horizon control continually updates the guidance, the error will naturally be lower.With receding-horizon control, the error after a horizon is not dependent on the final error state of the previous horizon, and does not tend to accumulate.Additionally, using a shorter control horizon will tend to further reduce the mean Lyapunov function, for the same reason.The error will be set to zero more often, and so will tend to be lower.This leads to lower mean Lyapunov function for cases using receding horizon control, which is reflected in Table 2, hence, the addition of spatio-temporal resolution reduces the deviation of the tracking relative trajectory from the guidance relative trajectory.
For the low activity case, with a 4-h control horizon, the maximum drag case represents the lowest mean Lyapunov function.Because of the short control horizon, this is expected to have a low mean Lyapunov function.However, the forecasting affects the accuracy of the atmospheric density, and so the trajectory generated with spatio-temporal resolution is imperfect.Since all the trajectories had such low mean errors, spatio-temporal resolution did not provide enough reduction in mean Lyapunov function to offset the effect of the uncertainty in the density forecast.The interpolation process is separate from the forecasting process, but both are required to produce the density estimate.If the forecast and interpolation are perfect, the density estimate will be perfect as well.Since the guidance generated is a result of the density estimate, a perfect density estimate will result in a perfect guidance, with a resultant mean error of zero.
However, although atmospheric density generally tends to decrease with increasing altitude, this does not necessarily hold true for any two individual points.This implies that the interpolation will generally, but not always, reduce the density error by accounting for altitude effects.Because the density error is nearly always reduced, the upper bound for the guidance error tends to be similar to that of the previous state of the art, with the least favorable comparison resulting when the atmosphere has an unusual distribution, and the guidance error is already small.

Conclusions
Spacecraft require precise relative maneuvering to maintain formation flight.Thrusters are not always an option for small spacecraft, so differential drag techniques present an alternative, propellant-less maneuvering method.Open-loop use of differential drag techniques require an accurate estimation of atmospheric density to allow for accurate determination of the resulting trajectory.
Atmospheric density is difficult to estimate accurately, as it is determined by solar and geomagnetic activity, the Earth's daynight cycle, and the position of the spacecraft in the atmosphere.Uncertainty in the atmospheric density forecast results in errors in guidance trajectories, hence, improving the atmospheric density forecast will make the guidance trajectories more realistic.
The addition of spatio-temporal resolution has been shown to reduce the density estimation error, compared to using a single forecasted trajectory (the existing state of the art).Spatio-temporal resolution can be achieved by propagating multiple orbits of spacecraft using a density forecast/estimate, varying the ballistic coefficient for each one.The density-location pairs result in the creation of a density field.This can be obtained on the ground prior to the maneuver, and uploaded at an opportune time.
Interpolating the uploaded density field allows the creation of a relative guidance trajectory.The resulting control inputs are then applied to the spacecraft, with the goal being to follow the guidance relative trajectory.By using spatio-temporal resolution with the simulated perfect forecasting, the effect of only adding spatiotemporal resolution is shown first.Next, after adding imperfect forecasting, the benefit of spatio-temporal resolution is still clearly seen.Finally, an algorithm for the implementation of this method onboard a spacecraft is presented, leveraging the computation speed of the ground station while still retaining autonomous control.

Fig. 1 .
Fig.1.LVLH and inertial reference frames and divergence of orbits with varied drag.

Fig. 4 .
Fig.4.Implementing algorithm onboard a spacecraft.Target spacecraft is shaded more lightly than chaser spacecraft.

Table 1
Initial orbital elements.D. Guglielmo et al. / Acta Astronautica 129 (2016) 32-43 ig. 5. Geomagnetic activity throughout January 2005, as seen in the Hourly Kp Index.Data taken from NASA GSFC OMNIWeb.(For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

Table 2
Control performance metrics.