Stationkeeping manoeuvres for geostationary satellites using feedback control techniques

Aerospace Science and Technology 11 (2007) 229–237 www.elsevier.com/locate/aescte

Stationkeeping manoeuvres for geostationary satellites using feedback control techniques Maniobras de mantenimiento en estación para satélites geoestacionarios utilizando técnicas de control feedback Pilar Romero a,∗ , Jose M. Gambi b , Esther Patiño c a Instituto de Astronomía y Geodesia (UCM-CSIC), Facultad de Matemáticas, Universidad Complutense de Madrid, E-28040 Madrid, Spain b Departamento de Matemática Aplicada, Escuela Politécnica Superior, Universidad Carlos III de Madrid, E-28911 Leganés – Madrid, Spain c Departamento de Matemática Aplicada, E.T.S. Arquitectura, Universidad Politécnica de Madrid, E-28040 Madrid, Spain

Received 6 March 2006; received in revised form 3 August 2006; accepted 3 August 2006 Available online 18 September 2006

Abstract In this paper we discuss the determination of the correction manoeuvres required to control geostationary satellites within a predefined space box over a place above the Earth using a feedback control technique. The control algorithm formulated here combines an accurate numerical propagator for the orbit evolution with an iterative process for the calculation of impulsive manoeuvres that gives numerical solutions convergent to the optimal values, both for the impulses to be applied and for the times of their implementations, to minimize fuel consumption. Numerical simulations are presented to illustrate this technique and to compare it with classical stationkeeping. The numerical results for this nonlinear approach show how the magnitude of the impulses can be reduced over 3 percent per year with respect to the classical approach. This means that the orbital life-time can be increased up to 5 months for a 15 years life-time mission. © 2006 Elsevier Masson SAS. All rights reserved. Resumen En este papel se discute la determinación de las maniobras de correción necesarias para el control de satélites geoestacionarios dentro de una ventana espacial predefinida sobre un punto de la Tierra utilizando una técnica de control feedback. El algoritmo de control aquí formulado combina un propagador numérico preciso para la evolución de la órbita con un proceso iterativo para el cálculo de las maniobras impulsivas que dan soluciones numéricas convergentes a valores óptimos, tanto para los impulsos que han de aplicarse como para los instantes de implementación, para minimizar el consumo de fuel. Se presentan simulaciones numéricas para ilustrar esta técnica y compararla con el mantenimiento en estación clásico. Los resultados numéricos de esta aproximación no lineal muestran cómo la magnitud de los impulsos puede reducirse más de un tres por ciento anual con respecto a la aproximación clásica. Esto significa que el tiempo de vida en órbita puede incrementarse hasta 5 meses para una misión prevista de 15 años. © 2006 Elsevier Masson SAS. All rights reserved. Keywords: Geostationary satellites; GEO; Stationkeeping; Optimal manoeuvres Palabras clave: Satélites geostacionarios; GEO; Mantenimiento en estación; Maniobras óptimas

1. Introduction

* Corresponding author. Tel.: +34 913944584.

E-mail address: [email protected] (P. Romero). 1270-9638/$ – see front matter © 2006 Elsevier Masson SAS. All rights reserved. doi:10.1016/j.ast.2006.08.003

A satellite in geostationary orbit is subjected to various forces which tend to move it from its assigned orbital position. This movement is primarily in the North/South and East/West

230

P. Romero et al. / Aerospace Science and Technology 11 (2007) 229–237

directions. Stationkeeping is therefore required and, therefore, implemented by on-board thrusters [1]. Whereas the aim of the (normal) North/South stationkeeping (NSSK) manoeuvres is to control the oscillations in latitude due to the lunisolar perturbations on the inclination of the satellite’s orbit, the (tangential) East/West stationkeeping (EWSK) manoeuvres are performed to control the satellite longitude evolution, which is due to the drift caused by the perturbing terms of the Earth potential and to the nonzero eccentricity because of the solar radiation pressure effects. In this paper the term “classical stationkeeping” is meant as the technique described in the ESA publication SP-1053 [20] and in the CNES publication [13], and that could be implemented in missions that are flying and/or have recently flown, such as those of Hispasat (30◦ W ). This technique uses analytical expressions to determine secular or long term variations in the orbit evolution, as well as linearized equations to compute the correction manoeuvres and the time of their implementations (see, e.g., [1,16,18]). But due to the increasing number of spacecrafts within the geostationary orbit, which often involves the collocation of several satellites in the same orbital window [11,12,15,22], and because also of the new trends to save cost in the use of electric propulsion for stationkeeping [7, 10,23] the control of geostationary satellites seems to require more precise and robust techniques as well as advanced nonlinear control system designs (for a complete review with detail of geosynchronous satellite control see [3]). This means that, as the dead-band limits for the stationkeeping became narrower, the reduction of the margins for modelling uncertainty, seasonal and short terms effects of satellite perturbations (which vary slightly from cycle to cycle of manoeuvres) should be accounted for. Following this line of thought, we can find several recent works. Thus, it can be found in [14] a variation of the classical stationkeeping that takes into account the actual orbit using a high precision orbit propagator, as well as curve fitting techniques to find out the target orbit in the manoeuvre planning with the impulse effects computed by linearizing equations. In [5], instead, the mean orbital elements are considered and a differential correction algorithm is used to compute the manoeuvres as in [4]. The stationkeeping technique of this paper takes into account the best characteristics of these two works with respect to accuracy and, therefore combines an accurate numerical integration to propagate the orbit with an iterative process for the calculation of the manoeuvres. To check this technique a simulation of this computation is carried out and the numerical results are compared with those corresponding to the classical stationkeeping. 2. Classical stationkeeping 2.1. Geostationary orbit evolution A satellite in geostationary orbit (equatorial and circular with an altitude of 35786 km, to keep as far as possible at rest with respect to the rotating Earth) is continuously perturbed by forces due to the Earth’s gravitational field, the solar radi-

Fig. 1. Annual evolution of the inclination vector for different years.

Fig. 2. Mean longitude evolution.

ation pressure and the gravitational attraction of the Sun and the Moon. In the analysis of the orbit evolution the following set of geosynchronous elements, with small values of eccentricity and inclination, are usually considered (see, e.g., [1] and [17]): the semi major axis, a; the two-dimensional inclination vector, i = (ix = i cos Ω, iy = i sin Ω); the eccentricity vector, e = (ex = e cos(Ω + ω), ey = e sin(Ω + ω)) and the mean longitude (from Greenwich), l = Ω + ω + M − ϑG ; i being the orbital inclination with respect to the equatorial plane, Ω the right ascension of the ascending node, e the dimensionless eccentricity, ω the argument of the perigee, M the mean anomaly, and ϑG the Greenwich sidereal time. The main perturbations on these parameters are illustrated in the subsequent figures. Fig. 1 illustrates the annual evolution, due to the luni-solar gravitation, of the mean inclination vector for a 10 years interval starting from 0◦ at the beginning of each year. The main effect is a secular drift (in fact, with a period of 18.61 years) in a direction, Ωsec , varying between 81.◦ 1 and 98.◦ 9, with periodic components superimposed: the six month periodic terms are the main solar perturbations and the two-week periodic waves correspond to lunar perturbations. The evolution of the mean longitude due to the tesseral harmonics of the Earth’s potential is shown in Fig. 2. Because the orbit is geostationary, the margin for the parabolic motion of the mean longitude, ±δl, for a time interval T between two corrections, is given by [21]: ¨ T 2 /16, δl = |l|

(1)

P. Romero et al. / Aerospace Science and Technology 11 (2007) 229–237

Fig. 3. Annual evolution of the centered eccentricity vector.

where l¨ denotes the mean drift rate corresponding to a given longitude ls . According to its specific localization the satellite finally tends naturally to shift east or westward. Fig. 3 shows the approximate circular annual evolution of the mean eccentricity vector due to the solar radiation pressure and, superimposed, the periodic small oscillations due to the lunisolar attraction. As is well known, the nonzero eccentricity induces daily oscillations in longitude. Therefore, if for a specified window band in longitude a margin for the mean periodic perturbation due to the Moon is considered, then the eccentricity vector evolution can be approximated by the following equation which represents that e moves pointing to the Sun in a circle, with center e0 and radius Re that depends on the geometrical and physical characteristics of the satellite, i.e.:    = e0 + Re cos s (t) , e(t) (2) sin s (t) where s (t) is the Sun right ascension at time t. The aim of the eccentricity control is to reduce the radius of this circle to a value, ec , which gives the limit circle for the mean eccentricity vector that avoids the violation of the E/W band, δλ = δl + 2 × ec ,

(3)

231

On the other hand, the burn direction in the East/West stationkeeping is made tangent to the orbit and makes both the drift rate and the eccentricity of the orbit adjust to maintain the satellite within the longitude band ±δλ given in (3). In the case that the solar radiation pressure effect is too large, which occurs most often, two tangential burns, the first manoeuvre being East (positive) or West (negative) according to the direction of the longitude drift, and separated half a sidereal day, have to be implemented to maintain the satellite within the specified limits in longitude. The effects on the corresponding geosynchronous elements from a two tangential burns V1 and V2 are given by the following linearized equations [21]:   2 cos sb , (6) e = (V1 − V2 ) sin sb V   3 l = − (V1 + V2 ) Ω⊕ (t − tb ) − 90◦ , (7) V where sb is the sidereal mean time at the satellite corresponding to the first thrust at the time tb , Ω⊕ is the Earth’s angular velocity, V is the geosynchronous velocity, and e = e + − e − , l = l(t) − l − ,

(8)

where the notations e − and l − depict the eccentricity vector and the mean longitude before the first manoeuvre; e + is the eccentricity vector just after the second manoeuvre and l(t) is the mean longitude for any time t after the second manoeuvre. 2.3. Optimal strategies The challenge in planning the manoeuvres is to design optimal strategies that, satisfying the mission constraints, minimize the propellant consumption. The Secular Mean Line (SML) strategy is usually chosen to plan the NSSK [2,19]. This means that with this strategy the correction is applied in the direction of the secular drift for the long term evolution of the inclination vector , which implies that the target inclination vector i + = iT must satisfy the condition

where δl is given by (1).

arg(i) = −Ωsec .

2.2. Linear equations for the stationkeeping manoeuvres

For the EWSK manoeuvres the optimal direction is obtained by pointing the perigee of the satellite orbit towards the sun, i.e., according to the strategy called Sun Pointing Perigee (SPP) [6, 9]. In practice, to carry out the SPP strategy it is sufficient the requirement that the direction of the perigee (given by e) points to the Sun at some instant during the cycle. In particular, taking into account Eq. (2) it can be seen that, by imposing them to coincide at the middle of the cycle, the target mean eccentricity to achieve, e+ = eT , must satisfy the condition     − cos s (tb + T /2) exT = e0 eyT − sin s (tb + T /2)   cos s (tb ) , (10) + Re sin s (tb )

As it has been said, the inclination vector correction is performed by means of a normal impulse to the orbital plane according to the equation given in [21]:   cos sb , (4) i = −Vn /V sin sb where the angle sb is the sidereal mean time of the satellite at the moment of the burn, and i = i + − i − . (The notations i −

(5) and i +

and after the manoeuvre.)

depict the inclination vector before

(9)

232

P. Romero et al. / Aerospace Science and Technology 11 (2007) 229–237

where e0 is the shifting away from the origin of the center of the eccentricity circle in the Sun direction at the time tb + T /2 obtained by an impulse at the time tb . By imposing in (10) that |eT | = ec , we obtain that e0 must have the value   e0 = Re cos s (tb + T /2) − s (tb )    − ec2 − Re2 sin2 s (tb + T /2) − s (tb ) . (11) It is characteristic of these optimization strategies that they define the direction of the manoeuvres (which, in turn, determines the time of the day for the thrusts) but do not define neither the day nor the size of each manoeuvre. Finally, it must be said that to ensure simplicity in the operations, the stationkeeping strategies generally used are based on repetitive and constant times cycles of recurrent sequences of orbit corrections, and thus, according to a fixed schedule previously selected, the determination of the thrust size to avoid violation of the specified margins in latitude and longitude of the orbital window is accomplished. To compute manoeuvres in the classical stationkeeping with linear equations the mean orbital elements are used. The evolution of these elements is obtained by means of the following linearized Lagrange equations where the perturbing function only contains those terms causing secular and long period perturbations: 2 ∂R1 dl = n − Ω⊕ − , (12) dt na ∂a da 2 ∂R1 = , (13) dt na ∂l where n is the mean motion and R1 is the terrestrial perturbing potential; dex 1 ∂R2 =− 2 , dt na ∂ey dey 1 ∂R2 = 2 , dt na ∂ex

(14) (15)

where R2 is the potential due to the lunisolar attraction, and dix 1 ∂R3 =− 2 , dt na ∂iy diy 1 ∂R3 = 2 , dt na ∂ix

with the precise numerical integration of the differential equations of the satellite motion has also been implemented in order to study the possibility of improving the accuracy and to minimize the fuel consumption. In order to compare the osculating elements σ (t) (obtained by transforming the cartesian positions and velocities previously derived by this numerical integration), with the mean orbital elements used in the classical stationkeeping, an intermediate set of centered elements are considered. These elements, σc (t), are obtained by applying to the osculating elements, first, the numerical filter 1 if t ∈ [tk , tk + Tc ], H (t) = (18) 0 if t ∈ / [tk , tk + Tc ], with Tc = 24 hours, to derive the filtered parameters σf (t) = H (t) ∗ σ (t); and then, by averaging σf (t) in Tc . The centered elements thus obtained, are assigned to the middle of each time interval [tk , tk + Tc ]. Therefore, σc (tk + Tc /2) =

N 1

σf (ti ), N

with N = Thc and ti = tk + ih, h being the integration step (h = 22m 51.s 3 ≡ 0.1 rad). Figs. 4, 5 and 6 show the comparison between the osculating and centered elements and Figs. 7, 8 and 9 show the comparison between the centered and mean elements. In this procedure the effects on the orbital elements, Y , produced by a set of thrusts Vi at times ti when these effects are calculated by means of the numerical integration of the satellite motion equations, give rise to a function, which is denoted by:  = Y , F (X) (20) and when the effects are calculated with the linear equations (4), (6) and (7) we have another function, denoted by f ( x ) = y. Then, to avoid the fact that Eq. (20) cannot be inverted to  that produce the specidetermine the vector of parameters X fied target values for the orbital elements, we use the iterative process   xn ) + Y − F ( xn ) (21) xn+1 = f −1 f (

(16) (17)

where R3 is the potential due to the solar radiation pressure. 3. Stationkeeping using feedback control 3.1. Iterative process for the manoeuvres optimization A high precision orbit propagator that uses the Gauss– Jackson method [8] has here been implemented for the numerical integration of the geostationary satellite motion equations. It includes the Sun/Moon attraction, the solar pressure and the Earth gravity up to the fifth order and fifth degree. Then, an iterative process that combines the standard linear model for the computation of the manoeuvres (summarized in Section 2.2)

(19)

i=1

Fig. 4. Osculating and centered longitudes.

P. Romero et al. / Aerospace Science and Technology 11 (2007) 229–237

Fig. 5. Osculating and centered eccentricity vectors.

Fig. 8. Centered and mean eccentricity vectors.

Fig. 6. Osculating and centered inclination vectors.

Fig. 9. Centered and mean inclination vectors.

233

Fig. 10. Iterative process.

3.2. North/South manoeuvres optimization Fig. 7. Centered and mean longitudes.

in order to solve the equation in (20) by means of a series of  as is illustrated in Fig. 10, [21]. approximations, xn , to X The function f ( x ) = y in (21) is therefore the analytical version of the function F in (20) but now calculated with the linear model cited previously. (Therefore, f has an inverse function, f −1 , that can be calculated easily.)

From this point of view, the problem consists in finding out  = (V , sb ) of Eq. (20) that gives the normal imthe solution X pulse V at a sidereal time sb to change the inclination vector according to the SML strategy given by (9). Using the same notation convention as in (5); denoting by yn = iTn = (ixTn , iyTn ) the set of approximations determined with the linear model (4), and denoting by Yn = iTNn = (ixNTn , iyNTn ) the series of results obtained by the numerical

234

P. Romero et al. / Aerospace Science and Technology 11 (2007) 229–237

integration, then, by applying (21), we have got the following recurrent relations: a) Vn2 2 Vn+1 = V + iT2 (sb ) + iTNn (sbn ) V2   + + − 2 ix+T (sb )ixNT (sbn ) + iy+T (sb )iyNT (sbn )   2Vn  ixTn ixNT+ (sbn ) − ix+T (sb ) + V iTn 1   2 (22) + iyTn iyNT+ (sbn ) − iyT (sb ) to obtain the (n + 1)-term of the series Vn in order to compute the impulse, and b)  sbn+1 = arctan ix2Tn + iy2Tn sin sbn  − iy+T (sb ) + iyNT+ (sbn )  × ix2Tn + iy2Tn cos sbn −1 (23) − ix+T (sb ) + ixNT+ (sbn ) to compute the time for the manoeuvres. 3.3. East/West manoeuvres optimization = This problem consists now in finding out the solution X (V1 , V2 , sb ) of Eq. (20) that gives the tangential impulses, V1 and V2 and sb (the sidereal time corresponding to the first impulse, V1 ) to modify both the longitude and the eccentricity. The correction for the eccentricity vector is made here according to the SSP strategy given in (10) to minimize |V1 | + |V2 | for each ec determined according to Eq. (3) to avoid the violation of the E/W band. However, to correct the longitude, two different criterions, described below, have been considered. 3.3.1. Target value In the first case, we look for the resulting true longitude, after applying two tangential burns, to reach the limit of the deadband in longitude at the end of the cycle, i.e., λT (tb + T ) = ls − δl (see (1)). Using the same notation convention as in (5), and denoting now by yn = (λTn , exTn , eyTn ) the series of approximation determined with the linear model (6) and (7) and by N N Yn = (λN Tn , extT n , eyTn ) the series of results obtained by the numerical integration, then, by applying (21) and by introducing the auxiliary parameters

(24)

which represent the change in velocity to correct the true longitude and the eccentricity vector, respectively, we have the following recurrent relations: a) N V λTn (tbn + T ) − λT (tb + T ) Un+1 = Un + · 3 Ω⊕ T − 90◦

to obtain the (n + 1)-term of series Un and Wn , and b)  ex2n + ey2n sin sbn sbn+1 = arctan  + eyT (tb + T /2) − eyNT (tbn + T /2)  × ex2n + ey2n cos sbn + exT (tb + T /2) −1 (27) − exNT (tbn + T /2) to compute the time for the first manoeuvre. 3.3.2. Symmetric criterion In the second case, as a result of the manoeuvres, the true longitude has to be symmetric about the nominal longitude with its maximum and minimum in the middle and at the end of the cycle, respectively, i.e., λ(t) − ls = ls −

max

tb ttb +T

min

tb ttb +T

λ(t).

(28)

The criterion corresponding to the eccentricity vector is the same used in the previous section and, therefore, the expressions for Wn+1 and sbn+1 are again given by (26) and (27). Here we consider only the effect of the manoeuvres on the true longitude, Un , and a new parameter, Cn , is needed to represent the symmetric criterion. Cn is given by: Cn =

max

tbn ttbn +T

= ls −

λn (t) − ls

min

tbn ttbn +T

λn (t).

(29)

Here, the vectors involved in the iterative process are  = (U, C), yn = (λMT , λmT ), and Yn = defined by X n n N (λMT , λN mTn ), with n

λMT =

max

λ(t) − λ− b

(30)

min

λ(t) − λ− b.

(31)

tb ttb +T

and λmT =

tb ttb +T

Finally, considering the auxiliary parameters

Un = V1n + V2n , Wn = V1n − V2n ,



4 2 2 W + eT2 (tb + T /2) + eTNn (tbn + T /2) V2 n   4 Wn  exn exT (tb + T /2) − exNT (tbn + T /2) + V en   + eyn eyT (tb + T /2) − eyNT (tbn + T /2)  − 2 exT (tb + T /2)exNT (tbn + T /2) 1  2 N + eyT (tb + T /2)eyT (tbn + T /2) (26)

V Wn+1 = 2

(25)

Umax =

V maxtbn ttbn +T λN n (t) − ls − Cn · 3 Ω⊕ (t − tbn ) − 90◦

(32)

Umin =

V mintbn ttbn +T λN n (t) − ls + Cn · , 3 Ω⊕ (t − tbn ) − 90◦

(33)

it is easy to verify that (see (7)): a) in the middle of the cycle, t = tbn + T2 ,

P. Romero et al. / Aerospace Science and Technology 11 (2007) 229–237

3 · (Ω⊕ T /2 − 90◦ ) V = λN n (tbn + T /2) − ls − Cn ;

235

Umax ·

(34)

and b) that at the end of the cycle, t = tbn + T , 3 · (Ω⊕ T − 90◦ ) = λN n (tbn + T ) − ls + Cn . V Now, introducing the function

Umin ·

(35)

 3 Ω⊕ (t − tbn ) − 90◦ , (36) V and, using (21) we have the following recurrent relations: a)

g(t) =

Un+1 = Un +

Umax g(tbn + T /2) + Umin g(tbn + T ) g(tbn + T /2) + g(tbn + T )

(37)

to obtain the (n + 1)-term of the series Un to compute the magnitude of the impulse, and b) Cn+1 = Cn + (Umax − Umin ) g(tbn + T /2)g(tbn + T ) × g(tbn + T /2) + g(tbn + T )

Fig. 11. NSSK with the SML strategy for the mean inclination vector.

(38)

to test that the symmetry criterion is reached. The values Cn can also be used to adjust optimally the evolution of the true longitude and, therefore, to maximize the value of ec . 4. Numerical simulation In order to compare the two different approaches herein described to compute the total cost of the impulses used to control, first, the mean orbital elements by using the linear formulas given in Section 2.2, and second, by using the iterative relations derived in Sections 3.2 and 3.3 applied to the osculating orbital elements, several numerical simulations have been carried out to evaluate the total V required in a year of North/South and East/West manoeuvres in the stationkeeping process. To this end, a satellite located at ls = 30◦ W, with reduced control bands of δλ = ±0◦ .05 for the longitude and δφ = ±0◦ .05 for the latitude, and a radius Re = 5.6 10−4 for the evolution of the mean eccentricity vector have been considered. Furthermore, to ensure simplicity in the operations, weekly manoeuvres cycles of 14 days have been chosen. Figs. 11 and 12 illustrate the NSSK with the SML strategy. The simulations have a total of 25 manoeuvres in a year and the inclination is maintained within a circle of radius ic = 0◦ .05. The classical approach for the mean inclination vector is shown in Fig. 11. A magnitude of 1.616 m/s for each impulse, computed by using formula (4), has been implemented. This corresponds to an annual cost of V = 40.416 m/s. The case for the osculating inclination vector that uses the complete dynamical model is shown in Fig. 12. Here, small waves appear within the control circle and impulses have been computed using the iterative formula (22). With 0 as initial value, the algorithm converges to the target values after 4 iterations, and the total amount of the impulses is V = 40.267 m/s. It is clear that in the NSSK the iterative method does not offer any significant improvement over the classical technique.

Fig. 12. NSSK with the SML strategy for the osculating inclination vector.

Fig. 13. EWSK with the SPP strategy for the mean longitude.

Fig. 13 shows the classical EWSK for a 3 month interval with the SPP strategy using target values at the end of the cycle for the mean longitude. With a mean drift rate of l¨ = −0.000887 deg /day2 for the longitude of ls = 30◦ W, and a margin of 0.◦ 004 for the mean periodic perturbations due to the Moon, according to (1) and (3), a control eccentricity value of ec = 1.6 10−4 has been considered to compute the impulses by using (6) and (7) Fig. 14 illustrates the control of the eccen-

236

P. Romero et al. / Aerospace Science and Technology 11 (2007) 229–237

Fig. 16. EWSK with the SPP strategy for the true osculating longitude and symmetrical criteria.

Fig. 14. EWSK with the SPP strategy for the centered eccentricity vector.

Fig. 15. EWSK with the SPP strategy for the true osculating longitude.

tricity vector for this case. The total V required in a 3 months stationkeeping is now 0.937 m/s. It is noteworthy that the linear prediction for the V cost is quite close to the iterated value V = 0.932 m/s obtained by means of (25) and (26) after 5 iterations if the value of ec = 1.6 10−4 is retained. Now, the EWSK with the SPP strategy using target values at the end of the cycle for the true longitude obtained by numerical integration is illustrated in Fig. 15. From the study of this simulation it results that, using the true longitude, the deadband for the periodic lunar perturbations can be removed so that the eccentricity control value may be increased up to a value of ec = 2.3 10−4 and, as a result, the V cost reduces to V = 0.810 m/s. In this respect, it should be noted that when ec is increased up to this value the mean longitude cannot be kept within the prescribed dead-band. However, as can be seen in Fig. 15, there still remains a deadband in longitude with this procedure. Fig. 16 shows how the true longitude fits into the total control band δλ by using the iterative procedure to compute thrusts that produce symmetric longitudes about a given longitude with the formula (38). The conclusion is that the eccentricity control can even be increased to a value of ec = 3.5 10−4 and, as a result, the V cost obtained after 9 iterations by (26) and (37) reduces to V = 0.601 m/s for a 3 months stationkeeping. On the other hand, the annual stationkeeping of the mean eccentricity vector within a circle of radius ec = 3.5 10−4 , is

Fig. 17. Annual stationkeeping for the centered eccentricity vector with the SPP strategy.

Fig. 18. Annual stationkeeping in longitude-latitude.

shown in Fig. 17. Its annual cost is V = 2.404 m/s. Finally, the annual stationkeeping of Fig. 18 shows how the satellite is maintained with this last strategy within its orbital window with the symmetric criteria. By comparison of both process the result is that the cost for the EWSK manoeuvres in a year reduces an

P. Romero et al. / Aerospace Science and Technology 11 (2007) 229–237

amount that goes from 8% (corresponding to the standard strategy), to 5% of the total annual amount. This means, in terms of life-time, that the orbital life-time could be increased up to 5 months for a 15 years life-time satellite. 5. Conclusions A feedback control technique has been applied to determine the optimal correction manoeuvres for the control of geostationary satellites. The numerical simulations carried out show how in the North/South stationkeeping there is not significant improvement if the impulses are computed by using this iterative procedure with the accurate numerical integration scheme proposed. However, the cost for the East/West stationkeeping manoeuvres according to the Sun Pointing Perigee strategy can be optimized if this technique is used to produce symmetric longitudes about a given longitude. In fact, it may be concluded that the application of this method outperforms the classical stationkeeping (which uses linear equations for the manoeuvres and mean elements derived with Lagrange equations) since it provides a way to extend geostationary mission life-time, which may go up to 5 months. References [1] P. Berlin, The Geostationary Applications Satellite, Cambridge University Press, Cambridge, USA, 1988. [2] J. Chan, Inclination control strategy for Intelsat VIII satellites, in: Proceedings of the 50 International Astronautical Congress, Amsterdam, The Netherlands 1999, paper IAF-99-A.1.07. [3] C. Chao, J. Baker, On the propagation and control of geosynchronous orbits, The Journal of the Astronautical Sciences 31 (1) (1983) 99–115. [4] T.A. Ely, K.C. Howell, East–West stationkeeping of satellite orbits with resonant tesseral harmonics, Acta Astronautica 46 (1) (2000) 1–15. [5] B.P. Emma, H.J. Pernicka, Algorithm for autonomous longitude and eccentricity control for geostationary spacecraft, Journal of Guidance Control and Dynamics 26 (3) (2003) 483–490. [6] D.J. Gantous, Eccentricity control strategy for geosyncronous communication satellites, in: Space Dynamics for Geostationary Satellites, Cépaduès Éditions, Toulouse, 1986, pp. 693–704.

237

[7] E. Gottzein, W. Fichter, A. Jablonsky, O. Juckenhofel, M. Mittnacht, C. Muller, M. Surauer, Challenges in the control and autonomy of communications satellites, Control Engineering Practice 8 (2000) 409–427. [8] S. Herrick, Astrodynamics, vols. I and II, Van Nostrand Reinhold Co, New York, 1972. [9] A.A. Kamel, C.A. Wagner, On the orbital eccentricity control of syncronous satellites, The Journal of the Astronautical Sciences 30 (1982) 61–73. [10] C.A. Kluever, Geostationary orbit transfers using solar electric propulsion with specific impulses modulation, Journal of Spacecraft and Rockets 41 (2004) 461–466. [11] B.S. Lee, J.S. Lee, K.H. Choi, Analysis of a station-keeping maneuver strategy for collocation of three geoestationary satellites, Control Engineering Practice 7 (9) (1999) 1153–1161. [12] B.S. Lee, K.H. Choi, Collocation of two GEO satellites and one inclined GSO satellite, Aerospace Science and Technology 4 (7) (2000) 507–515. [13] G. Legendre, Le maintien à poste de satellites géostationnaries I. Evolution de L’orbite, in: Cours de Technologie Spatiale: Le mouvement du véhicule spatial en orbite, Centre National D’etudes Spatiales Toulouse, 1980, pp. 583–609. [14] T.S. No, Simple approach to East–West station keeping of geosynchronous spacecraft, Journal of Guidance Control and Dynamics 22 (5) (1999) 734– 736. [15] B.K. Park, M.J. Tahk, H.C. Bang, S.B. Choi, Station collocation design algorithm for multiple geostationary satellites operation, Journal of Spacecraft and Rockets 40 (6) (2003) 889–893. [16] B. Pattan, Delta velocities and propulsion requirements for East–West station keeping for spacecraft in GSO, Space Communitacions 17 (2001) 313–319. [17] J.J. Pocha, An Introduction to Mission Design for Geostationary Satellites, Space Technology Library, Reidel Publishing Company, Dordrecht, 1987. [18] M.J. Sidi, Space Dynamics and Control. A Practical Engineering Approach, Cambridge University Press, Cambridge, USA, 1997. [19] D. Slavinskas, H. Daddaghi, W. Bendent, G. Johnson, Efficient inclination control for geoestationary satellites, Journal of Guidance, Control and Dynamics 11 (6) (1988) 584–589. [20] E.M. Soop, Introduction to geostationary orbits, in: W.R. Burke (Ed.), ESA Scientific and Technical Publication SP-1053, ESTEC, Noordwijk, The Netherlands, 1983. [21] E.M. Soop, Handbook of Geostationary Orbits, Space Technology Library, Kluwer Academic, Dordrecht, 1994. [22] K.N. Srinivasamurthy, N.S. Gopinath, Strategy analysis for collocation of INSAT-2 satellites, Acta Astronautica 50 (6) (2002) 343–349. [23] N.A. Titus, Optimal station-change maneuver for geostationary satellites using constant low thrust, Journal of Guidance, Control and Dynamics 18 (5) (1995) 1151–1155.