Spacecraft Formation Flying: Relative Orbit Control based on Euler Orbital Elements
Preprint, Inst. Appl. Math., the Russian Academy of Science
( : , )

G.V.M.Gaias, M.R.Lavagna, A.R.Golikov, M.Y.Ovchinnikov
( ..., .., .., ..)

Russian Academy of Science, Keldysh Institute of Applied Mathematics
. ..

Moscow, 2008


This work addresses a nonlinear Lyapunov-like feedback control to establish and maintain a given relative orbit. Traditional feed-back of cartesian position and velocity tracking error vectors results in an accurate, but expensive, control action. Working with orbit elements eventually allows to reduce such cost by accomplishing corrective manoeuvres on particularly convenient conditions. Moreover it allows to define convenient target relative motions. Hence, by pursuing this approach, here errors in Euler orbit elements are fed-back, thus joining to the previously mentioned benefits the intrinsic more exhaustive exploitation of the dynamics that the Euler set carries with him. These elements, in fact, arise from the solution of the Generalized Problem of Two Fixed Centres, hence already take into account the contributions of J2, J3 and partially J4 zonal harmonics. A critical analysis of the design process of the controller is carried out.

( ) . ( ) , . , . , . , , ( ) . , J2, J3 , , J4 . .

1. Introduction


Formation Flying FF arises when a satellite, orbiting with respect to a massive body, depicts a close relative trajectory around a point which moves under the effect of the same attracting centre. The first satellite is referred as deputy whereas the one that plays the role of the reference point is called chief. Depending on each space application, the chief could be a real satellite or just a virtual point. In cluster missions there are more than one deputies.

From a space system design point of view, FF represents an extremely rich and exploitable working environment. It allows to deal with either space-borne distributed instruments or cooperating systems. Moreover, the further degree of freedom embodied by the relative dynamics allows to perform system reconfigurations to meet either scientific objectives or safety recovering. Nevertheless, this increase of versatility and reliability at the system level involves a set of extremely demanding practical issues from an engineering point of view. Such topics spread over different research fields such as the aim of developing hardware and algorithms for measuring the relative state (position and velocity); the problem of efficiently controlling the relative motion and the search of natural bounded relative orbits. Both the two last features are connected to the dynamical aspect of the formation and affect fuel consumption, mission lifetime and collision avoidance.

By referring only to the dynamics of the centre of mass, the first studies regarding the relative motion are related to the docking problem described through the Clohessy-Wilthshire equations in the Hill frame HCW [1]. They were derived for the Keplerian motion, circular reference orbit and small relative radius and they posses an analytical solution. When dealing with proper FF missions, however, such some severe hypothesis set up a model which is far from the real phenomenon. As explained in [2], denying non linear and eccentricity effects turns out into a really expensive control action. Another source of modelling error is committed when considering the Earths gravity potential as the one produced by a point mass. In [3] are defined families of bounded relative orbits that already take into account the effect of the 2nd order zonal harmonic of the gravity potential. They are described in terms of differences in mean orbital elements. In [4] it is performed a comparison between that mean orbital element approach and the traditional cartesian description in the Hill frame. For both of them a control scheme is provided. In [5] and [6] it is highlighted how the physical insight gained by the orbital element description can be exploited to set up cheap control policies for some given applications. Therefore it can be stated that, the more accurate the model, the less the manoeuvres required for orbit maintenance. Nevertheless, as equations can become more complex, designing an efficient controller can become more challenging.

An alternative approach for achieving a light control amount is to track orbits that do not suffer from harsh disturbances: the idea is to establish some natural bounded motion as it would provide a longer lifetime. In [7] and [8] are performed numerical investigations in order to identify convenient initial conditions for the relative motion.

In the work here presented the idea of exploiting the physical insight that orbital elements provide is kept; though using Euler orbital elements instead of the mean Keplerian ones. According to their definition, Euler elements carry with them information till the 3rd and partially 4th order of zonal harmonics of the gravity potential as they arise from the analytical solution of the Intermediary Motion [9]. As a result the equations that describe the Euler elements time variations already include all the main disturbances that effect FF of small and similar satellites in Low Energy Orbits LEO. As a drawback, however, the controller design becomes more complex, due to the structure of such those equations. Here a Lyapunov-like control were errors in Eulerian elements are fed-back is developed. Some devices for handling the coupling effects embedded in the Intermediary description are needed.

In this report first the main issues concerning the Generalized Problem of Two Fixed Centres GP2FC [10] and the consequent Intermediary Motion IM [9] are recalled. Then a comparison between the equations in the orbital elements for the Keplerian and Eulerian cases is carried out. Finally the control problem is formulated together with the gains tuning process supported by the results obtained by the simulations run.


2. Euler Orbital Elements

The GP2FC problem was first formulated in [11]. It consists in studying the motion of a point mass P under the gravity attraction of 2 fixed points P1 and P2 placed on the z axis of a rectangular reference frame at a complex distance aj. The origin is set in the systems centre of mass and to each centre is assigned a complex mass Mj according to:


P in (x, y, z) experiences the following potential:


By writing W in terms of spherical harmonic series it can be possible to represent the gravitational potential field of an arbitrary rigid body by imposing the conditions (3). The first requirement is satisfied by letting the origin of the frame to coincide with the Planets centre of mass. The following two constraints reintroduce the physics of the problem as they state that all the final quantities shall be real.


In order to meet (3) the quantities aj and Mj must be pairs of complex conjugate numbers obeying to the following structure:


By substituting them back into W and equating the resulting expression with the spherical harmonic expansion of the Earths potential, k-constraints are derived:


Though, due to (4), there are only three degrees of freedom available. Hence after having fixed the one concerning the mass, (5) can only be met for k equal to 2 and 3.


As a result, the potential function W, together with the numerical values computed in (6), perfectly reproduces the effect of the Earth gravity potential till the 3rd zonal harmonics. The level of agreement of the next terms depends on the actual planets mass distribution. For the Earth the 4th term is met up to the 70% of its real value. With respect to the whole Earths potential W denies terms of magnitude less than 1e-6; furthermore it does not take into account sectorial and tesseral contributions. The asymmetrical behaviour with respect to the equatorial plane is described by s. By setting this parameter to zero only the J2 effect would be accounted for.

Being the XYZ frame the Geocentric Inertial Equatorial frame, then the dynamics of a point-mass in the most general case is described by (7): W is the potential of the GP2FC, R collects all the other disturbances written as potentials whereas F represents all the remnant types of accelerations acting on the system, i.e. non conservative contributions and control action. In the pure Keplerian motion none of these three inputs appear. It is defined as Intermediary Motion the dynamics gained when only R and F are zero.


The IM turns out to be a Hamiltonian Autonomous system, hence it possesses 6 constants of motion, among which the first integral of the energy a1 and the conservation of the projection of the angular momentum on the Equatorial plane a3. In order to integrate the IM by quadrature, it is necessary to introduce a set of spheroidal-type coordinates related to XYZ according to:


The motion of P remains near the attracting planet if its energy a1 is negative, that implies, as stated by (8), that the orbit never exceeds the region delimitated by the two ellipsoids of xmax and xmin. The asymmetrical behaviour inducted by s affects the maximum excursions in z; whereas w has the meaning of a longitude.

Starting from the set of constants of motion Euler Orbital Elements EE can be defined in such a way to take into account the qualitative considerations mentioned before; to be the closets of possible to their correspondents in the Keplers set and to possess the cleanest dependence with respect to time. In Table 1 the complete EE set is reported together with their main characteristics.

Tab.1. Euler Orbital Elements




Time Dependence

a > 0







The first element a carries an indication of the distance of the spacecraft from the origin; e is related to the separation between the two ellipsoids bounding the region of motion. The third element i is associated to the maximum height that the spacecraft can reach with respect to the equatorial plane and to the conservation of the projection of the angular motion on the same plane. As a result, for i greater than 90 the orbit is retrograde.

The shape of the region of motion is fully determined by these three elements. It shall be emphasized that, due to their definition, a, e and d involve the computation of x1, x2 and h from the 1st integrals of motion. To achieve such an outcome expansion till the 4th order in e (9) and s can be used for defining a set of mapping relations. The final error gained is of magnitude of 1e-9.


The remnant elements w0, W0 and M0 map the initial position of the satellite on the intermediary orbit. They are worked out from the canonical equations, which involve some elliptic integrals of the 1st type in the amplitude variables q and y, that respectively play the role of the true latitude and true anomaly of the satellite on the intermediary orbit. Once known the initial value of the angular EE, the osculating terms are computed by first solving for the true anomaly y through equation (10) and subsequently by using the relations reported in the last column of Table 1.



4. Satellites Equations for the Perturbed Motion

The closed form solution of the IM already takes into account the secular effects due to the first zonal harmonics. As described in Table 1, they affect the last three EE but they do not modify the shape of the region of motion. Nevertheless, when dealing with the dynamics in the most general case all the disturbing contributions of the right side of equations (7) are present. As a consequence, all the EE would experience some variations with time with respect to the osculating nominal values of the IM solution. In order to evaluate such actions an equivalent set of the Gauss Variational Equations GVE for the Keplerian motion can be written by exploiting the variational methodology and by expressing the disturbances accelerations into a convenient local frame.

Besides, in order to avoid the numerical issues that arise for circular or polar chiefs orbits, a transformation of the EE set into a more convenient one, is performed as shown in (11). To the Eulerian eccentricity value it is summed a vanishing and always positive quantity. Then the forthcoming semilatus rectum is employed. For what concerns the Eulerian inclination the same idea is exploited: from s it is subtracted a vanishing always positive quantity. The rest of the working set is constituted by the angular osculating EE.



From now on, this working set of variables would be referred as the state of the system. Though, to have a straightforward comprehension of the orbit shape, chiefs orbit will be given by a Kepler set. Hence at the initial instant of time a transformation between that set and the working one is performed.

The equations that express how the state varies due to the action of the disturbances not already accounted for in the IM problem can be written as follows:


where A represents the nominal IM solution whereas B describes how the rest of the accelerations acting on the system affect each component of the state. Therefore B is the control influence matrix for this set of equations and by analyzing its elements it can be understood when is convenient to correct for each component of the state, i.e. orbital element. To this extent the set of Keplerian elements listed in Table 2 is assumed and it is performed a comparison between the behaviour of each element of B for GVE and equations (12). Trends are evaluated over 2p radians of the true anomaly y; the eccentricity is varied in the range written down in the table.

Tab.2. Keplerian elements set employed for the evaluation of each term of B

a [km]

e [ad.]

I [deg]

W [deg]

w [deg]

M [deg]








Figure 1 and 2 illustrate the plots obtained, dashed lines refer to the minimum value of eccentricity here considered. In some subplots it is useful to mark, by a square on the axis of abscises, the moment when the true latitude q is 90; respectively the north and south poles of the positive orbits. Every picture is called according to the position of the element that represents into the control influence matrix. Unit dimensions are coherent according to such dimensional equations.

For what concerns the Keplerian motion, i.e. Figure 1, it can be noted that B11 and B21 vary proportionally to sin(y); B12 oscillates around a mean positive value depending on 1/r(y); in B51, B52, B61, B62 small values of eccentricity amplify the nominal cos(y) trend. Such an eccentricity effect can be appreciated also in B22, B33, B43, B52, B53 and B62 where r(y) multiplies the oscillating factor. Finally B33 is proportional to cos(q) whereas B43 and B53 vary with r(y)sin(q). As each colon of B introduces the action of an acceleration into one of the local-frame directions, it is easily derived that in order to correct on semi-axes, eccentricity and mean anomaly in-plane manoeuvres must be used. Changes in inclination and longitude of ascending node are performed by out-of-plane actions, more efficiently respectively at the Equator and at poles. The perigee anomaly can be affected by all the accelerations. Further information is carried by the order of magnitude of each component of B.

Figure 2 reports the control influences components for the Euler problem. Due to the intrinsically more complex problem carried out and to the different choice of state variables, the subplots present different behaviours with respect of the previous ones. In particular the second and third equations now regard p and s and on them accelerations in every direction can produce effects. Together with the strengthening due to small eccentricity values already pointed out before, here there are amplifications due to the dependence of the semilatus rectum on e. This effect is clearly visible in B23, B33, B43, B53. As a consequence, in order to correct produce changes on p and s accelerations either in-plane either out-of-plane can be provided. Still the order of magnitude completes the information required.

Fig.1. Kepler problem: control influence matrix B on one orbit. The dashed trends refer to e = 0.001.

Fig.2. Euler problem: control influence matrix B on one orbit. The dashed trends refer to e = 0.001.


5. Relative Orbit Control Scheme


Satellites relative motion can be fruitfully described through differences in the orbital elements associated to the chief and the deputy spacecrafts. Generally this methodology gives a straight physical insight of the shape of the relative trajectory and a relative orbit control scheme can be established on such an approach. To this extent, if EE are employed, equations (12) need to be used to compute how the state varies when control and disturbance accelerations act on the deputy satellite. The reference relative orbit is defined starting from the chiefs EE. As relative control is sought, no further perturbing actions work on the reference spacecraft:


The controller receives as input the error with respect to the reference state and elaborates which accelerations to provide in order to contrast it:


To avoid to insert at every time step errors connected to the mappings between EE and (11), proper initial conditions are given to the reference state, i.e. dashed block in Figure 3, and is also propagated. Due to the fact that an intermediary orbit is assumed as the reference to be tracked, the controller counteracts only the remnant disturbances and errors in the initial conditions.

Fig.3. Control scheme structure: feed-back based on information in terms of EE.


The disturbance actions considered in this study are the remnant zonal harmonics of the Earth gravitational field, till the 6th order. In the actuation block it is taken into account the maximum acceleration provided by the engines, here assumed as 0.01 m/s2.


6. Design of the Controller

According to (14) the dynamics of the error is given by:


By defining V as a positive defined function of de and by taking its time derivatives on the systems trajectories, u and the matrix P can be designed in order to have an asymptotic stable closed-loop error dynamics, according to:


The real accelerations the controller shall provided are computed as follows:


hence the pseudo-inverse of the 6x3 matrix B is involved. The control is composed by a part proportional to the tracking error, one that compensates the disturbances and one that accounts for the discrepancies between dynamics of the reference and the deputy. The first one is the responsible for recovering the initial error and gains, i.e. terms of the diagonal positive defined matrix P, shall be properly tuned. The second term cancels the effects of the disturbances, hence turns out in an oscillating contribution still present after that the intermediary reference orbit is established. The last component, due to the FF description in terms of orbital elements, is small: differential elements of order 1e-3 can give rise to a large formation geometry.

If were the identity matrix, then P could be simply chosen from the time duration of the transient for recovering the initial error. However, both involves some numerical issues and carries with itself the natural couplings intrinsic in the nature of the phenomenon. Figure 4 and 5, respectively depict the element-by-element trends over 2p radians of those two matrixes. The first picture deals with the eccentricitys span of Table 2, whereas the weighting matrix is computed only for e = 0.05. The pseudo-inverse presents some picks when the determinant of the quantity decreases to quite zero, and especially for small eccentricity values. As a consequence the real control commands turn out to be amplified and edgier than the ideal ones. Figure 5 allows to evaluate the coupling effects among the control actions. In order to accomplish it the diagonal terms should be monitored: they should be dominant with respect to the extra-diagonal ones, to be likely the ideal situation represented by the identity matrix. The first two rows, and columns due to symmetry, satisfy this requirement. All the angular terms, on the contrary, do not fulfil it. And, even whenever for almost the whole orbit the diagonal term is greater than the correspondent extra-diagonal ones, they all are of the same order of magnitude anyway. To fully understand this behaviour it is useful to go back to Figure 2. There, each column can be viewed as how the correspondent control component acts, at a fixed true anomaly, on every element of the state. For what concerns the radial and transversal accelerations, respectively first and second columns, it can be noted that the most dangerous couplings arise between pericentre and mean anomaly and between s and p. In this second case, however, the order of magnitude of the different possible control actions fixes what is more convenient, hence compensates the before mentioned coupling. The out-of-plane action generates a dangerous combined effect between pericentre anomaly and longitude of ascending node. Actions on s and W are naturally decoupled.

Taking into account all these considerations the structure of the gain matrix P can be employed to limit some of those issues. In [4], where an approach on orbital elements was also exploited, it was suggested how a time varying gain matrix turns out to be useful to help in supporting the natural dynamics of the system instead of fighting it. Together with such a heritage, here it is pursued the idea of creating some differences in the order of magnitude among some of the tracking error contributions to unbalance the equilibrium in the structure of that would lead to a steady-state error.

Tab.3. Variable gains: structure and values assumed for the incoming simulations.


P constant


Speed N




















In Table 3 the final choice of structure and numerical values of the gains is resumed. The first two gains are constant as their correspondent equations are decoupled. They are mainly corrected by the in-plane accelerations and the magnitude of the gain was derived by trading off the time for recovering the initial error and the maximum acceleration that the engines can provide. All the other gains are composed by a constant term plus one varying along the orbit. The idea of the constant contribution is to switch off some noise information when it is not needed. On the contrary, when the satellites reaches convenient positions corrections on same components of the state are encouraged by the time varying part. Then the amplitude gives the magnitude effect whereas the exponent, i.e. velocity, gives how fast the switching is activated. For what regards s and W, it is

Fig.4. Euler problem: pseudo-inverse on one orbit. The dashed trends refer to e = 0.001.


Fig.5. Euler problem: matrix that weights the control commands evaluated on one orbit for e = 0.05. True anomaly is on the abscisses, whereas state-component weights are on the ordinates. Sectors divide the rows; diagonal terms are highlighted.

simply encouraged the natural tendency of correcting them respectively at Equator and poles. Their amplitude is the same. Regarding the last two equations, again it is exploited some complementary scheduling of the corrections. Nevertheless, due to the structure of B, when acting on one also the other is affected. The only option available is to emphasized some difference in order of magnitude, together with the one performed when the command is switched off by the angle dependency. For all the time varying quantities, the action is not too impulsive: a smooth behaviour is preferred.


7. Simulations and Results

In the simulations here presented, as a benchmark, it is considered the establishment of a relative motion with respect to a chief satellite on an eccentric orbit, as performed in [4]. Usually the eccentricity of the reference satellite involves some issues in modelling and controlling the relative dynamics. All the Keplerian initial elements of the chiefs orbit are reported in Table 4. They are provided in order to give a straighter insight of the problem. According to the control scheme depicted in Figure 3, this information is first translated into the correspondent Eulerian set and subsequently in the working one.

Tab.4. Chief initial conditions: Keplerian orbital elements.

a [km]

e [ad]

i [deg]

W [deg]

w [deg]

M [deg]








The reference orbit is defined from the chief EE set through a difference in orbital elements. Generally the approach of defining the relative motion through differences between such constants of the two satellites is convenient for at least the following reasons. First there is no need to solve the equations of relative motion to derive information about the spacecraft position, as for the dynamics expressed in a local frame. Secondly it doesnt rely on any hypothesis embedded in the modelling scheme, such as small relative radii for linearization purposes or small eccentricity values of the reference orbit like for Hill-modelling extensions. In the introduction paragraph it was pointed out how the choice of the target relative motion is a crucial aspect for real FF missions. Coherently, when possible, invariant orbits are sought as they do not require strong and expensive control for position maintenance during their lifetime. By exploiting Eulerian orbital elements it is possible to define relative motions which are insensitive up to the effects caused by the presence of J3. Such orbits arise from the observation that, for the nominal IM, the elements connected to the shape of the region of motion do not vary with time. Moreover all the coefficients related to the time dependences in the remnant elements are function of exactly those shape-describing elements. Hence, in order to have a frozen-bounded relative orbit, taking into consideration the deformations due to the first zonal harmonics of the gravity potential, it is simply necessary to ask that all the satellites have matching deviations from the nominal behaviour, i.e. same time derivates of Eulerian longitude of ascending node, pericentre anomaly and mean anomaly. Nevertheless, by asking this for all the angular elements simultaneously turns out into a poor freedom in designing the shape of the relative motion. In this case, in fact, only extensions of Leader-Follower configuration, or oscillations in the out of plane direction are achievable together with a displacement into the along-track direction fixed by a difference in mean anomaly. However, if it can be accepted to perform some corrections on the relative orbit, one of those three constraints can be relaxed, hence establishing some circular motion of the deputy satellite with respect to the reference one. It shall be emphasized that, as here EE are employed, the before mentioned further corrections are tiny because they have only to compensate for the disturbances not already accounted for in the IM.

The relative orbit used in the simulations is defined in Table 5 and it is a J2-J3 invariant relative orbit of the type of a General Circular Orbit GCO. The designing parameter assigned consists in a difference in the Eulerian inclination; differences in a and e are consequently derived to meet the constraints of matching time variations in true latitude and longitude of ascending node between the two satellites.

Tab.5. Definition of the relative orbit to be tracked, through differences of EE.

da [m]

de [ad]

di [deg]

dW [deg]

w [deg]

M [deg]








Simulation depicted in Figure 6 refers to the establishment of the target orbit given in Table 5, when the deputy initial error is given in terms of differences in EE as follows:


On the secondary satellite disturbances are acting as well. In the subplots are respectively shown the transients of the tracking errors in a and p, in the left-top corner, of s and of all the remnant angular orbital elements. In the last view it is sketched the trajectory the deputy performs while trying to achieve the target orbit, plotted in the Hill frame centred in the chief satellite. As expected from the structure of the control weight matrix, a and p present an ideal-like behaviour, and quickly recover their initial errors. Also s slowly moves to the requested value. The other angular elements, on the contrary, present a transient phase and then stabilize on a steady-state value different from the required one, i.e. zero. Such a behaviour was predicted when analyzing the couplings among such terms, and, in the Hill frame, it turns out into a relative orbit of the proper shape but displaced from the target one of a quantity composed by the combination of the steady-state errors on pericentre and mean anomaly. Except for the variables shown in the first subplot, for whom gains were set to be constant, it can be easily recognizable how corrections are performed with an impulse-like approach, coherently with the exploitation of the dynamics of the system. The total DV imparted for performing this reconfiguration is 15.73 m/s, see Table 6, hence greater than what reported in [4] for a similar manoeuvre. However, it shall be remarked that, such a performance index depends on the gains design, the consequent thrust profile, the type of relative orbit sought and, finally, the initial mean anomaly value. Some optimization of the control features shall be accomplished focused on each specific application. For what concerns the final goal, however, a finer result is here achieved as osculating elements are used instead of mean ones, and target orbit already takes into account also the J3 effect.

Fig.6. 1st simulation: Lyapunov-like control. Initial condition and gains are summarized respectively in Tables 4, 5 and 3. Initial errors are given by (18).

In order to avoid the steady-state error, the action of the controller can be enriched by adding a contribution proportional to the error with respect to the aimed value of the pericentre anomaly as follows:


The idea is to force some further action on the system, hence spending something more, to impose some artificial ranking on the corrective manoeuvres to be performed. The gain of this extra contribution is fixed trading off the time needed for recovering that component of the tracking error, the maximum acceleration provided by the engines and the subsequent increase of DV required. According to the dynamics, the proportional part (19) shall be introduced in (16) in order to compute the real commands that the actuators shall provide:


The matrix BP is a design parameter and it regulates which control acceleration, i.e. direction, to use to impart such a further action. Finally, the contribution is reinserted in the equation through the control influence matrix B computed at the current state of the deputy. In Figure 7 there are presented the results gained by a simulation in which the new control action is employed. All the rest of the environment is kept as before. What immediately catches ones attention is the different behaviour of the tracking error of s and of the pericentre anomaly, of course.

Fig.7. 2nd simulation: Lyapunov-like control and proportional part on the error in PA. Initial condition and gains are summarized respectively in Tables 4, 5 and 3. Initial errors are given by (18).

Now the out-of-plane action is more effective and the inclination of the deputys orbit changes more, till the last orbit insertion in the targets plane is accomplished. Because of the artificial ranking imposed on the corrections, the tracking errors on every component of the state, with the exception of the mean anomaly go to the expected value. Yet the steady-state error is decoupled, hence according to the physics of the problem this means that, by simply accomplishing the same manoeuvre some time, i.e. final mean anomaly error, before or later, the expected target trajectory is established.

Fig.8. 3rd simulation: all the features employed in the 2nd simulation are kept. However here a correction in the initial time of the manoeuvre is performed.

Fig. 9. Control accelerations profiles needed for the 3rd simulation. They are expressed in the [S, T, B] local frame.


Fig.10. 4th simulation: Lyapunov-like control and proportional part on the error in PA given by only the radial component S. Initial condition and gains are summarized respectively in tables 4, 5 and 3. Initial errors are given by (18). A proper correction of the initial time of the manoeuvre is performed.

Fig.11. Control accelerations profiles needed for the 4th simulation. They are expressed in the [S, T, B] local frame.


Figure 8 presents the result gained when, together with the complete control law, correction in difference of Eulerian mean anomaly on the definition of the target relative orbit is set. Both the amount of the correction and of the final errors gained are reported in Table 6. Figure 9 shows the control accelerations that the actuators shall provide. Due to the further action the DV required increases with respect to the first simulation performed. Despite this, the final result is much more accurate.

So far the action of the proportional part was carried out by all the three components of the control acceleration, and all of them play a role according to Bd. Nevertheless by considering that the proportional contribution was added to compensate the error in pericentre anomaly, which, according to the dynamics of the system is primarily affected by the in plane corrections, it is convenient to demand all that work to the component S of this supplementary control. The other two components B and T would only insert respectively, a strong disturbing coupling action with the out of plane dynamic and a fighting strategy with respect to the approach fixed for the gains, see Table 3. To achieve this, it can be used the matrix BP which compares in (20). Figure 10 presents the result gained when the correction proportional to the error in pericentre anomaly is performed by the radial-in-plane component S. The profile of s reveals how the manoeuvre less involves changes in inclinations, with great benefits on the total DV spent, see Table 6. Figure 11 shows the correspondent control accelerations, hence it is possible to compare the different behaviour of B in the two cases: now it only works for part of the control law. Finally, as a result, the real trajectory covered by the deputy presents a very different profile, coherently with the fact that small changes in the orbital elements reflect into great modifications of the relative dynamics of the spacecrafts, according to the way they are defined.


Tab.6. Final errors gained in terms of EE and correspondent DV spent by the control action performed in the different simulations presented before.

Final Error

1st simulation

2nd simulation

3rd simulation

4th simulation

da [m]





de [ad]





di [deg]





dW [deg]





dw [deg]





dM [deg]





DV [m/s]





Final Remarks

In this work a relative orbit control scheme for formation flying missions is addressed. The study relies on the exploitation of orbital elements, as they provide a clear physical insight in the description of the dynamics of a point with respect to a massive body. By supporting this approach, relative target orbits and relative motion are expressed respectively by differences in Euler orbital elements and in their variations in time. Eulerian orbital elements represent a convenient mathematical tool as they already involve all the most effective phenomena that occur in the formation flight of small satellites in low energy orbits. However, the intrinsic complexity of the problem moves from the modelling level to the design of the controller: the equations that rule the time variation of the Eulerian elements for the perturbed motion are highly nonlinear and intrinsically coupled.

The proposed control law is based on the theory of Lyapunov. Besides gains that vary across the orbit are used. In order to avoid a final steady-state error a further control action is introduced: it consists in a contribution proportional to the error in pericentre anomaly. Due to such a complex design process, for each special application it is suggested to accomplish an optimization study to work out what are the best gains values and starting time of the reconfiguration manoeuvre, taking into account the available thrusters. This means that, the great enlargement of possibilities that such a comprehensive modelling tool offers, both in terms of accuracy and generality of the situations that can be described, is paid by a certain lack of robustness which affects the control scheme. From a higher level point of view, the exploitation of such a rich modelling tool, together with the proposed control scheme, requires a parallel investigation of the management of the operations during the whole formation flying mission: scheduling of different working modes.

Despite of these features, the proposed approach can be employed either for performing reconfigurations of the formation, either for accomplishing fine control phases, i.e. during scientific segments. The first chance is based on the idea of exploiting the natural dynamics of the phenomenon to achieve expensive manoeuvres; the second one gets benefits from the fact that Eulerian elements already take into account zonal effects till J3 and partially J4. For what concerns tight formations, it is suggested to carry out also a check for collision avoidance: as mentioned when analysing the results gained, small changes in the orbital elements can reflect into great modifications of the relative motion.

Finally, beneath this work an introduction on the design of convenient relative orbits is presented. The idea is to match the time variations of some orbital elements, caused by the environment, in the definition of the orbit to be tracked. As a result the control has to balance only the further disturbances not already taken into account in the model. The simulations run deal with the establishment or reconfiguration of a formation into such frozen-bounded motion.


1. W. H. Clohessy and R. S. Wiltshire, Terminal Guidance Systems for Satellite Rendezvous, Journal of the Aerospace Sciences, pp. 653-658, Sept. 1960.

2. S. S, Vaddi, Modelling and Control of Satellite Formations, Ph.D. thesis c/o Texas A&M University, May 2003.

3. H. Schaub and K. T. Alfriend, J2 Invariant Relative Orbits for Formation Flying", International Journal of Celestial Mechanics and Dynamical Astronomy, Vol. 79, 2001, pp. 77-95.

4. H. Schaub, S. R. Vadali, J. L. Junkins and K. T. Alfriend, Spacecraft Formation Flying Control using Mean Orbit Elements, Journal of the Astronautical Sciences Vol. 48, No. 1, Jan.March, 2000, pp. 6987.

5. S. R. Vadali, H. Schaub and K. T. Alfriend, Initial Conditions and Fuel Optimal Control for Formation Flying of Satellites, AIAA Guidance, Navigation, and Control Conf., JulyAug. 1999.

6. S. S. Vaddi, K. T. Alfriend, S. R. Vadali and P. Sengupta, Formation Establishment and Reconfiguration Using Impulsive Control, Journal of Guidance, Control, and Dynamics Vol.28. No. 2, March-April 2005, pp. 262.

7. J.D. Biggs, V.M. Becerra, S.J. Nasuto, V.F. Ruiz and W. Holderbaum, A Search for Invariant Relative Satellite Motion, ARIADNA Study 04/4104.

8. M. Sabatini, D. Izzo and G. Palmerini, Analysis and control of convenient orbital configuration for formation flying missions, 16th AAS/AIAA Space Flight Mechanics Conference, Tampa, Florida January 22-26, 2006.

9. E. P. Aksenov, Theory of the Motion of Artificial Satellites of the Earth, ed. Nauka, Moscow, 1977 (in Russian).

10. V. G. Demin, Motion of an artificial satellite into a noncentral gravitational field, ed. Nauka, Moscow, 1968 (in Russian).

11. E. P. Aksenov, E. A. Grebenikov and V. G. Demin, Generalized Problem of the Two Fixed Centres and its application in the theory of motion of artificial satellites of the Earth, Russian Astronomical Journal, Vol.40, No.2, pp.363, 1963.