Small-satellite constellations promise increased ground coverage, higher re-visit rates, and improved sensing resolution. However, a key enabling technology for these constellations is effective, autonomous, relative navigation and control. Because of the size, weight, and power constraints of small satellites, non-impulsive control methods using low-thrust propulsion systems [1] and manipulation of perturbation forces through techniques such as differential drag [2] and solar-sails [3] are gaining traction. A large number of models for the relative motion between spacecraft exist; however, these models aren't well-equipped for the constant low acceleration produced by non-impulsive control systems. In particular, the low accelerations result in long time horizons 978-1-6654-9032-0/23/$31.00 ©2023 IEEE over which the control action must be planned and executed. When formulated using Earth-centered inertial coordinates, the resulting trajectory optimization problems require tens to hundreds of thousands of timesteps-much too large to perform autonomously on an embedded flight computer. In contrast, orbital-element-based models are not generalizable as they must be developed to handle the specific control inputs and perturbations a spacecraft encounters.
Recently, the size of nonlinear trajectory optimization problems for long-horizon orbital maneuvers has been reduced by transforming the orbital dynamics using the Kustaanheimo-Stiefel (KS) transformation [4]. The KS transformation lifts the three inertial position coordinates of the spacecraft into a four-dimensional representation in which the unperturbed Keplerian dynamics become linear and time invariant (LTI). We refer to the four-dimensional KS-lifted position coordinates as the "KS space."
We apply the KS transformation to relative orbital maneuvers between spacecraft in low-Earth orbit. We modify the KStransformed dynamics to include perturbation forces due to non-spherical gravity and low-thrust control inputs. These modifications result in nonlinear dynamics; however, since only the perturbation terms are nonlinear, the KS dynamics are better approximated by linearization than other relative dynamics formulations. This "near linearity" allows for significantly longer step sizes during numerical integration. We solve the relative orbital maneuver problem by linearizing the perturbed KS dynamics with respect to a reference orbit in KS space.
• A novel optimization-based method for smoothly transforming Cartesian states into lifted KS states • A KS formulation of relative orbital dynamics that includes J2 perturbations and low-thrust control inputs • Accuracy comparisons between the KS-based relative dynamics and several existing state-of-the-art relative orbital dynamics models • A convex-optimization formulation of the orbital rendezvous problem using our KS-based relative-orbit dynamics • An example computation of an optimal rendezvous trajectory for a small spacecraft in low-Earth orbit with very low thrust capability
The paper proceeds as follows: We describe related work on relative-orbit models and previous applications of the KS transform to orbital maneuvers in section 2. In section 3 we provide a description of Cartesian and KS transformed orbit dynamics. Section 4 describes our method for transforming smooth trajectories from Cartesian to KS coordinates. In section 5 we derive a relative orbital dynamics model using the KS transform, and in section 6 we compare this model with other relative-orbit models found in the literature by computing the trajectory prediction error versus a numerically integrated ground truth. In section 7 we incorporate the KS relative-orbit model into a convex trajectory optimization formulation, and solve a low-thrust orbital rendezvous problem. We conclude in section 8 by summarizing our results and suggesting future research directions.
There is an extensive literature on relative-orbit models. Sullivan, Grimberg, and D'Amico [5] provide a survey of these models and perform extensive comparisons between them. In section 6, we compare our KS relative-orbital model with the Clohessy-Wiltshire (CW) [6]; Yamanaka-Ankersen (YA) [7]; and Koenig, Guffanti, and D'Amico (KGD) [8] relative-orbit models. These models are developed by linearizing and integrating either the nonlinear Cartesian equations of motion or the Gauss Variational Equations for the orbital elements.
The CW relative-orbit model has been used extensively since the 1960s. It assumes a Keplerian circular reference orbit, is linear-time-invariant, and is parameterized by time. The CW model has been developed for both Cartesian and curvilinear relative coordinate frames [9]; in our comparisons we use the Cartesian coordinates.
The YA relative-orbit model extends the CW model to Keplerian eccentric orbits. It is parameterized by the true anomaly and uses a normalized Cartesian relative state representation. It is considered the state of the art Cartesian relative state representation for arbitrary Keplerian orbits [5]. In the circular case, the YA model reduces to the CW model.
The KGD model uses relative orbital elements and reflects the state of the art in state transition matrices for perturbed elliptical orbits. It provides a significant increase in accuracy over the CW and YA models and has similar or better accuracy to other models in the literature [10], [5].
The Kustaanheimo-Stiefel transformation was originally introduced as a method for regularizing the numerical integration of perturbed two-body motion [11]. It extends the planar Levi-Civita transformation [12] to three dimensions, and provides exact linear-time-invariant equations of motion for unperturbed Keplerian orbits in three dimensions. To our knowledge, the first work applying the KS transform to the relative motion between spacecraft is by Eldin, who studied the KS transform in the context of unconstrained planar rendezvous maneuvers [13]. Thorne and Hall [14] use the planar KS transform to develop analytic expressions for minimum-time continuous-thrust orbit transfers. Hernandez and Akella [15] use the Levi-Civita transformation to find a Lyapunov control policy for finite-thrust orbital rendezvous from arbitrary orbital positions, illustrating the power of working in the Levi-Civita or (more generally) the KS coordinates. Perturbation forces were not considered in these previous works.
Recently, Tracy and Manchester [4] used the KS dynamics and nonlinear trajectory optimization to perform low-thrust transfers from a geostationary transfer orbit (GTO) to a geostationary orbit (GEO). The difference between the approach we present here and the approach in [4] is that we linearize the relative KS dynamics in the presence of perturbations, and perform convex optimization to compute rendezvous maneuvers between multiple spacecraft. Liu and Lu [16] approach the satellite-rendezvous problem using successive convexification methods to approximate the nonlinear relative dynamics and to satisfy safety constraints. In contrast, the linear KS dynamics allow us to solve a single convex optimization problem.
In inertial Cartesian coordinates, the unperturbed Keplerian dynamics of a satellite orbiting a massive body are
where x ∈ R 3 is the position vector relative to an inertial frame centered on the massive body, r = ∥x∥ 2 , and µ is the standard gravitational parameter of the massive body. In low-Earth orbit, the perturbation of these dynamics is dominated by atmospheric drag and the non-spherical shape of the Earth.
We capture the dominant non-spherical gravitational effects by including the J 2 acceleration [17],
where J 2 is the normalized J 2 spherical-harmonic coefficient for the Earth's gravitational field, R E is the radius of the Earth, and x 1 , x 2 , x 3 are the components of x along the axes of the Earth-centered inertial (ECI) coordinate frame. The J 2 perturbed Cartesian dynamics are then
In this work, we focus on the effects of eccentricity and the J 2 perturbation, so we neglect atmospheric drag. However, the model we present can be readily extended to include drag forces.
We now consider the transformation of eq. ( 3) into the fourdimensional KS space [18]. The transformation is not unique when transforming from Cartesian (R 3 ) to KS space (R 4 ), so we first define the transform from KS space to Cartesian.
Let y ∈ R 4 be the KS variable corresponding to the Cartesian position x ∈ R 3 , and define the matrix
The KS transformation from y to x is then
The fourth row of L(y)y is always zero. The matrix L(y) has some useful properties. In particular,
where I is the identity matrix. It follows that
The KS transformation also introduces a scaled fictitious time s that relates to real time by the inverse of the radius,
We denote variables differentiated with respect to real time with a dot, ẋ = dx/dt, and variables differentiated with respect to the fictitious time with a prime, y ′ = dy/ds.
The KS transformed velocity is
Under the KS transformation, the Keplerian two-body dynamics in eq. ( 1) become
where
is the total energy of the orbit. Because h is constant for unperturbed orbits, eq. ( 10) is a four-dimensional simpleharmonic oscillator. In the KS space, the dynamics of any Keplerian orbit are linear and time invariant.
Arbitrary Cartesian disturbance accelerations d(x, ẋ) ∈ R 3 and control accelerations u ∈ R 3 , can be transformed to the KS dynamics. Under perturbation, eq. ( 10) becomes
With acceleration inputs, the energy h is no longer constant:
To account for the energy dynamics, we define an augmented state:
and write the perturbed dynamics
The disturbances d(x, ẋ) can be written as a function of y and y ′ as follows:
For the remainder of this paper, we let d be the J 2 acceleration in eq. ( 3).
For each Cartesian position, x ∈ R 3 , there is a onedimensional submanifold of R 4 such that any point y on that manifold satisfies the KS transform in eq. ( 5). For this reason, the inverse of eq. ( 5), transforming from Cartesian to KS space, is not unique. A single solution y can be found by algebraically inverting eq. ( 5) and arbitrarily choosing the value of one of the elements of y [18], [4]. Unfortunately, while the transformation between the Cartesian and KS spaces is smooth, this method exhibits singularities in the KS space and the transformation of a smooth trajectory into the KS space will not necessarily be smooth. Additionally, when computing the relative position between two KS states, this non-uniqueness leads to two degrees of freedom in the relative position and there may be a relative KS position with smaller norm.
To ensure that transformed trajectories are smooth, and to compute a relative position of minimum norm, we formulate the inverse KS transform as an optimization problem,
where x is the Cartesian position we wish to convert and ȳ is the KS position we desire y to be close to. The solution, y * , of this optimization problem is the KS position closest to ȳ that satisfies the KS transform. To convert points along a trajectory, we let ȳ be the transform of the previous point. If there is no logical ȳ, we let ȳ = [1, 0, 0, 0] T . We solve (17) efficiently using Newton's method.
Figure 1 shows the difference between our proposed nearest state method and the common method. The lines shown are the trajectory of an orbit with unit amplitude and period. The trajectory was originally computed in Cartesian space and has been transformed to the KS space using the common method of fixing an arbitrary element of the state vector and our nearest-state method. The plot shows the transformed KS coordinates of the trajectory. In this case, the common method arbitrarily assigns y 4 = 0 to invert eq. ( 5). This results in the discontinuity at t = 0.5. Our nearest-state method, which minimizes the difference between each point and the previous one, produces a smooth trajectory.
To define the relative KS dynamics, we let z and u be the state and control of the deputy satellite, and we define the chief state z and control ū. The relative state is δz = z -z, the relative control is δu = u -ū, and the relative dynamics are:
Since the Keplerian dynamics for y ′′ are already linear, the higher-order terms in eq. ( 18) are due to the control inputs, perturbations, and the difference in energy between the deputy and chief orbits. The effect of these are orders of magnitude smaller than the Keplerian dynamics, so it is a very good approximation to drop the higher-order terms. This yields the linear-time-varying relative dynamics
To find the discrete-time linear-time-varying relative dynamics, we numerically integrate the controlled state-transition matrix dynamics,
along a given trajectory z, ū. We then define the discrete-time linear-time-varying state-space system
where A k ∈ R 9×9 contains the first nine rows and columns of Φ(s k+1 , s k ) and B k ∈ R 9×3 contains the first 9 rows and last 3 columns of Φ(s k+1 , s k ). Since Φ is computed with the Jacobian of the J2-perturbed KS dynamics, the LTV dynamics include both periodic and secular effects of the J2 perturbation.
We now compare our linearized KS relative-orbit state transition matrix (eq. ( 20)) with the CW, YA, and KGD linear relative-orbit state transition matrices. Table 1 summarizes the differences between these models. The procedure of section 5 is not unique to the KS dynamics, so we also use it to linearize the nonlinear Cartesian dynamics by substituting eq. ( 3) for f and x for z. We refer to the resulting model as the "LIN" model. In [8], the KGD model is developed for three different relative orbital element (ROE) states, the singular ROEs, quasi-nonsingualar ROEs, and the nonsingular ROEs; the nonsingular ROEs are the most general of the representations, so we use them for our comparisons.
Figure 2 shows the root-mean-square (RMS) position error measured along trajectories propagated for one orbital period. The ground-truth orbits are a numerical integration of the J2perturbed nonlinear equations of motion in eq. ( 3) using a high-accuracy adaptive Runge-Kutta method [19], [20]. The deputy initial conditions are perturbed with a range of offsets in mean anomaly, inclination, eccentricity, and semi-major axis while the other orbital parameters are held constant. We use the same reference orbit as in scenario 1 of Sullivan, et al. [5] for comparison with the relative-orbit models they present. Table 2 shows the reference orbit initial conditions used for the mean anomaly, inclination, and semi-major axis variation experiments.
As in [5], to compare performance on eccentric reference orbits, both the reference and deputy orbits are initialized with the same variation of eccentricity. All other initial orbital elements for the reference orbit match Table 2. The deputy orbit is offset from the reference orbit by 0.001 degrees in both mean anomaly and inclination. This corresponds to a distance of approximately 125 meters.
Figure 2 shows that our KS model exhibits significantly less propagation error than any of the other relative-orbit models.
Since the reference orbit is circular for the mean anomaly, inclination, and semi-major axis plots, the CW, YA, and LIN models perform identically. On the eccentricity plot, the CW model exhibits higher error than the YA model, which again matches the LIN model. On the mean anomaly and inclination models, the KGD model exhibits higher error than the KS model within the small angles shown on the plot. The astute reader will notice that extrapolating the KS and KGD mean anomaly and inclination plots, the KS error does grow faster than the KGD error. While not shown, the difference between the KS and KGD error at large mean anomaly and inclination separation angles remains within an order of magnitude of each other. This should not detract from the excellent small-angle performance of the KS model, since relative-orbit models are most commonly used at deviations of less than 10 degrees in mean anomaly or inclination. For the semi-major axis variations, the KGD model exhibits higher error than any of the other models. As an additional note, the extensive use of orbital elements in the YA and KGD models leads to numerous degeneracies and singularities, significantly complicating their practical use. In contrast, the Cartesian and KS dynamics are globally smooth and well behaved.
The LTV relative dynamics given by eq. ( 20) allow us to construct a convex trajectory optimization formulation of the orbital rendezvous problem. With the discretized dynamics in eq. ( 21), trajectories of length N can be computed by solving a convex optimization problem over δz 1:N , and δu 1:N -1 :
where J(δz 1:N , δu 1:N -1 ) is a convex cost function, and Z k and U k are convex sets.
The time steps in eq. ( 22) are scaled fictitious KS times.
Once the optimal trajectory is found, the real times at which δu * 1:N -1 should be applied are found by integrating eq. ( 8) along the chaser states,
We demonstrate the convex trajectory optimization with KS dynamics by solving a low-thrust rendezvous maneuver. The orbits and relative states for this scenario are similar to the International Space Station final approach performed by Soyuz and SpaceX Dragon spacecraft. The target and chaser initial orbital elements, as well as the initial RTN state of the chaser with respect to the target, are given in Table 3. To formulate this problem in the context of eq. ( 22), we assume the target spacecraft is not producing thrust, but is experiencing perturbations, and integrate the target states with eq. ( 15) to compute A k and B k . We additionally right-multiply B k by a rotation matrix which maps vectors in the target spacecraft RTN frame to Earth-centered inertial vectors. This allows us to compute the controls in the target RTN frame, which is a typical choice for formation flying problems [21]. The quadratic cost is,
where Q * { 0, R * { 0. To demonstrate the long optimization horizon possible with KS dynamics, we use a maximum thrust acceleration constraint of 20µm/s 2 . This maximum thrust falls in the range of low-thrust, high specific impulse propulsion systems currently available for small satellites [1]. The optimization uses 20 timesteps per orbit and a 100 orbit horizon, resulting in 2000 knot points. A solution is computed once per orbit, and the controls from that solution are applied to the J2-perturbed nonlinear dynamics over the following orbit in a receeding-horizon fashion. We solve these trajectory optimization problems using the convex quadratic program solver OSQP [22]. It takes approximately 6 seconds to integrate the discrete dynamics, set up, and solve this trajectory optimization on a MacBook Pro with an Apple M1 Pro processor.
The results of this maneuver are shown in fig. 3. The topleft plot shows that the position and velocity errors do not converge monotonically, but do converge to zero over time.
The top-right plot shows the thrust control inputs over time.
The thrust constraints are active for much of the first 50 orbits. The bottom two plots show the chaser trajectories on the radial-tangential plane and radial-normal plane.
We have shown that a relative-orbit-dynamics model based on the Kustaanheimo-Stiefel transformation that includes linearized J2 perturbations and control inputs is highly accurate and achieves higher accuracy than other state-of-the-art models in the literature. The KS relative dynamics model provides a linear-time-varying dynamics formulation that can be incorporated into standard estimation and control tools.
Because the KS relative dynamics are very accurate, longhorizon prediction and trajectory-planning problems can be solved.
Our rendezvous demonstration provides one application of the KS relative dynamics. Many other scenarios are possible, including complex maneuvers with differential drag and solar sails. Additionally, safety constraints are an essential consideration in rendezvous or proximity operations problems that we will investigate in future work.
The code used to produce the results in this paper is available at https://github.com/RoboticExplorationLab/ KSRelativeOrbits.