This section presents the main aspects of the developed constellation model. Specifically, the orbital dynamics is a special form of the Hill–Clohessy–Wiltshire equations [
22,
23]. Such formulation, described in detail in [
24], is based on the definition of a peculiar formation reference frame (the Formation Local Orbital Frame, FLOF) and the formation triangle virtual structure [
21]. These, in turn, lead to the derivation of a novel set of equations describing the spacecraft’s orbit and constellation dynamics. The rationale behind our novel set of CW-type equations is the integration, into an unique model, of the orbit and constellation dynamics through the formation triangle concept [
25], defining the so-called Triangle Dynamics (TD) model. In this context, the term “CW-type equations” means equations based on the perturbation of S/C position, in the LVLH frame, with respect to a chief reference frame in the case of a formation of two satellite. Moreover, with stochastic disturbance dynamics, similarly to [
26,
27], it is possible to account for the orbit’s eccentricity and
effects, thus making the state equations quasi-periodic.
2.1. Introduction and FLOF Definition
The TD is in charge of describing the orbital altitude of the two spacecraft, as well as their relative distance. From this, we obtain an accurate and robust model that allows the inter-spacecraft or inter-satellite distance
d to be easily maintained within a prescribed interval centered at a given nominal value
. In this model, the constellation dynamics is derived by the definition of the FLOF reference frame, described in detail below, with the result of achieving the formulation of a special variation of the Hill–Clohessy–Wiltshire (HCW) equations [
23], describing the S/C orbit and constellation behavior [
25].
Let us now review the rationale behind the TD model design. As already mentioned, it is based on the definition of a specific orbital reference frame, i.e., the Formation Local Orbital Frame (FLOF)
and the triangle virtual structure, allowing a common description of the constellation altitude and the inter-satellite distance. This triangle is a virtual structure for which its vertices are the two S/C centers of mass (CoMs),
and
in
Figure 1a, in addition to Earth CoM,
. The FLOF frame is the constellation common reference frame. Hence, the FLOF axes in Equation (
1) are defined by the relative satellites position
and by the mean radius
, whereas
and
. The FLOF frame axes are defined as follows:
Figure 1 shows a sketched representation of the FLOF axes, together with the triangle virtual structure. The first axis
is referred to as tangential (or along-track), the second one
is reffered to as lateral (or cross-track), and the third one
is reffered to as radial. Other than the FLOF and the formation triangle concepts, for the derivation of the TD architecture, common assumptions for Earth gravity sensing missions (such as GOCE) are taken into account. The satellite nominal orbits are assumed to be polar, with an altitude between 300 and
, with the same right ascension of the ascending nodes (in-line constellation). In addition, the satellite-to-satellite distance variations are measured along the satellite-to-satellite line (SSL). SSL is defined as the line connecting the center of mass of the two satellites, and the first FLOF axis
sets its direction. In a low-Earth-orbit, SSL can be materialized by a differential Global Navigation Satellite System (GNSS).
Further assumptions concern the additional control unit tasks. Indeed, both the drag-free control and the attitude and pointing controls are properly working. As a matter of fact, the TD model assumes that the high-frequency (
) forcing accelerations are only due to gravity periodic components that are related to the scientific objective of the gravimetric mission. This can be considered true if the short-term non-gravitational accelerations are canceled by an appropriate control action. In other words, if it is a drag-free control, both linear and angular settings are available. Note that, despite these assumptions (reasonable in this context), the model remains valid even under different conditions than those defined above (as proved in
Section 4.2).
2.2. Constellation Dynamics
Before detailing the TD model derivation, a brief account of spacecraft and constellation dynamics is provided. The model should be built with respect to the chief frame of reference, namely the FLOF frame defined in Equation (
1). Indeed, the orbit and constellation model requires the preliminary definitions of the orbit and constellation perturbations, which are needed for the definition of triangle dynamics. Let us start by considering the standard satellite CoM inertial dynamics with a GNSS range and rate measurements,
and
, and accelerometer measurements
:
where
is the satellite number, and
and
are the model errors, thus including measurement errors and the effect of neglected dynamics. On the other hand,
and
are, respectively, the low-frequency and high-frequency accelerometer errors, while
denotes accelerometer dynamics. In (
2), we define
that accounts for the non-gravitational accelerations acting upon the
k-th S/C, where
and
are the disturbance and thruster command forces, respectively. Furthermore,
m is the spacecraft mass, while
is the gravity tensor to be measured in order to retrieve the scientific gravimetry signal, namely the object of an Earth-gravity monitoring mission. The dynamics in Equation (
2) is equivalent to the constellation mean and differential dynamics (r and
), which involves, starting from (
2) for each S/C, dropping subscript
k and holding the following.
We also define the following.
Let us recall that the drag-free setting aims at reaching, in (
3), the following conditions.
From (
5), it follows that the drag-free constellation time-evolution is the free response of the Equation (
3). This allows the correlation of the interferometric measurements of
with gravity tensor
to be measured to retrieve the scientific observable of interest. As a matter of fact, the condition in (
5) is achieved by forcing the drag-free command
to cancel the real-time prediction of the environment term
obtained from accelerometer wide-band measurements. Notwithstanding, this cancellation occurs for less than some residuals mainly due to the limited bandwidth of the drag-free control law, causality, and high-frequency accelerometer errors.
The effect of such residuals, as stated before, would impinge on the constellation’s (and attitude’s) stability, without the proper position/velocity feedback ensured by a formation control.
Generally speaking, the constellation and orbit perturbations, needed to build the constellation model, can be expressed in two ways: (i) the standard Hill–Clohessy–Wiltshire perturbations [
23], defined as the 3D S/C CoM displacement
from a Keplerian reference orbit
that is either circular or elliptical [
27] or (ii) the S/C radial and angular perturbations approach, which is adopted in this study. Specifically, the radial and angular perturbations consist in the following: (i) the radial perturbation
,
being the local-vertical-local-horizontal (LVLH) radial axis, and (ii) the angular rate
of the actual orbit LVLH frame. It can be proven that
,
has only two degrees-of-freedom, adding up to the radial perturbation
to a total of three degrees of freedom (DoFs) for the Cartesian perturbations. As a matter of fact, a specific combination of Cartesian and angular perturbations can be extended to a two-satellite constellation with the help of the FLOF frame, defined in Equation (
1). As sketched in
Figure 1b, a nominal reference sphere of radius
(i.e., the nominal constellation altitude) can be associated to the constellation CoM r, defining the nominal orbit, while a nominal constellation SSL, with a length of
(i.e., the mean reference inter-satellite distance) and tangent to the reference sphere, can be associated to the relative position
. Indeed, the constellation radius r and the satellite radii (
and
) hold the following:
where the SSL distance is
d, and the triangle’s height is
, which is the component of the constellation radius r along the third FLOF axis
(see
Figure 1). Furthermore,
is the mean radius component along the SSL direction
, defining the radius’s orthogonality with respect to the triangle base vector
(see
Figure 1). As a matter of fact,
is proportional to the satellite radii difference. Indeed, its nominal condition, namely
, is a target of the control action, which tries to align radius r to
, thus making the triangle isosceles. Consequently, three Cartesian perturbations in FLOF coordinates,
,
, and
, can be defined as follows.
The perturbations defined in Equation (
8) will be directly leveraged to build the TD model based on the FLOF frame. Specifically, perturbation
is referred to as ‘differential altitude’ or ‘differential radius’. Indeed, from
Figure 1b, it is possible to observe that
is proportional to the mean radius difference of the satellite orbits. In addition, three DoFs are provided by the 3D nonzero components of the FLOF angular rate vector
. Therefore, we have a total of six DoFs describing the constellation dynamics.
It is worth point out that the three perturbations in (
8), being defined with respect to a reference sphere of radius
, perfectly fit the in-line constellation type and is characterized by a circular reference orbit.
2.3. TD Model Derivation
In this section, the derivation of the continuous-time constellation dynamics model, namely the Triangle Dynamics model, is considered. Specifically, in the following, a brief development from the triangle differential dynamics to the continuous-time final models is provided. Of course, the model developed here is defined with respect to the FLOF frame (see Equation (
1)).
The constellation dynamics expresses the triangle dynamics based on the differential equations of the six DoFs defined above in
Section 2.1: (i) the three Cartesian perturbations, namely the distance, the differential and mean altitudes (
) defined in Equation (
8), and (ii) the three FLOF angular rate perturbations defined by
.
Hence, by starting from the relative satellite position vector
in (
8), which is aligned to the FLOF axis
by construction, as sketched in
Figure 1, it is possible to derive the differential equation of the inter-satellite distance variable,
d. Let us observe that, since
has three DoFs, we obtain three scalar second-order differential equations: (i) One expressing the second-derivative of the distance, namely the distance variation acceleration
, while two further differential equations detail the pitch and yaw angular accelerations of
, namely (ii)
and (iii)
. As a result, the differential radius
kinematic equation in FLOF coordinates reads as follows:
where
is the FLOF angular rate, while
is the angular acceleration. Consequently, by wisely comparing components with respect to the differential radius kinematics in Equation (
10) with the differential dynamic equation in (
3) (second line), we have the following:
which describes the FLOF dynamics of the first three DoFs (
d,
, and
) out of the six needed to build the complete triangle model (
Figure 1b). In (
11), the FLOF coordinates of the differential external acceleration
in (
4) are highlighted. Moreover, the gravity gradient term
can be expressed in FLOF coordinates and through the identified triangle variables (r and
) in terms of the spherical gravity terms only, while ideally arranging the higher-order terms in the external acceleration contribution
:
where
and
. For the sake of completeness, in this context, with spherical gravity, a gravity model defined by a central potential (spherical symmetry) is intended, defined as
, where
is Earth’s gravitational parameter and r is the distance from Earth’s CoM to the orbiting spacecraft’s CoM. Moreover, we can say that the model simplification in (
12) could impinge on the model capability’s in describing the constellation dynamic behavior only in case that distance
d is not a negligible fraction of the mean constellation radius
. Indeed, in that case, where the dimensional coefficient
, which will be introduced later (see (
18)), is greater than, let us say,
, the
gravity term (i.e., the second order harmonics of the gravity gradient expansion) contribution to the gravity gradient should be accounted for to ensure a more reliable control action. Thus, from (
11) and (
12), the following holds:
where
,
is the mean orbital rate, while
and
come from the definitions in (
8).
As a further step, let us detail the second set of equations completing the six constellation degrees of freedom. For this purpose, let us focus on the mean and the differential altitude,
and
, by computing the kinematic equation of the constellation mean radius r and then applying the same component-wise comparison procedure between the computed radial kinematics (
14) and the mean dynamic equation in (
3) (first line). Thus, if the following constellation radial kinematics holds:
from (
14) and (
3), the radial dynamics in FLOF coordinates read as follows.
Equation (
15) describes the FLOF dynamics of the further three DoFs (
,
, and
) needed to build the complete triangle model (
Figure 1b). In (
15), the mean gravity term,
in FLOF coordinates, is still described through the identified constellation triangle variables (r and
) in terms of the spherical gravity terms only, as specified above. Indeed, the higher-order terms are again ideally accounted for in the three components of the external mean acceleration contribution a from (
4). As a result, the second set of differential equations, concerning the triangle DoFs
,
(
Figure 1b), and
, reads as follows:
where
is defined above.
Provided with the two equation sets from (
13) and (
16) that fully describe two-spacecraft constellation dynamics by means of triangle dynamics and the perturbation variables from (
8), we proceed by adopting the perturbation equations modeling strategy. Indeed, since the control is expected to keep the selected constellation variables close to their nominal value, some equilibrium points and trajectories can be defined for the control design purpose. Hence, we seek the equilibrium by setting the constellation variable derivatives in (
13) and (
16) to zero under zero input signals, thus finding the following equilibrium components and making up the equilibrium state vector
.
Accordingly, by means of (
17), we can define the triangle perturbed state variables, namely
, as well as the external non-gravitational acceleration input vector
:
where the perturbed state vector
is split into its sub-vectors
and
, accounting for the length and the rate state variables all in length unit, whereas the adimensional scale factor
is defined as above. Hence, the orbit and constellation dynamics is based on the perturbation differential equations, for which its equations were built by combining the kinematics equations involving the six selected perturbations with the triangle dynamics, as explained above. The final perturbation equations were used for the control design after having been linearised around the equilibrium point defined in (
17). As a result, by rewriting model Equations (
13) and (
16) in the equilibrium conditions through the perturbed state variables and inputs (
18), the following linear time-invariant (LTI) continuous-time state equations hold follows:
for which matrices
,
,
, and
are parameter free. As a design choice, the angular rate components,
,
, and
, are normalized in (
19), through the definition of the state variables
(roll),
, (pitch), and
(yaw) in (
18). By doing so, all the state variables in
are length increments, with the same unit of measurements, i.e., the meter, and the constellation model looks more homogeneous and is easier to interface with the model’s inputs and outputs. Furthermore, from (
19), we can observe how the roll and yaw state variables,
and
, are decoupled from the other state variables. Likewise, decoupling applies also to input signals
and
, which are orthogonal to the formation plane. The state matrix of the model (
19) is of ninth order, and the nine eigenvalues are as follows.
Hence, the three eigenvalues
need to be stabilized, while the other ones are imposed by the nature. As a matter of fact, by including also the angular kinematics of the triangle, further zero eigenvalues would arise. On the other hand, the out-of-plane constellation dynamics, encompassing the (normalized) roll and yaw angular rates,
and
, are decoupled in (
19) and likewise its input signals
and
, and it will be neglected from here. Thus, from the completed model in (
19), let us consider only the three Cartesian perturbations, namely the differential height
, the mean height
, and the distance
in
and their normalized rates, in
, in addition to the normalized perturbation of orbital rate
. Moreover, the input vector is restricted to the four in-plane input variables in
. It is worth noting that the restricted seventh-order model from (
19) can be proven to be observable and controllable when assuming
, which implies
and
. Hence, after some treatment, we have the following:
where the following is the case.
The eigenvalues of the system are as follows.
From the derived LTI model, it follows that only seven states are required to control a formation of two S/C. This can be a benefit from a control design point of view with respect to the use of the HCW model. Indeed, with the HCW equations, six states are used to provide the approximate dynamics of a single spacecraft with respect to an orbiting frame. This results in the implementation of two different HCW models and then in the use of twelve states for controlling a formation of two satellites.