# Low-Thrust Out-of-Plane Orbital Station-Keeping Maneuvers for Satellites

Authors: Vivian M. Gomes; Antonio F. B. A. Prado
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

This paper considers the problem of out of plane orbital maneuvers for station keeping of satellites. The main idea is to consider that a satellite is in an orbit around the Earth and that it has its orbit is disturbed by one or more forces. Then, it is necessary to perform a small amplitude orbital correction to return the satellite to its original orbit, to keep it performing its mission. A low thrust propulsion is used to complete this task. It is important to search for solutions that minimize the fuel consumption to increase the lifetime of the satellite. To solve this problem a hybrid optimal control approach is used. The accuracy of the satisfaction of the constraints is considered, in order to try to decrease the fuel expenditure by taking advantage of this freedom. This type of problem presents numerical difficulties and it is necessary to adjust parameters, as well as details of the algorithm, to get convergence. In this versions of the algorithm that works well for planar maneuvers are usually not adequate for the out of plane orbital corrections. In order to illustrate the method, some numerical results are presented.


---

## Body

## 1. Introduction

The present paper studies the orbital maneuvers required by a spacecraft that needs to perform an orbital correction maneuver in order to remain with its orbital elements inside a region where it can achieve the goals of the mission. It is assumed that the orbit of the satellite is given, as well as a nominal orbit for this satellite that allows it to be useful for the planned activities. Then, it is necessary to maneuver this satellite from its current position to the nominal specified orbit. Emphasis will be given in out-of-plane maneuvers, because this situation was not much explored in previous studies, like in [1, 2], and it has a high cost in terms of fuel consumption. Another contribution of the present paper is that it studies the problem of station-keeping maneuvers. In terms of theoretical definition, orbital transfers and station keeping are quite similar. But, when facing practical problems, from the engineering point of view, they are different. This fact can be seen in the literature that presents many methods for station-keeping and orbital transfers, considering  them as two different problems. For station-keeping, it is necessary much lower values for the parameters that control the convergence of the method. It means that convergence is much more difficult in the present case. To solve this problem, it was necessary to divide the convergence of the method in several steps. First a larger number for the tolerance was used and then, using the solution obtained in this step as first guess, we searched for new values of the parameters α1 and β1, in order to make the steps in the direction of the convergence smaller than in the previous case. This technique is repeated until the goal defined for the convergence is reached. So, the present paper shows that this method is also valid for the important application of station keeping maneuvers.The control available to perform this maneuver is the application of a low thrust, and the objective is to perform this maneuver with minimum fuel consumption. An optimal approach will be used. There is no time restriction involved here, and the spacecraft can leave from any point in the initial orbit and arrive at any point in the final orbit. The stochastic version of the projection of the gradient method is used. This version allows us to include the fact that the constraints do not need to be exactly satisfied (see [1, 2]). This is done to realistically treat the numerical inaccuracies and/or flexibilities in terms of tolerance in mission requirements, leading to situations where the final state is constrained to lie inside a given region, instead of having an exact value.


## 2. Review of Orbital Maneuvers

A very important result in this field was obtained by Hohmann [3]. He solved the problem of minimum ΔV transfers between two circular coplanar orbits. The Hohmann transfer would later be generalized to the elliptic case (transfer between two coaxial elliptic orbits) by Marchal [4].After that, can be found in the literature the problem of a two-impulse transfer, where the magnitude of the two impulses are fixed, like in Jin and Melton [5] and Jezewski and Mittleman [6].Then, the three-impulse concept was introduced in the literature by Hoelker and Silber [7] and Shternfeld [8]. They showed that a bielliptical transfer between two circular orbits has a ΔV lower than the one required by the Hohmann transfer, for some combinations of the initial and final orbits. Continuing those studies, Ting [9] showed that the use of more than three impulses does not lower the ΔV, considering only situations where impulsive maneuvers are used.Some other researchers worked on methods where the number of impulses is a parameter to be optimized, and not a value fixed in advance. It is the case of the papers written by Jezewski and Rozendaal [10] and Eckel [11]. Most of the research done in this situation is based on the “Primer-Vector” theory, developed by Lawden  and showed in [12, 13].Later, two more techniques were introduced in the orbital maneuvers field, using the concepts of swing-by and gravitational capture. Those techniques are based on the use of the gravitational forces of a third body to increase or decrease the energy of the spacecraft, so reducing the fuel consumption of the maneuver. References [14–20] describe this problem and show some applications of those techniques in more details.


## 3. Definition of the Problem

The objective of the problem considered here is to change the orbit of a spacecraft. The amplitude of this orbital maneuver has to be compatible with the requirement that the spacecraft needs to return to its nominal orbit that was affected by some perturbations [21, 22], but this is done in a small amplitude framework. Note that the third-body perturbation, when the perturbed body is not in the same plane of the perturbing body, can change the inclination of the perturbed body, as shown in [21]. Then, an initial and a final orbit around the Earth are completely specified. The problem is to find how to transfer the spacecraft between those two orbits and make this transfer in a form that minimizes the fuel consumption. Time restrictions will not be considered in the present investigation, and the spacecraft can leave and arrive at any point in the given initial and final orbits. The maneuver is made by an engine that can deliver a thrust with constant magnitude and variable direction. Although the algorithm is generic regarding the orbital elements changed, some adaption is made in order to search for faster convergence in the case of a maneuver that changes the inclination of the orbit.The spacecraft is assumed to be travelling in a Keplerian orbit without perturbation forces during the maneuver. This orbit is controlled by the thrusts, so the orbit is no longer keplerian when the engines are turned on. In a situation like this, there are two types of motion:(i)
a Keplerian orbit, assuming that the Earth is a point of mass and that it dominates the motion of the spacecraft. This motion occurs when the engine is not being applied to the satellite;(ii)
the motion governed by the two forces: the gravity of the Earth and the force delivered by the thrusts. This motion occurs when the engine is working.The thrusts have the following characteristics:(i)
fixed magnitude, so intermediate values are not allowed and the thrusts can deliver the full power or are not working;(ii)
constant ejection velocity that has to occur to be consistent with the assumption that the velocity of the gases ejected from the thrusts is constant;(iii)
free angular motion: the direction of the thrusts is not constant during the transfer. This direction is specified by the anglesu1 and u2, called pitch (the angle between the direction of the thrust and the perpendicular to the line Earth-spacecraft) and yaw (the angle with respect to the orbital plane);


## 4. Optimal Control Formulation

This approach is based on optimal control theory. First-order necessary conditions for a local minimum are used to obtain the adjoint equations the Pontryagin's maximum principle, which allow us to obtain the control angles at each instant of time (in fact, we search for those control angles as a function of the range angle that is an independent variable that replaces the time and is equivalent to the true anomaly). These assumptions lead us to a Two-Point Boundary Value Problem (TPBVP), where the difficulty is to find the initial values of the Lagrange multipliers. The approach used in the present research is the hybrid approach of guessing a set of values, numerically integrating all the differential equations, and then searching for a new set of values, based on a nonlinear programming algorithm. With this approach, the problem is reduced to a parametric optimization, where the Lagrange multipliers are variables to be optimized. This idea follows the same principles used by Biggs in [23, 24].The method also showed by Biggs [24], where the “adjoint-control” transformation is performed and, instead of the initial values of the Lagrange multipliers, one guesses the control angles and their rates at the beginning of the “burning arc,” is used here. Then, it is easier to find a good initial guess, and the convergence is faster. This hybrid approach has the advantage that, since the Lagrange multipliers remain constant during the “ballistic arcs,” it is necessary to guess values of the control angles and their rates only for the first “burning arc.” This transformation reduces very much the number of variables to be optimized and, as a consequence, the time of convergence.The next paragraphs show how to write this problem using an optimal control approach.Objective Function
LetMf be the final mass of the vehicle. It has to be maximized with respect to the control u̲(·);
Subject to(4.1)ẋ̲̲=f̲(x̲,u̲,s),Ce̲(x̲,u̲,s)=Ee̲,Cd̲(x̲,u̲,s)≤Ed̲,h̲(x̲(tf),(tf))=Eh̲,t0,x̲(t0)given,
where x̲ is the state vector, f̲(·) is the right-hand side of the equations of motion, in the same way that was used by Biggs [23, 24] and by Prado and Rios-Neto [25] and Bryson and Ho [26] s is the independent variable that replaces the time (s0≤s≤sf),Ce̲(·) and Cd̲(·) are the algebraic dynamic constraints on state and control that have dimensions me and md,h̲(·) are the boundary constraints of dimension mh, and Ee̲,Ed̲,Eh̲ are error vectors satisfying
(4.2)|Eei|≤EeiT,i=1,2,3,…,me,|Edi|≤EdiT,i=1,2,3,…,md,|Ehi|≤EhiT,i=1,2,3,…,mh,
where the fixed given values EeiT,EdiT, and EhiT characterize the region around zero within which errors are considered tolerable.
To avoid singularities or close approaches to them, we use the following variables:(4.3)X1=[a(1-e2)μ]1/2,X2=ecos⁡(ω-ϕ),X3=esin⁡(ω-ϕ),X4=FuelConsumedm0,X5=t,X6=cos⁡(i2)cos⁡(Ω+ϕ2),X7=sin⁡(i2)cos⁡(Ω-ϕ2),X8=sin⁡(i2)sin⁡(Ω-ϕ2),X9=cos⁡(i2)sin⁡(Ω+ϕ2),
where: a = semi-major axis, e = eeccentricity, i = inclination, Ω = argument of the ascending node, ω = argument of perigee, f̲ = true anomaly, s = range angle, ϕ=f̲+ω-s,μ = gravitational constantl, m0 = initial mass of the spacecraft.
In those variables, the equations of motion are the following  ones:(4.4)dX1ds=f1=SiX1F1,dX2ds=f2=Si{[(Ga+1)cos⁡(s)+X2]F1+GaF2sin⁡(s)},dX3ds=f3=Si{[(Ga+1)sin⁡(s)+X3]F1-GaF2cos⁡(s)},dX4ds=f4=SiGa̲F(1-X4)X1W,dX5ds=f5=SiGa(1-X4)m0X1,dX6ds=f6=-SiF3[X7cos⁡⁡(s)+X8sin⁡(s)]2,dX7ds=f7=SiF3[X6cos⁡(s)-X9sin⁡(s)]2,dX8ds=f8=SiF3[X9cos⁡(s)+X6sin⁡(s)]2,dX9ds=f9=SiF3[X7sin⁡(s)-X8cos⁡(s)]2,
where
(4.5)Ga=1+X2cos⁡(s)+X3sen(s),Si=μX14[Ga3m0(1-X4)].F,F1,F2,F3  are the forces generated by the thrust and they are given by
(4.6)F⃗=F⃗1+F⃗2+F⃗3,|F⃗|=F,F1=Fcos⁡(α)cos⁡(β),F2=Fsen(α)cos⁡(β),F3=Fsen(β),
where α is the angle between the direction of the thrust and the line that is perpendicular to the radius vector and β is the angle between the direction of the thrust and the orbital plane.The equations for the Lagrange multipliers are given as follows:(4.7)dp1ds=-4∑j=19pjfj+p1f1-p4f4-p5f5X1,dp2ds=cos⁡(s)Ga[3∑j=19pjfj-p4f4-p5f5]-Sip2F1-Sicos⁡2(s)(p2F1-p3F2)-Sicos⁡(s)sen(s)(p2F2+p3F1),dp3ds=sen⁡(s)Ga[3∑j=19pjfj-p4f4-p5f5]-Sip3F1-Sicos⁡(s)sen(s)(p2F1-p3F2)-Sisen2(s)(p2F2+p3F3),dp4ds=-[∑j=19pjfj-p4f4-p5f5m0(1-X4)],dp5ds=0,dp6ds=-SiF3[p7cos⁡(s)+p8sen(s)]2,dp7ds=SiF3[p6cos⁡(s)-p9sen(s)]2,dp8ds=SiF3[p6sen(s)+p9cos⁡(s)]2,dp9ds=-SiF3[p8cos⁡(s)-p7sen(s)]2.
The control that is applied to the spacecraft also needs a substitution of variables, with the objective of solving numerical problems. In this situation we use the following set of variables:(4.8)u1=s0,u2=(sf-s0)cos⁡(β0)cos⁡(α0),u3=(sf-s0)cos⁡(β0)sin⁡(α0),u4=(sf-s0)sin⁡(β0),u5=α′,u6=β′.
The first-order necessary conditions of the optimal problem, which are the conditions that allow us to obtain the optimal control, can be obtained, at every instant of time, by extremizing the Hamiltonian of the system. The equations are given as follows: (4.9)sin⁡(α)=q2S′,sin⁡(B)=q3S′′,cos⁡(α)=q1S′,cos⁡(B)=S′S′′,
where
(4.10)S’=±[q12+q22]1/2,S′′=±[q12+q22+q32]1/2,q1=p1X1+p2[X2+(Ga+1)cos⁡(s)]+p3[X3+(Ga+1)sin⁡(s)],q2=p2Gasin⁡(s)-p3Gacos⁡(s),q3=-p6[X7cos⁡(s)+X8sin⁡(s)]2+p7[X6cos⁡(s)-X9sin⁡(s)]+p8[X6sin⁡(s)+X9cos⁡(s)]+p9[X7sin⁡(s)-X8cos⁡(s)].
It is possible to include constraints in this type of problem. They can be represented by the following equations:
(4.11)S(⋅)≥0,(a-a*)|a0-a*|=0,[a(1+e)-a*(1+e*)]|a0(1+e0)-a*(1+e*)|=0,(i-i*)|i0-i*|=0,(Ω-Ω*)|Ω0-Ω*|=0,(ω-ω*)|ω0-ω*|=0.
In those constraints, the first equation represents a generic inequality constraint and the other five represent the important constraint of having an orbit specified.


## 5. Numerical Method

To solve this problem, the stochastic version of the projection of the gradient method (Rios-Neto and Pinto [1]) can be used. It is showed in the next paragraphs.We start by assuming an initial value forp̲̅, the vector of parameters p̲ that represents the variables that we are searching. It may come from an initial guess or from a value that belongs to an immediately previous iteration. Then, a first-order direct search approach is adopted in a typical iteration to obtain an approximate solution for the increment Δp̲:(5.1)Minimize:J(p̲̅+Δp̲)(5.2)Subjectto:Ce̲(p̲̅+Δp̲)=α1Ce̲(p̲̅)+Ee̲,(5.3)Cd̲(p̲̅+Δp̲)=β1Cd̲(p̲̅)+Ed̲,where J(p̲) is the objective function, Ce̲(p̲) is the equality constraints, Cd̲(p̲) is the active inequality constraints at p̲, and 0≤α1<1,0≤β1<1 are parameters that have to be chosen close enough to one to make the increments Δp̲ be of the first order of magnitude.Linearizing the left-hand sides of (5.2) and (5.3) and using a stochastic interpretation for the errors Ee̲ and Ed̲, we have(5.4)(α1-1)Ce̲(p̲̅)=(d[Ce(p̲̅)]dp̲)Δp̲+Ee̲,(5.5)(β1-1)Cd̲(p̲̅)=(d[Cd(p̲̅)]dp̲)Δp̲+Ed,
where Ed̲ and Ee̲ are zero mean uniformly distributed errors, given by(5.6)E[EeEe̲T]=diag⁡[ei,i=1,2,…,me],E[EdEd̲T]=diag⁡[di,i=1,2,…,md],
where E[·] indicate the expected value of its argument.The condition shown in (5.1) is approximated by(5.7)-g⋅∇̲JT(p̲̅)=Δp̲+n̲,
where g≥0 needs to be adjusted to guarantee that Δp̲ is small enough to permit the use of the linearized representation of J(p̲̅+Δp̲), and n̲ is taken as a zero mean uniformly distributed random vector, modeling the a priori searching error in the direction of the gradient ∇̲J(p¯̲), with(5.8)E[nn̲T]=P̲̅
as its diagonal covariance matrix. The values of the variances in P̲̅ are chosen such as to characterize an “adequate order of magnitude” for the dispersion of n̲. The diagonal form adopted is to model the assumption that it does not impose any a priori correlation between the errors in the gradient components.The simultaneous consideration of conditions of (5.4), (5.5), and (5.7) characterizes a problem of parameter estimation, which can be written as(5.9)U̲̅=U̲+n̲,(5.10)Y̲¯=MU̲+V̲,
where U̲̅Δ̲-g·∇JT(p̲̅) is the “a priori information”; U̲Δ̲Δp̲;Y̲Δ̲[(α1-1)CeT(p̲̅):(β1-1)CdT(p̲̅)] is the observation vector; M̲TΔ̲[(d(Ce̲(p̲̅))/dp̲)T:(d(Cd̲(p̲̅))/dp̲)T];V̲T=[Ee̲T:Ed̲T].Adopting a criterion of linear, minimum variance estimation, the optimal search increment can be determined using the classical Gauss-Markov estimator, which in Kalman form gives(5.11)Û̲=U̲̅+K̲(Y̲-M̲U̲̅),(5.12)P̲=P̲̅-KM̲P̲̅,(5.13)K̲=P̲̅M̲T(M̲P̲̅M̲T+R̲)-1,
where P̲̅ is defined as before, R̲Δ̲E[VV̲T] = diag [Rk,k=1,2,…,me+md], and P̲ has the meaning of being the covariance matrix of the errors in the components estimates of U̲, that is,(5.14)P̲=E[(U̲-Û̲̲)(U̲-Û̲̲)T].To build a numerical algorithm using the proposed procedure, the following types of iterations are considered:(i)
initial phase of acquisition of constraints: when starting from a feasible point that satisfies the inequality constraints, the search is done to capture the equality constraints, including those inequality constraints that eventually became active along this phase;(ii)
search of the minimum: when from a point that satisfies the constraints in the limits of the tolerable errorsV̲ in (5.10), the search is done to take the objective function (5.1) to get closer to the minimum; this search is conducted by relaxing the order of magnitude of the error bounds around the constraints;(iii)
restoration of the constraints: when from a point that resulted from a type (ii) iteration, the search is done to restore constraints satisfaction, within the limits imposed by the errorV̲ in (5.10).Rios-Neto and Pinto [1] suggest how to choose good values for the numerical parameters that must be different for each type of iteration.


## 6. Simulations and Numerical Tests

To verify if the algorithm proposed is useful for this particular type of problem, three maneuvers of station keeping were simulated, all of them having modifications in the inclination. These results were compared with the ones obtained by the deterministic version, without flexibility in the constraint’s satisfaction. Similar problems can be found in [27–29]. The first maneuver will occur with the data given in Table 1. The thrust level is 2.0 N. Table 2 shows the errors allowed in the final Keplerian elements of the orbit. The goal is to combine a modification in semimajor axis, eccentricity, and inclination in a single maneuver.Table 1
Data for Maneuver 1.
OrbitsInitialFinalSemimajor axis7060.007090.00Eccentricity0.030.01Inclination (degrees)1.000.50Ascending node (degrees)0.00FreeArgument of perigee (degrees)0.00FreeMean anomaly (degrees)0.00FreeTable 2
Errors allowed for final keplerian elements for Maneuver 1.
Semimajor axis1.0 KmEccentricity0.005Inclination0.001 degSeveral numbers of “burning arcs” were used for the same maneuver, in order to obtain some information about the importance of this parameter.The consumptions found are showed in Table3, as well as comparisons with deterministic methods that we also simulated for comparison.Table 3
Fuel expenditure comparisons (kg) for Maneuver 1.
ApproachStochasticDeterministic2 Arcs0.3010.3033 Arcs0.3000.3014 Arcs0.2980.2995 Arcs0.2970.2986 Arcs0.2960.2977 Arcs0.2950.2968 Arcs0.2940.296The second maneuver will occur with the data given in Table4. Now, only the inclination is changed. The thrust level is again 2.0 N. Table 5 shows the errors allowed in the final Keplerian elements of the orbit.Table 4
Data for transfer Maneuver 2.
OrbitsInitialFinalSemimajor axis8000.008000.00Eccentricity0.020.02Inclination (degrees)3.002.50Ascending node (degrees)0.00FreeArgument of perigee (degrees)0.00FreeMean anomaly (degrees)0.00FreeTable 5
Errors allowed for final keplerian elements for Maneuver 2.
Semimajor axis0.50 KmEccentricity0.005Inclination0.001 degThe choice of the number of “burning arcs” was done for several different values, in the same way that was made in the first maneuver.The consumptions found are showed in Table6, as well as comparisons with deterministic methods.Table 6
Fuel expenditure comparisons (kg) for Maneuver 2.
ApproachStochasticDeterministic2 Arcs0.1300.1323 Arcs0.1280.1304 Arcs0.1270.1295 Arcs0.1270.1286 Arcs0.1270.1287 Arcs0.1270.1288 Arcs0.1270.128The third maneuver occurs with the data given in Table7. Now, inclination and eccentricity are changed. The thrust level is again 2.0 N. Table 8 shows the errors allowed in the final Keplerian elements of the orbit.Table 7
Data for transfer Maneuver 3.
OrbitsInitialFinalSemimajor axis8000.008000.00Eccentricity0.050.02Inclination (degrees)3.002.50Ascending node (degrees)0.00FreeArgument of perigee (degrees)0.00FreeMean anomaly (degrees)0.00FreeTable 8
Errors allowed for final keplerian elements for Maneuver 3.
Semimajor axis0.50 KmEccentricity0.005Inclination0.001 degThe number of “burning arcs” was again varied from two to eight, in the same way that was made in the other maneuvers. The consumptions found are showed in Table9, as well as comparisons with deterministic methods.Table 9
Fuel expenditure comparisons (kg) for Maneuver 3.
ApproachStochasticDeterministic2 Arcs0.1700.1723 Arcs0.1670.1694 Arcs0.1650.1675 Arcs0.1630.1666 Arcs0.1620.1607 Arcs0.1620.1608 Arcs0.1620.160Next, a new maneuver is considered to emphasize the effects of the errors allowed for the final Keplerian elements. This maneuver considers the data shown in Table1 for the initial and final orbits, but increases the tolerable errors to the values shown in Table 10.Table 10
Errors allowed for final keplerian elements for Maneuver 4.
Semimajor axis2.0 KmEccentricity0.05Inclination0.01 degSeveral numbers of “burning arcs” were also  used for this maneuver, to be compatible with the previous studies. The consumptions found are showed in Table11, as well as the usual comparisons with the deterministic methods. Note that the difference in savings increased, when larger values for the errors are allowed.Table 11
Fuel expenditure comparisons (kg) for Maneuver 4.
ApproachStochasticDeterministic2 Arcs0.2950.3033 Arcs0.2940.3014 Arcs0.2920.2995 Arcs0.2910.2986 Arcs0.2900.2977 Arcs0.2900.2968 Arcs0.2890.296


## 7. Conclusions

Optimal control was used to build an algorithm to search for solutions for the problem of minimum fuel consumption to make orbital maneuvers for a satellite that needs to return to its nominal orbit after deviations caused by perturbations forces.In this research, emphasis was given in the problem considering station keeping with out-of-plane maneuvers. The adjustments made in the algorithm, as well as in the parameters used, allow us to get convergence in most of the cases for this version of the problem.This problem took into account the accuracy in the constraint's satisfaction by using the nonlinear programming algorithm proposed by Rios Neto and Pinto [1].The results showed that it is possible to reduce the costs by exploring tolerable errors shown in Tables2, 5, 8 and 10 in the constraint's satisfaction. The amount saved can be important in many cases.It is also clear that increasing of the number of propulsion arcs decreases the fuel costs. It can be seen from Tables3, 6, 9 and 11, because the fuel consumption is smaller for larger numbers of arcs. The reason is that increasing this number causes an increase in the degrees of freedom available for the optimization technique. Since the maneuvers have small amplitudes, this increase in the burning arcs has a limit in the savings, so the fuel consumption reaches a constant value after a certain value.



---

*Source: 532708-2012-03-22.xml*