Relative State Transition Matrix Using Geometric Approach for Onboard Implementation

ABSTRACT

In the present study relative state transition matrix was obtained. It matches the relative motion of 2 satellites while including the oblate perturbation. The used formulation applies the geometric approach and is in Cartesian frame. The relative state transition matrix uses absolute state transition matrix of individual satellites. Thereon, simplification in computing methods and on-board implementation at controls are explored in a leader/follower coordination method. Numerical experiments illustrate the accuracy when the baseline separation is equal to 2 km and eccentricity is 0.005 and 0.05 for all the inclinations.

KEYWORDS:
State transition matrix; Formation flying; Orbit; Oblate

INTRODUCTION

Formation flying is about harnessing and controlling the relative orbital dynamics. When using multiple observations obtained from a formation of the low Earth orbiting satellites, the combined data is enhanced. The Clohessy-Wiltshire (CW) equation that describes a relative motion is valid for a small relative distance between the satellites compared to the radius and it assumes that the reference orbit is circular besides using a spherically uniform gravitational field at the centre (Alfriend et al. 2010). In practice this is not true for a low Earth orbiting satellite. Maintaining a formation with respect to the circular reference orbit requires prohibitive fuel consumption. It needs frequent orbit corrections to reduce the deviations caused by perturbations. The first significant perturbation is due to the effect of oblate Earth denoted as J₂. When the reference orbit is initially nominally circular, a small eccentricity usually occurs due to the J₂ effect, which also causes a regression in the orbital plane. It may be noted that an advanced scheduling of the satellite payload data collection on the Earth is usually based on a ground system. This mapping system uses an orbit model that includes mean oblate effect. Thus, in general, it is more important to include the J₂ effect in the on-board reference orbit of the satellite.

State transition matrix (STM) plays an important role in the orbit control. It occurs in the system dynamics equation that is used to control the satellite. The STM that is implemented on-board a satellite needs to be simpler, computationally efficient, and yet scaled close to match the accompanying reference on-board orbit model. Here, the reference orbit on each of the satellites is a first-order motion that includes the J₂ effect (Roy 1982). In this paper, it was obtained an approximate computationally simpler relative STM in Cartesian frame to match this relative motion.

There are various formulations of relative STM and each assumes certain orbital characteristics (Alfriend et al. 2010). Here a rectangular Cartesian frame is used since it helps in the direct use of GPS navigation systems data. Incrementing the maneuver velocity is easier to understand with respect to the orbit frame of reference. Cartesian frame has been found to perform satisfactorily under propagation for formation types that are not largely separated along-track (de Brujin et al. 2001).

The first STM is that of the circular CW equation. The assumptions in the CW model, as mentioned earlier, are the limitations in its use. The STM is derived in Cartesian frame from the analytical solution of the CW differential equation (Alfriend et al. 2010). When the eccentricity of the reference orbit is small, a STM can be derived in Cartesian frame, as carried out by Melton (2000). This involves a lengthy expression as a function of time and is a series with terms which are powers of the eccentricity (Melton 2000). When the orbit is more elliptical then the Keplerian one, it is transformed following the approach of Tschauner and Hempel (1965) to yield a closed form relative motion. This STM given in Yamanaka and Ankersen (2002) is an improvement made to Carter (1998). Here the independent variable is the true anomaly. The time-explicit STM is easier to adopt in the control system.

The method of Gim and Alfriend (2003) has solved the problem of STM, having matched the relative motion about an elliptical orbit with J₂ effect. Though the time-explicit expression is analytical and accurate, it is laborious for on-board implementation. The model includes short and long periodic effects and is in terms of non-singular elements. It is important to mention that there is no need to include the long-period effects under J₂ perturbation. In Hamel and de Lafontaine (2007), this is explored, and an alternative STM was presented, being valid for all cases when the reference orbit is not circular. It is based on mean orbit element differences between 2 satellites. When the eccentricity is small, then a lesser complex version is given in Alfriend and Yan (2005). The STM formulated in Yan et al. (2004) is based on the unit-sphere method (Vadali 2002). The approximate STM in Tsuda (2011) is in Cartesian frame and involves a series for on-board implementation.

Here the approach is geometric, that is, it applies a geometric transformation between the relative states of the deputy with respect to the chief and the orbital frame of the chief satellite. This was adopted to describe the relative motion in Balaji and Tatnall (2003) which is similar to the unit-sphere approach (Vadali 2002) and appeared at the same time. The relative motion described in Vadali (2002) as well as in Balaji and Tatnall (2003) includes the J₂ effect, is applicable for eccentric orbits, and is valid for larger relative distances, using Keplerian elements. The present geometric approach is in Cartesian frame. The relative STM obtained here is approximate and explores optimized computation. Navigation plan suggested here is suited for on-board real-time usage of STM in leader/follower coordination (Alfriend et al. 2010). The STM is applicable for both osculating and mean elements of the 2 satellites. It also provides a scope of including secular effects due to higher-order perturbation. Finally, numerical experiments bring out the accuracy of state updated by the STM in terms of varying eccentricity and inclination.

REFERENCE ORBIT AND STATE TRANSITION MATRIX

The orbital motion of the deputy and chief satellite with J₂ effect is discussed next. The absolute STM of the 2 satellites is then derived. This is required in each of the satellites in controlling their reference orbits.

The equation of motion of an individual satellite is (Roy 1982):

(1) $\ddot{\mathbf{R}}=\nabla \mathrm{U}$

where $\ddot{\mathbf{R}}$ denotes the second derivative with respect to time of R = (X,Y,Z).

The disturbing potential is:

(2) $\mathrm{U}\left(\mathrm{R}\mathrm{,}\mathrm{\psi }\mathrm{,}\mathrm{Z}\right)=\left(\mathrm{\mu }/\mathrm{R}\right)+F$

with

(3) $F\mathrm{=}\mathrm{-}\left(\mathrm{\mu }/\mathrm{R}\right){\left({\mathrm{R}}_{\mathrm{e}}/\mathrm{R}\right)}^2{\mathrm{J}}_2\left(1/2\right)\left\{{3\mathrm{s}\mathrm{i}\mathrm{n}}^2\mathrm{\delta }\mathrm{-}1\right\}\mathrm{,}$

where: µ is the Earth's gravitational parameter, equal to 398,600.64 km³/s²; J₂ = 0.0010826; δ is the instantaneous declination; R_e is the Earth's radius, which is 6,378.458 km; R is the magnitude of the satellite position vector, R.

The position vector of the deputy and the chief satellites is denoted by R_d and R_c, respectively, in the Earth Centred Inertial Frame (ECIF). The potential is axi-symmetric about the Z-axis and is independent of the azimuth angle ψ.

F has 2 types of expressions (Roy 1982):

(4) $\mathit{F}={\mathrm{F}}_1\mathrm{or}{\mathrm{F}}_2$

where: F₁ contains only the secular component while F₂ includes the short-period perturbation.

(5.1) ${\mathrm{F}}_1=\left(\mathrm{\mu }/\mathrm{R}\right){\left({\mathrm{R}}_{\mathrm{e}}/\mathrm{R}\right)}^2\left({\mathrm{J}}_2/2\right)\left(1-\left(3/2\right){\sin }^2\left(\mathrm{i}\right)\right)$

(5.2) ${\mathrm{F}}_2=\left(\mathrm{\mu }/\mathrm{R}\right){\left({\mathrm{R}}_{\mathrm{e}}/\mathrm{R}\right)}^2\left({\mathrm{J}}_2/2\right)\left[\left(1-\left(3/2\right){\sin }^2\left(\mathrm{i}\right)\right)+\left(\left(3/2\right){\sin }^2\left(\mathrm{i}\right)\cos \left(2\left(\mathit{f}+\mathrm{\omega }\right)\right)\right)\right]$

where f is the true anomaly.

The Lagrange's planetary equation of motion (Roy 1982) is:

(6) $\begin{array}{c}\mathrm{d}\mathrm{a}/\mathrm{d}\mathrm{t}=\left(2/\mathrm{n}\mathrm{a}\right)\mathrm{\partial }\mathit{F}/\mathrm{\partial }\mathit{M}\hfill \\ \mathrm{d}\mathrm{e}/\mathrm{d}\mathrm{t}=\left(\left(1-{\mathrm{e}}^2\right)/\left({\mathrm{na}}^2\mathrm{e}\right)\right)\mathrm{\partial }\mathit{F}/\mathrm{\partial }\mathit{M}-\left({\left(1-{\mathrm{e}}^2\right)}^{1/2}/\left({\mathrm{na}}^2\mathrm{e}\right)\right)\left(\mathrm{\partial }\mathit{F}/\mathrm{\partial }\mathrm{\omega }\right)\hfill \\ \mathrm{d}\mathit{i}/\mathrm{d}\mathrm{t}=\left(\cos \left(\mathit{i}\right)/\left({\mathrm{na}}^2{\left(1-{\mathrm{e}}^2\right)}^{1/2}\sin \left(\mathit{i}\right)\right)\right)\left(\mathrm{\partial }\mathit{F}/\mathrm{\partial }\mathrm{\omega }\right)-\left(1/\left({\mathrm{na}}^2{\left(1-{\mathrm{e}}^2\right)}^{1/2}\right)\sin \left(\mathit{i}\right)\right)\mathrm{)}\mathrm{\partial }\mathit{F}/\mathrm{\partial }\mathit{\Omega }\hfill \\ \mathrm{d}\Omega /\mathrm{d}\mathrm{t}=\left(1/\left({\mathrm{na}}^2{\left(1-{\mathrm{e}}^2\right)}^{1/2}\sin \left(\mathit{i}\right)\right)\right)\mathrm{\partial }\mathit{F}/\mathrm{\partial }\mathit{i}\hfill \\ \mathrm{d}\mathrm{\omega }/\mathrm{d}\mathrm{t}=-\left(\cos \left(\mathit{i}\right)/\left({\mathrm{na}}^2{\left(1-{\mathrm{e}}^2\right)}^{1/2}\right)\sin \left(\mathit{i}\right)\right)\mathrm{\partial }\mathit{F}/\mathrm{\partial }\mathit{i}+{\left(1-{\mathrm{e}}^2\right)}^{1/2}/\left({\mathrm{na}}^2\mathrm{e}\right)\mathrm{\partial }\mathit{F}/\mathrm{\partial }\mathit{e}\hfill \\ \mathrm{d}\mathrm{M}/\mathrm{d}\mathrm{t}=\mathrm{n}-\left({\left(1-{\mathrm{e}}^2\right)}^{1/2}/\left({\mathrm{na}}^2\mathrm{e}\right)\right)\mathrm{\partial }\mathit{F}/\mathrm{\partial }\mathit{e}-\left(2/\mathrm{n}\mathrm{a}\right)\mathrm{\partial }\mathit{F}/\mathrm{\partial }\mathit{a}\hfill \end{array}$

where:

(7) ${\mathrm{n}}^2{\mathrm{a}}^3=\mathrm{\mu }$

where: a, e, and i are, respectively, the semi-major axis, eccentricity, and inclination; ω, Ω, and M are, respectively, argument of perigee, longitude of the ascending node and mean anomaly.

There are 2 kinds of reference orbits for on-board implementation on both satellites. It is either (a) osculating and obtained by substituting F₂ in Eq. 6 or (b) the mean motion obtained by substituting F₁ in Eq. 6. The latter is the more familiar first-order secular perturbation with a, e, and i remaining constant, as it is at the start epoch. When mapping the orbit elements to the mean ones, only 3 orbit elements are found to exhibit secular growth due to the J₂ effect and are

(8) $\begin{array}{c}\mathrm{\omega }={\mathrm{\omega }}_0+\left(3/2\right){\mathrm{J}}_2{\left({\mathrm{R}}_{\mathrm{e}}/\mathrm{p}\right)}^2\mathrm{ň}\left\{2-\left(5/2\right){\sin }^2\left(\mathit{i}\right)\right\}\Delta \mathrm{t}\hfill \\ \Omega ={\Omega }_0-\left(3/2\right){\mathrm{J}}_2{\left({\mathrm{R}}_{\mathrm{e}}/\mathrm{p}\right)}^2\mathrm{ň}\mathit{cos}\left(\mathit{i}\right)\Delta \mathrm{t}\hfill \\ \mathrm{ň}={\mathrm{n}}_0\left[1+\left(3/2\right){\mathrm{J}}_2{\left(\mathrm{R}\mathrm{e}/\mathrm{p}\right)}^2\left[1-\left(3/2\right){\sin }^2\left(\mathit{i}\right)\mathrm{)}{\left(1-{\mathrm{e}}^2\right)}^{1/2}\right]\right]\hfill \\ \mathrm{M}={\mathrm{M}}_0+\mathrm{ňt}\hfill \end{array}$

where p = a(1 − e²); ň is the perturbed mean motion.

This has been used in the unit-sphere approach (Vadali 2002). It is more meaningful to use reference orbit with periodic terms for close formations to avoid any collision. On the other hand, having the reference to be a mean motion fuel is saved by having impulse maneuver used once or twice appropriately in the orbit besides maximizing data collection duration (Gill et al. 2007).

Let X_j(t) = (X_j,Y_j,Z_j,X¹_j,Y¹_j,Z¹_j) be in orbital state in the ECIF. At the epoch (t₀), the state is X_j(t₀), where j = 1 denotes the state of the chief and j = 2 that of the follower.

The absolute STM formulation follows the method of Markley (1986). Equation 9 describes the potential and is in Cartesian frame when using periodic and secular terms:

(9) $\mathrm{U}=\left(\mathrm{\mu }{/\mathrm{R}}_{\mathrm{j}}\right)-\left(\mathrm{\mu }{/\mathrm{R}}_{\mathrm{j}}\right){\left({\mathrm{R}}_{\mathrm{e}}{/\mathrm{R}}_{\mathrm{j}}\right)}^2{\mathrm{J}}_2\left(1/2\right)\left(\left({3\mathrm{Z}}_{\mathrm{j}}^2{/\mathrm{R}}_{\mathrm{j}}^2\right)-1\right)$

where: R_j = (X_j²+ Y_j²+ Z_j^2)1/2.

On the other hand, the potential for secular effects alone is:

(10) $\mathrm{U}=\left(\mathrm{\mu }{/\mathrm{R}}_{\mathrm{j}}\right)-\left(\mathrm{\mu }{/\mathrm{R}}_{\mathrm{j}}\right){\left({\mathrm{R}}_{\mathrm{e}}{/\mathrm{R}}_{\mathrm{j}}\right)}^2{\mathrm{J}}_2\left(1/2\right)\left\{\left(3/2\right){\sin }^2\left(\mathrm{i}\right)-1\right\}$

The acceleration derived from the potential for the case (a) with secular and short periodic motion is presented as in Chairadia et al. (2012):

(11) $\begin{array}{c}{\mathrm{a}}_{\mathrm{x}}=-\left(\mathrm{\mu }{\mathrm{X}}_{\mathrm{j}}{/\mathrm{R}}_{\mathrm{j}}^3\right)\left[1+\left(3/2\right){\left({\mathrm{R}}_{\mathrm{e}}{/\mathrm{R}}_{\mathrm{j}}\right)}^2{\mathrm{J}}_2\left(1-\left({5\mathrm{Z}}_{\mathrm{j}}^2{/\mathrm{R}}_{\mathrm{j}}^2\right)\right)\right]\hfill \\ {\mathrm{a}}_{\mathrm{y}}=\left({\mathrm{Y}}_{\mathrm{j}}/{\mathrm{X}}_{\mathrm{j}}\right){\mathrm{a}}_{\mathrm{x}}\hfill \\ {\mathrm{a}}_{\mathrm{z}}=-\left(\mathrm{\mu }{\mathrm{Z}}_{\mathrm{j}}{/\mathrm{R}}_{\mathrm{j}}^3\right)\left[1+\left(3/2\right){\left({\mathrm{R}/\mathrm{R}}_{\mathrm{j}}\right)}^2{\mathrm{J}}_2\left(3-\left(5{\mathrm{Z}}_{\mathrm{j}}^2{/\mathrm{R}}_{\mathrm{j}}^2\right)\right)\right]\mathrm{.}\hfill \end{array}$

On the other hand, for (b), when only secular effects are considered, the accelerations can be seen in Ramachandran (2015) as:

(12) $\begin{array}{c}{\mathrm{a}}_{\mathrm{x}}=\left(\mathrm{a}{\mathrm{X}}_{\mathrm{j}}{/\mathrm{R}}_{\mathrm{j}}^5\right)-\mathrm{\mu }{\mathrm{X}}_{\mathrm{j}}{/\mathrm{R}}_{\mathrm{j}}^3\hfill \\ {\mathrm{a}}_{\mathrm{y}}=\left(\mathrm{a}{\mathrm{Y}}_{\mathrm{j}}{/\mathrm{R}}_{\mathrm{j}}^5\right)-\mathrm{\mu }{\mathrm{Y}}_{\mathrm{j}}{/\mathrm{R}}_{\mathrm{j}}^3\hfill \\ {\mathrm{a}}_{\mathrm{z}}=\left(\mathrm{a}{\mathrm{Z}}_{\mathrm{j}}{/\mathrm{R}}_{\mathrm{j}}^5\right)-\mathrm{\mu }{\mathrm{Z}}_{\mathrm{j}}{/\mathrm{R}}_{\mathrm{j}}^3\hfill \end{array}$

where the constant α = −µR_e²J₂ (3/2) (1 − (3/2) sin²(i)).

It can be noted that the inclination i is invariant in the motion described in Eq. 8. These accelerations are then used to derive the partial derivative and obtain the absolute STM matching the 2 types of reference motions: (a) or (b). Details are given in Markley (1986), as well as in Chairadia et al. (2012).

RELATIVE STATE TRANSITION MATRIX

FORMULATION

Next, it is described the relative motion of 2 satellites by the geometric principle. The motion of the deputy satellite is in the chief satellite's local-vertical-local-horizontal coordinate frame (LVLH). The relative position vector r, with respect to the chief, has 3 components: x, y, and z, as illustrated in Fig. 1. Here, (x) is the component in the chief's radial direction, (z) is in the direction of the chief instantaneous angular momentum, and (y) is along the direction that completes the triad (Alfriend et al. 2010). The relative distance is:

Figure 1
Relative motion of deputy satellite in the frame of the chief.

(13.1) $\mathbf{r}={\mathrm{C}}_{0\mathrm{i}}\left[{\mathbf{R}}_{\mathrm{d}}-{\mathbf{R}}_{\mathrm{c}}\right]$

(13.2) ${\mathrm{C}}_{0\mathrm{i}}=\left(\begin{array}{c}{\mathbf{R}}_{\mathrm{c}}/\left|{\mathbf{R}}_{\mathrm{c}}\right|\\ \left({\mathbf{H}}_{\mathbf{c}}\mathrm{x}{\mathbf{R}}_{\mathbf{c}}\right)\left(\left|{\mathbf{H}}_{\mathbf{c}}\mathrm{x}{\mathbf{R}}_{\mathbf{c}}\right|\right)\\ {\mathbf{H}}_{\mathbf{c}}/\left|{\mathbf{H}}_{\mathbf{c}}\right|\end{array}\right)$

The transformation matrix C_oi (3 × 3) relates the LVLH frame of the satellite to the ECIF. The position and velocity vectors are, respectively, R_c = (X_j,Y_j,Z_j) and V_c = (X¹_j,Y¹_j,Z¹_j), where j = 1 denotes the chief satellite while R_d and V_d for j = 2 are used to denote the deputy satellite. The magnitudes are |R_c| and |R_d| for the chief and deputy satellites, respectively. The instantaneous angular momentum vector of the chief is H_c = (R_c × V_c), where (x) is the usual vector cross-product. In Eq. (13.2), the second row is the normalized cross-product of the first and third rows.

Taking the derivative of Eq. 13.1, the relative velocity is:

(14) $\mathbf{v}={\mathrm{C}}_{0\mathrm{i}}\left[{\mathbf{R}}_{\mathrm{d}}^1-{\mathbf{R}}_{\mathrm{c}}^1\right]+{\mathrm{C}}_{0\mathrm{i}}^1\left[{\mathbf{R}}_{\mathrm{d}}\mathrm{-}{\mathbf{R}}_{\mathrm{c}}\right]$

Combining to get

(15) $\begin{array}{c}{\mathrm{x}}_{\mathrm{rel}}=\left[\mathrm{P}\right]\mathrm{S}\hfill \\ {\mathrm{x}}_{\mathrm{rel}}={\left[\mathbf{r}\mathbf{,}\mathbf{v}\right]}^{\mathrm{T}}\hfill \\ \mathbf{r}=\left[\mathrm{x}\mathrm{,}\mathrm{y}\mathrm{,}\mathrm{z}\right]\mathrm{and}\mathbf{v}=\left[{\mathrm{x}}^1{\mathrm{,}\mathrm{y}}^1{\mathrm{,}\mathrm{z}}^1\right]\hfill \\ \left[\mathrm{P}\right]=\left(\begin{array}{cccc}-{\mathrm{C}}_{0\mathrm{i}}& \mathbf{0}& {\mathrm{C}}_{0\mathrm{i}}& \mathbf{0}\\ -{\mathrm{C}}_{0\mathrm{i}}^1& -{\mathrm{C}}_{0\mathrm{i}}& {\mathrm{C}}_{0\mathrm{i}}^1& {\mathrm{C}}_{0\mathrm{i}}\end{array}\right)\hfill \end{array}$

where S = [R_c R_c¹ R_d R_d¹]^T.

Here [P] is a (6 × 12) matrix. In case the mean elements of the state Š are used, then another transformation needs to be carried out as in Gim and Alfriend (2003), in addition to that of using Cartesian transformation.

(16) $\mathbf{S}\mathbf{=}\mathbf{D}\mathbf{Š}$

The transformation matrix D is reversible and has been sufficiently discussed by Gim and Alfriend (2003).

The following discussion uses osculating element. The formulation is analogous while using mean elements after applying the transformation as in Eq. 16. The relative STM Φ_r (t,t₀) follows:

(17) $\Delta {\mathrm{x}}_{\mathrm{rel}}\left(\mathrm{t}\right)={\mathbf{\Phi }}_{\mathbf{r}}\left({\mathrm{t}\mathrm{,}\mathrm{t}}_0\right)\Delta {\mathbf{x}}_{\mathrm{rel}}\left({\mathrm{t}}_0\right)$

with

(18) ${\mathbf{\Phi }}_{\mathbf{r}}\left({\mathrm{t}\mathrm{,}\mathrm{t}}_0\right)=\left[\mathrm{P}\right]\left[\mathrm{Q}\right]\left[\mathrm{T}\right]$

The matrix P is known in Eq. 15. The second term on the right-hand side of Eq. 18 is the matrix Q and it uses the absolute STM derived earlier for both satellites. If the input is the reference state S(t₀), then the update is obtained by using

(19.1) $\mathbf{S}\left(\mathrm{t}\right)=\left[\mathbf{Q}\right]\mathbf{S}\left({\mathrm{t}}_0\right)$

Using Φ_a to individually denote the STM of the chief and deputy satellites, we have

(19.2) $\left[\mathbf{Q}\right]=\left(\begin{array}{cc}{\mathbf{\Phi }}_{\mathbf{a}}^{\mathbf{chief}}\left({\mathrm{t}\mathrm{,}\mathrm{t}}_0\right)\left(6\mathrm{x}6\right)\hfill & \mathbf{0}\left(6\mathrm{x}6\right)\hfill \\ \mathbf{0}\left(6\mathrm{x}6\right)\hfill & {\mathbf{\Phi }}_{\mathbf{a}}^{\mathbf{deputy}}\left({\mathrm{t}\mathrm{,}\mathrm{t}}_0\right)\left(6\mathrm{x}6\right)\hfill \end{array}\right)$

The matrices are the absolute STM obtained either (a) using periodic accelerations or (b) admitting only secular accelerations as mentioned in Eqs. 11 and 12, respectively. In case the relative STM for mean elements is needed, then the transformation for the osculating to mean, as in Eq. 18, is required.

In order to derive the elements of C_oi¹ in P, a general relation is used

(20) $\frac{\mathit{d}}{\mathit{dt}}\left(\frac{\mathit{u}}{\left|\mathbf{u}\right|}\right)=\frac{{\mathit{u}}^{\prime }}{\left|\mathbf{u}\right|}-\frac{\mathit{u}\left(\mathit{u}·{\mathit{u}}^{\prime }\right)}{{\left(\left|\mathbf{u}\right|\right)}^3}$

Substituting the 3 rows of C_oi in place of (u/|u|) in the above equation, it is obtained

(21) ${\mathbf{C}}_{0\mathrm{i}}^1=\left(\begin{array}{c}{\mathbf{R}}_{\mathbf{c}}^{\mathbf{1}}/\left|{\mathbf{R}}_{\mathbf{c}}\right|-{\mathbf{R}}_{\mathbf{c}}\left({\mathbf{R}}_{\mathbf{c}}·{\mathbf{R}}_{\mathbf{c}}^{\mathbf{1}\mathrm{T}}\right)/\left({\left|{\mathbf{R}}_{\mathbf{c}}\right|}^3\right)\\ \left[{\left({\mathbf{H}}_{\mathbf{c}}\mathrm{x}{\mathbf{R}}_{\mathbf{c}}\right)}^1/\left|{\left({\mathbf{H}}_{\mathbf{c}}\mathrm{x}{\mathbf{R}}_{\mathbf{c}}\right)}^1\right|\right]-\left[\left({\mathbf{H}}_{\mathbf{c}}\mathrm{x}{\mathbf{R}}_{\mathbf{c}}\right)\left(\left({\mathbf{H}}_{\mathbf{c}}\mathrm{x}{\mathbf{R}}_{\mathbf{c}}\right)·{\left({\mathbf{H}}_{\mathbf{c}}\mathrm{x}{\mathbf{R}}_{\mathbf{c}}\right)}^{1\mathrm{T}}\right)/\left({\left|{\mathbf{H}}_{\mathbf{c}}\mathrm{x}{\mathbf{R}}_{\mathbf{c}}\right|}^3\right)\right]\\ {\mathbf{H}}_{\mathbf{c}}^{\mathbf{1}}/\left|{\mathbf{H}}_{\mathbf{c}}\right|-{\mathbf{H}}_{\mathbf{c}}\left({\mathbf{H}}_{\mathbf{c}}·{\mathbf{H}}_{\mathbf{c}}^{\mathbf{1}\mathrm{T}}\right)/\left({\left|{\mathbf{H}}_{\mathbf{c}}\right|}^3\right)\end{array}\right)$

Again the second row is computed as a cross product of the first and third rows instead of the algebraic expression.

The matrix [T] in Eq. 18 is the inverse of matrix [P], and its new form is

(22) $\mathbf{T}=\left(1/2\right)\left(\begin{array}{cc}-{\mathrm{C}}_{0\mathrm{i}}^{\mathrm{T}}& 0\\ -\mathrm{N}& -{\mathrm{C}}_{0\mathrm{i}}^{\mathrm{T}}\\ {\mathrm{C}}_{0\mathrm{i}}^{\mathrm{T}}& 0\\ 0& -{\mathrm{C}}_{0\mathrm{i}}^{\mathrm{T}}\end{array}\right)$

where: N = −C_oi^TC_oi^1TC_oi^T.

The relative STM Φ_r (t,t₀) in Eq. 19 is approximate due to truncation in Q (Markley 1986).

COMPUTATIONAL REDUCTION

The computation reduction aspect is discussed next. The absolute STM are anyway needed on-board for individual orbit control. The additions for relative control are the matrices [P] and [T] and subsequently the matrix multiplication. It is now examined how this computation of the relative STM in Eq. 18 can be decomposed for computation efficiency. It may be noted that there are 4 matrices.

(23) $\begin{array}{c}\left[\mathrm{\alpha }\right]={\mathbf{C}}_{0\mathrm{i}}\left(3\mathrm{x}3\right)\hfill \\ \left[\mathrm{\beta }\right]={\mathbf{C}}_{0\mathrm{i}}^1\left(3\mathrm{x}3\right)\hfill \\ \left[\mathrm{\gamma }\right]={\mathbf{\Phi }}_{\mathbf{a}}^{\mathbf{chief}}\left({\mathrm{t}\mathrm{,}\mathrm{t}}_0\right)\left(6\mathrm{x}6\right)\hfill \end{array}$

The relative STM Φ_r (t,t₀) in Eq. 18 can be rewritten as

(24) ${\mathbf{\Phi }}_{\mathbf{r}}\left({\mathrm{t}\mathrm{,}\mathrm{t}}_0\right)=\left(\begin{array}{cccc}-\mathrm{\alpha }& 0& \mathrm{\alpha }& 0\\ -\mathrm{\beta }& -\mathrm{\alpha }& \mathrm{\beta }& \mathrm{\alpha }\end{array}\right)\left(\begin{array}{cc}\mathrm{\gamma }& 0\\ 0& \mathrm{\lambda }\end{array}\right)\left(\begin{array}{cc}-{\mathrm{\alpha }}^{\mathrm{T}}& 0\\ \mathrm{N}& -{\mathrm{\alpha }}^{\mathrm{T}}\\ -{\mathrm{\alpha }}^{\mathrm{T}}& 0\\ 0& -{\mathrm{\alpha }}^{\mathrm{T}}\end{array}\right)$

The matrix multiplication can be further decomposed as:

(25) ${\mathbf{\Phi }}_{\mathbf{r}}\left({{\displaystyle \mathrm{t}\mathrm{,}\mathrm{t}}}_0\right)=\left(\begin{array}{cc}{}^{\left[\mathrm{f}\right]}\left(\begin{array}{c}{\alpha }^{\mathrm{T}}\\ \mathrm{N}\end{array}\right)& {}^{\left[\mathrm{-f}\right]}\left(\begin{array}{c}0\\ {\alpha }^{\mathrm{T}}\end{array}\right)\\ {}^{\left[\mathrm{-g}\right]}\left(\begin{array}{c}{\alpha }^{\mathrm{T}}\\ 0\end{array}\right)& {}^{\left[\mathrm{g}\right]}\left(\begin{array}{c}0\\ {\alpha }^{\mathrm{T}}\end{array}\right)\end{array}\right)$

where: [f] = [−α 0 ] [γ] and [g] = [β α] [λ].

The repetition of matrices and 0 (3 × 3) matrices involved in matrix multiplication helps to optimise the computation. The following navigation architecture further eases the computation.

The chief or target reference orbit is available in the deputy or follower satellite due to inter-satellite communication (Gill et al. 2007). The thruster on the deputy then controls the orbit to maintain the formation in a leader/follower coordination (Alfriend et al. 2010). The chief reference orbit is undisturbed during this operation. Using the states at t₀ and t₁, the elements of the matrix in Eq. 25 can be computed, ahead of that time. This management ensures that the elements of the first 3 rows of the matrix can be obtained beforehand. The orbit of the deputy or main satellite and the matrix [λ] alone varies due to orbit control thrusters and needs to be updated. This computational reduction in real time is the advantage of the present approach. On the other hand, in the case of unit-sphere approach (Vadali 2002), the transformation matrix in Eq. 15 is constructed using both chief and deputy states. Hence, all the elements of the relative STM need to be computed in real time and are task-intensive. The STM in Gim and Alfriend (2003) is more difficult for implementation.

While using mean motion, the secular accelerations given in Eq. 11 are used to construct mean relative STM. This, however, calls for conversion of the measurement transformed to first-order mean elements using Eq. 16 and that further to Cartesian frame. The orbit model in Eq. 8 can be enhanced to accommodate secular effects of other perturbation for near circular orbits. Towards this, the orbit model suggested in Ramachandran (2015) can be adopted. It means that the motion is more realistic to account for secular variation which is along the track between the satellites in formation, and the STM is useful to match the realistic motion of the satellite.

NUMERICAL ANALYSIS

The experiment carried out here intended to estimate the accuracy of the STM given in Eq. 17. The mean motion is selected as input for each equation. From Eq. 8, the initial states at times t₀ and t₁ are achieved, and the matrix [Q] is obtained. Using the known relative position at time t₁ as input and the approximate STM in Eq. 17, the update of the relative state at time t₂ is obtained. This appears on the left-hand side of Eq. 17. The time step selected is 1 s. Subsequently, it is obtained the difference between this output ∆x_rel (t₂) and the actual relative state at that instant computed using Eqs. 13 and 14. The error in position and velocity in the LVLH frame is related to the reference orbit. This error is gauged against variations in inclination and eccentricity. It is known that the J₂ effect brings an out-of-plane (z) effect which is absent in CW equations. In the first experiment, the orbital eccentricity is 0.05, and the semi-major axis a is 7,000 km. The deputy satellite is at a distance of 2 km along the track. The plot in Fig. 2 shows the error while varying inclination from 0 up to 90 degrees when retaining the eccentricity as 0.05. The numerical accuracy of the methods from Gim and Alfriend (2003) and Yan et al. (2004) shall be no different as they are mathematically similar and differ in formulation. Also, the comparison of propagation error will be limited by the frame of reference. A case of higher eccentricity when e = 0.1 is carried out. The error in the position and velocity was observed to be greater than that shown in Fig. 2, though the variation with inclination is similar. The errors in position and velocity are −0.68 m and −0.41 m/s, respectively. This increase in the error is as expected for higher values of e and due to the approximation made in Q of Eq. 19.2.

Figure 2
Error variation with respect to inclination when e = 0.05.

It is seen that x and x¹ (see Eqs. 13 and 14) have more perceivable error than (y, y¹, z, z¹). The position error is in the top plot while the velocity error is that at the bottom. The (y, z) values agree to 4 decimal places while (y¹, z¹) ones agree beyond 5 decimals. This implies that the update error along these axes is negligible.

To understand the effect of eccentricity, the experiment is repeated by changing the eccentricity to 0.005. The units now are in cm and cm/s - yet again the error is large in (x, x¹). The error in the other axes agrees in 2 decimal places in position and 3 in velocity. When the distance of the satellites is 10 km, then the error gets proportionally scaled to that in Fig. 3. This is due to the approximation in the absolute STM.

Figure 3
Error variation with respect to inclination when e = 0.005.

The next test is carried out by varying the eccentricity from 0.0006 to 0.06. The inclination is fixed and equal to 20 degrees. In Fig. 4, the error along (x, x¹) is greater again and increases with eccentricity. It can be seen in Markley (1986) that the error in absolute STM depends upon the position in the non-circular orbit. In this example, the argument of perigee is retained while the eccentricity varies.

Figure 4
Error variation with respect to eccentricity when inclination is 20 degrees.

As a summary, it is clear that the error, using Eqs. 13 and 14, is the least along the track and across the track (y, z) - these are very important in the control of the satellite. The radial error along (x, x¹) was found more agreeable for smaller eccentricity.

CONCLUSION

Based on the geometric approach, a relative STM was obtained, and it includes the oblate Earth effect. It is simpler for on-board implementation as it is decomposable in terms of the orbital states of the 2 satellites in formation. The navigation plan, as suggested in this paper, reduces the real-time computation when applying the proposed relative STM for orbit control purposes. A numerical experiment was carried out to measure the accuracy of the proposed method for all inclinations. It shows that the approach is more accurate when the distance between the satellites is below 2 km and when the eccentricity is 0.005. For larger baselines and eccentricity, the accuracy of the relative STM is seen limited due to the first-order approximation in the absolute STM.

REFERENCES

• Alfriend KT, Vadali SR, Gurfil P, How JP, Breger LS (2010) Spacecraft formation flying. Oxford: Butterworth-Heinemann.
• Alfriend KT, Yan H (2005) Evaluation and comparison of relative motion theories. J Guid Contr Dynam 28(2):254-261. doi: 10.2514/1.6691
  » https://doi.org/10.2514/1.6691
• Balaji SK, Tatnall A (2003) System design issues of formation-flying spacecraft. Proceedings of IEEE Aerospace Conference; Big Sky, USA.
• Carter TE (1998) State transition matrices for terminal rendezvous studies: brief survey and new example. J Guid Contr Dynam 21(1):148-155. doi: 10.2514/2.4211
  » https://doi.org/10.2514/2.4211
• Chairadia APM, Kuga HK, Prado AFBA (2012) Comparison between two methods to calculate the transition matrix of orbital motion. Mathematical of Problems in Engineering 2012(2012):Article ID 768973. doi: 10.1155/2012/768973
  » https://doi.org/10.1155/2012/768973
• De Brujin F, Gill E, How J (2001) Comparative analysis of Cartesian and curvilinear Clohessy-Wiltshire equations. Journal of Aerospace Engineering, Sciences and Applications 3(2):1-15.
• Gill E, Montenbruck O, D'Amico S (2007) Autonomous formation flying for PRISMA Mission. J Spacecraft Rockets 44(3):671-681. doi: 10.2514/1.23015
  » https://doi.org/10.2514/1.23015
• Gim DW, Alfriend KT (2003) State transition matrix of relative motion for perturbed noncircular reference orbit. J Guid Contr Dynam 26(6):956-971. doi: 10.2514/2.6924
  » https://doi.org/10.2514/2.6924
• Hamel J, de Lafontaine J (2007) Linearised dynamics of formation flying spacecraft on a J2 - perturbed elliprical orbit. J Guid Contr Dynam 30(6):1649-1658. doi: 10.2514/1.29438
  » https://doi.org/10.2514/1.29438
• Markley FL (1986) Approximate Cartesian state transition matrix. J Astronaut Sci 34(2):161-169.
• Melton RG (2000) Time-explicit representation of relative motion between elliptical orbits. J Guid Contr Dynam 23(4):604-610. doi: 10.2514/2.4605
  » https://doi.org/10.2514/2.4605
• Ramachandran MP (2015) Approximate state transition matrix and secular orbit model. International Journal of Aerospace Engineering 2015(2015):Article ID 475742. doi: 10.1155/2015/475742
  » https://doi.org/10.1155/2015/475742
• Roy AE (1982) Orbital motion. Bristol: Adam Hilger Ltd.
• Tschauner J, Hempel P (1965) Rendezvous with a targets in an elliptical orbit.Astronautica Acta 11(2):104-109.
• Tsuda Y (2011) State transition matrix approximation with geometry preservation for general perturbed orbits. Acta Astronautica 68(7-8):1051-1061. doi: 10.1016/j.actaastro.2010.09.016
  » https://doi.org/10.1016/j.actaastro.2010.09.016
• Vadali SR (2002) An analytical solution for relative motion of satellites. Proceedings of the 5th dynamics and control of systems and structures in space conference; Cranfield, UK.
• Yamanaka K, Ankersen F (2002) New state transition matrix for relative motion on an arbitrary elliptical orbit. J Guid Contr Dynam 25(1):60-66. doi: 10.2514/2.4875
  » https://doi.org/10.2514/2.4875
• Yan H, Sengupta P, Vadali SR, Alfriend KT (2004) Development of a state transition matrix for relative motion using the Unit-Sphere Approach. Proceedings of the AAS/AIAA Space Flight Mechanics Conference; Maui, USA.

Publication Dates

• Publication in this collection
  20 July 2017
• Date of issue
  Jul-Sep 2017

History

• Received
  12 Aug 2016
• Accepted
  24 Jan 2017