{'fraction_non_alphanumeric': 0.029315162507966078, 'fraction_numerical': 0.010980930437766557, 'mean_word_length': 4.7174887892376685, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimal Control for Minimum-Fuel Geostationary Station Keeping of Satellites Equipped with Electric Propulsion

: In this paper, a control scheme is elaborated to perform the station keeping of a geostationary satellite equipped with electric propulsion. The use of electric thrusters imposes to take into account some additional non linear constraints that make the overall station keeping optimal control problem diﬃcult to solve directly. That is why we choose here to decompose the station keeping problem in two control problems. The ﬁrst one consists in solving a classical optimal control problem with an indirect method initialized by a direct method without using the thrusters operational constraints. The second problem deals with the thrusters operating constraints, that are taken into account by two diﬀerent ways. Simulation results validate the eﬀect of the optimal control thrusts obtained with these methods.


INTRODUCTION
Due to orbital disturbing forces, any satellite in Geostationary Earth Orbit (GEO) drifts outside its station keeping (SK) window (a rectangular box of a given geographical longitude and latitude range). Performing an accurate SK strategy is therefore necessary to compensate for the induced environmental secular and periodic disturbances and GEO satellites are equipped with electric and/or chemical thrusters.
Chemical propulsion systems have been and are widely used. For these propulsion systems with high thrust capabilities, SK control laws are usually designed assuming an impulsive idealisation of the thrust, as described for example in Soop (1994). The idea of using electric propulsion for station keeping dates back to the sixties (see, e.g. Barrett (1967) and Hunziker (1970) ) and some theoretical developments have been presented in the eighties by Anzel (1988) and Eckstein (1980). Nowadays the electric propulsion is a viable alternative to the chemical one, in particular in the case of SK of GEO satellite, to take thrust operations constraints into account: large onboard power needs, mission requirements restricting the duration of use of the electric power system, impossibility to perform SK maneuvers at eclipse epochs. Moreover, its bigger specific impulsion and the consequent savings in fuel consumption, leads to a reduction of the satellite mass and enables increased payload capacity and/or improved satellite longevity.
Considering these technological and operational features, optimal control strategies for electric SK, taking various constraints into account (minimum elapsed time between two consecutive firings, on-off profile of the thrusters, thrust allocation) have to be carefully designed. The problem of station keeping is in general expressed as an optimal control problem even if the above mentioned constraints, inherent to the use of electric thrusters, prevent us to solve it with classical methods.
Several types of methods may be used to solve optimal SK control problems. When simple models are used to describe the disturbing forces, analytical models provide sufficient control laws, as in Sukhanov and Prado (2012). Otherwise, it is necessary to resort to numerical methods, such as direct collocation based methods as described in Hull (1997) and Betts (1998). For this family of approaches, the state and the control variables are discretized to produce a non linear programming problem and get an optimal open loop control. To deal with on-off models of thrusts, the Pulse Width Modulation technique may be used to generate rectangular profiles from a continuous one (see Vazquez et al. (2015) and the references therein). Losa et al. (2005) has formulated a method based on differential inclusion, and a first avenue for the use of decomposition methods to solve the problem is given in Losa et al. (2006).
In this paper, the idea is to decompose the overall SK optimal control problem into two sub-problems for which appropriate solution methods are designed. In a first step, an indirect method based on the application of Prontryagin Maximum Principle (PMP) with mixed controlstate constraints is applied to solve a simplified optimal SK control problem, i.e. without considering some hard constraints on the control law (thrust constraints such as latency between two bursts of the same thruster and no simultaneous thrusting for instance). The solution of the Two Point Boundary Value Problem (TPBVP) derived from the optimal necessary conditions is initialized by an approximate solution given by a collocation based direct method. In a second step, a numerical approach is used to enforce all the thrust constraints left apart at the first step. A realistic numerical example illustrates the efficiency of the proposed approach.

Minimum-Fuel Station Keeping Problem
Let consider a satellite equipped with 4 electric thrusters mounted on the anti-nadir face. The position of the satellite on its orbit is described with the equinoctial orbital elements as defined in Battin (1999): where a is the semi-major axis, (e x , e y ) the eccentricity vector components, (i x , i y ) the inclination vector components, M Θ = ω + Ω + M − Θ is the mean longitude where Ω is the right ascension of the ascending node, ω is the perigee's argument, M is the mean anomaly and Θ(t) is the right ascension of the Greenwich meridian. The dynamics of the satellite may be represented by the following non linear state-space model: where f L ∈ R 6 is the Lagrange contribution part of the external force model described by the CNES ORANGE model (cf. Campan and Brousse (1994)) and f G ∈ R 6×3 is the Gauss contribution part.
is the control vector expressed in the local orbital RT N frame (also written RSW ) defined in Soop (1994) by: • N is the unit vector along the kinetic momentum; • R is the unit vector along the direction Earth's center -satellite; • T makes orthogonal direct basis.
In order to deal with the station keeping problem, the relative state of the satellite with respect to the station keeping state is rather used. The station keeping state is defined as follows: where a sk is the synchronous semi-major axis and M Θ sk is the station mean longitude.
The relative dynamics equations are developed by linearization of equation (2) about the station keeping point (3). By denoting x = x eoe − x sk the relative state model for the SK problem reads: where matrices A ∈ R 6×6 , B ∈ R 6×3 , C ∈ R 3×6 and D ∈ R 6 are obtained from the linearization of functions f L and u → f G u.
Remembering that 4 thrusters are available to realize the control, the transcription of the station keeping problem can be expressed considering the 4 thrusts provided by the 4 engines in the satellite dynamics. The control u(t) is a linear combination of the 4 thrusts such that u = The thrust direction matrices Γ j ∈ R 3 are defined such that: where the angles θ j and α j are defined exactly as in Anzel (1988).
Two other transformations are proposed. Firstly, for the sake of simplicity, we will normalize the thrust vector: Secondly, as the thrusts are on-off, the thrust profile is modeled as a rectangular signal that is parametrized by the date t i,j corresponding to the middle instant of the thrust and by its half width duration denoted ∆t i,j as shown on Figure 1: where P i is the number of thrusts of thruster i.
The main goal of the station keeping problem is to maintain the longitude and the latitude of the satellite in a box defined by its size δ by acting on the orbital parameter via the 4 thrusters. The associated optimal control problem is in general defined over a fixed horizon for the computation of optimal open loop control laws. In this context, optimality means that a minimum fuel-solution is looked for to extend the operational life time of the satellite. Therefore, performing minimum-fuel station keeping is minimizing the performance index: Constraints have to be imposed so that the satellite stay in the station keeping geographical box. This box is defined in the plane (latitude,longitude) of width 2δ × 2δ centered on the station keeping geographical position y sk . These constraints on the output variables are transformed into constraints on the state variables: The initial position is chosen to be at the center of the station keeping box, i.e. x(0) = 0. The station keeping has to be performed on the time interval [0, T ]. In order to use the same control law on the intervals [kT, (k + 1)T ], k ∈ N * , it is convenient to add the terminal condition : x(T ) = 0.

Operational Constraints on Actuation
Beside the station keeping geographical constraints and the usual bounds on the maximum thrust, some additional technological operational constraints on the actuation have to be taken into account: (i) thrusters cannot have simultaneous thrusts; (ii) a thrust must last at least T l : 2∆t i,j T l ; (iii) two successive thrusts of a given thruster must be separated of an interval of latency equal to T s ; (iv) two thrusts of two different thrusters must be separated by an interval of latency equal to T d .
In order to give a tractable mathematical expression for constraints (iii) and (iv), let us define for the thruster i the ordered sequence of firing times (t i,k ) k=1...Pi in increasing order. The constraint for the time latency between the thrust k of thruster i and the thrust l of thruster j is thus mathematically expressed as: The operational constraint (10) raises some difficult mathematical issues. Firstly, due to the parametrization of the rectangular functions by t i,j and ∆t i,j , the firing times and durations are both optimization variables and intrinsic variables of the problem to solve. Secondly, this constraint is non convex and logical, and is therefore difficult to be tackled within an optimal control problem.
In addition, some other convenient constraints are enforced in order to make sure that the thrusts will not begin before 0 and end after T :

Decomposition of the Overall Problem
Considering all the operational constraints described above, the minimum-fuel SK problem to solve may be summarized as the following optimal control problem: It is possible to solve the whole problem with direct methods such as collocation methods as described in Hull (1997) and Betts (1998). However, to obtain a solution with this method, the number of bursts per thruster has to be known beforehand. To solve this problem in its whole generality, it is necessary to solve a non linear mixed integer optimization problem with respect to the thrusts number for each thruster. Moreover, the solution is very sensitive to the number of collocation points and the initial guess of the optimization variables.
An alternative strategy for solving the Station Keeping Problem 1 consists in splitting it into two different steps: (A) The thrusters constraints (10) and (11) are removed so that only an optimal control problem (OCP) with state and control constraints remains. This particular OCP is tackled via a hybrid approach relying on an indirect method initialised by a direct method solution dedicated to the search of adjoint variables as described in Grimm and Markl (1997) and von Stryk and Bulirsch (1992). The Pontryagin Maximum Principle (PMP) is applied to derive necessary conditions of optimality and the associated Two Points Boundary Value Problem (TPBVP). The TPBVP is then solved via a collocation method for which the results of the direct collocation step serves as first guess. (B) As the result of the first step produces a control law that does not necessarily respect the thruster operational constraints, a second part is needed in order to obtain modified results of the first part fitting the thrusters constraints.
To sum up, instead of solving the overall problem at once, this problem is split into two smaller and simpler problems that are more easily solved in general as in Losa et al. (2006).

SOLUTION OF THE SIMPLIFIED OCP (STEP (A))
Removing the operational actuation constraints (10) and (11) allows to redefine the optimal control problem. In particular, the thrust functionsF i are not a priori modeled as rectangular functions parametrized by t i,j and ∆t i,j and the simplified OCP to be solved reads as : Problem 2. The simplified OCP to be solved is finding: Problem 2 is a minimum-fuel linear OCP with constraints both on the state and the control vectors and defined on a fixed horizon. To solve this problem, a direct collocation method finds a first approximation of the optimal solution. This approximated solution is then used to initialize the solution of the TPBV problem obtained by the application of Pontryagin Maximum Principle (PMP) and the derivation of the necessary optimality conditions. Using an indirect method to solve an OCP leads to a very precise solution which is very sensitive to the initial conditions. A good guess for this initial condition is therefore needed. It is provided by the direct collocation method.

Finding an Initial Guess via a Direct Method
Suboptimal solution of the Problem 2 is found by applying a classical collocation method. The time interval [0, T ] is divided into N equidistant points: 0 = τ 1 < τ 2 < · · · < τ k < · · · < τ N −1 < τ N = T, (14) with τ k+1 − τ k = ∆τ so that a grid of N points defines the discretization of Problem 2, with x k = x(τ k ) and F k =F (τ k ). The collocation method aims at finding the value of the control at the points τ k of the grid. The continuous control is then obtained by interpolation.
The collocation problem is mathematically formulated as follows: Problem 3. The collocation problem to be solved to initialize the optimal control Problem 2 is finding: where the defects R k are obtained by discretizing the differential constraint of Problem 2 using the 4th order Simpson algorithm described in Hull (1997). The reference Grimm and Markl (1997) shows that these Lagrange parameters are an approximation of the adjoint vector for the OC Problem 2. Therefore, (x 1 , λ 1 ) can be used as an initial condition for the TPBV problem obtained from the PMP and from the derivation of the necessary optimality conditions.

Solving the Simplified OCP via an Indirect Method
The idea is now to tackle the usual issues of the direct methods (sensitivity of the solution to the initial guess) by combining it with an indirect approach based on the PMP. To handle the state geographical constraints, the technique is based on a penalty function as described in Naidu (2002). The following penalty function is added to the cost function of Problem 2 (see Equation (18)) and is chosen such that: The Hamiltonian of the system becomes:

t) t A(t) x(t) + D(t) +B(t)F (t) + µ g(x(t)), (17)
where λ(t) ∈ R 6 is the adjoint vector, and µ is a constant parameter that can be chosen as large as needed.
Applying the PMP described in Naidu (2002), the minimization condition:

H x(t),F (t), λ(t) = min u∈[0,1] 4 H (x(t), u(t), λ(t))
(18) leads to the following switching conditions: It is assumed that singular arcs are not possible optimal solutions for this problem. Another way to take the state constraint into account in the derivation of the optimality conditions with the PMP is presented in Hartl et al. (1995).
By applying the first order necessary conditions: and the transversality conditions as boundary conditions, a two points boundary value problem is obtained and formulated as Problem 4.
• To solve this problem, the initial guess used is the one obtained as a solution of Problem 3.
The initial conditions for the state and adjoint vectors being well chosen thanks to the solution of Problem 3, Problem 4 is easily solvable. However, this solution does not satisfy in general the operational actuation constraints described in Section 2.2. The next section proposes an additional complementary step for which a direct method is applied on an auxiliary problem enforcing the actuation constraints on a new and equivalent control law, while preserving the structure and overall effect of the thrust.

ENFORCING THE OPERATIONAL THRUSTERS CONSTRAINTS (STEP (B))
The solution of the TPBV Problem 2 consists in a series of rectangle signals naturally exhibited as a result from application of the PMP. However, this control law does not respect the actuation constraints left apart in the first part. For instance, simultaneous activation of two different thrusters may occur. The aim of this second part is to find an equivalent control law satisfying the required non linear operational constraints. To do so, it is first necessary to clearly define the meaning of equivalent control laws in the sequel. Two different notions of equivalent control laws will be used hereafter.
The first notion of equivalence between two control laws relies on the fuel consumption argument: the goal is to compute a raw control profile that has the same fuel consumption as the profile obtained by solving the TPBV problem. LetF BV P be the control obtained by solving Problem 4. Finding a consumption based equivalent (CBE) control for the satellite is then equivalent to solve Problem 5 defined as follows: Problem 5. Find: The second way to obtain an equivalent control respecting the actuation constraints is to define an "effect-based" equivalent (EBE) control. To get a control profile that has the same effect at time T as the solution of Problems 4 and respecting the actuation constraints, Problem 6 has to be solved and is defined as follows: Problem 6.
such that the constraints (23) Antsaklis (2003). This problem is a non linear optimization problem. Both Problem 5 and 6 can be solved by classical non linear optimization solvers.

NUMERICAL RESULTS
In this section, simulation results obtained with the proposed methodology are presented. Let consider a satellite of mass 4850 kg equipped with 4 electric thrusters oriented in the directions North-East, North-West, South-East and South-West. This satellite has to be controlled in order to remain close to its geostationary position at a fixed longitudeλ. The decomposition presented in Sections 3 and 4 is now illustrated numerically. For the second part, parameter µ has to be chosen. µ = 1.10 4 for the CBE and µ = 1.10 12 for the EBE.
On Figure 2, the geographical parameters obtained by solving Problem 3, and the geographical parameters obtained by solving the TPBV Problem 4 are drawn. Furthermore, Figure 2 shows the latitude and longitude obtained often the CBE and the EBE optimization. It is recalled that the solution of Problem 3 is used as an initial guess for solving Problem 4. Both equivalence schemes of step (B) do not lead to identical results. The effect-based one leads to a trajectory that remains close to the center of the station keeping window whereas the consumption based one leads to a trajectory that goes close to the border of the window. Figure 3 shows a detailed view of the the control profile obtained after solving Problem 4 and the firing sequence for both equivalence scheme. As explained above, the thrusters constraints are not satisfied after step (A), as thrusters NE and SE are simultaneously active, as well as thrusters NW and SW. This issue requires to perform the second step to make the control profile feasible with respect to the thrusters operational constraints. The consumption is very different for both equivalence schemes: the control profile computed with the CBE requires more fuel than the control profile computed with the EBE. The relative difference is of 81, 6%. Despite these differences, the two control profiles solve the station keeping problem, as can be seen in Figure 4.

CONCLUSION
In this paper, a decomposition of the overall station keeping optimal control problem under many operational constraints is used to take into account some difficult constraints inherent to the use of electric propulsion. Firstly, a classical optimal control problem is solved with state constraints using a precise indirect method initialized by a collocation based direct method. Secondly, two ways of dealing with the thrusters operational constraints are proposed resulting in two different fuel consumption results. Despite the positive results presented in this paper, different issues remain open: the optimality of possible singular arcs has to be studied and alternative formulation of the PMP for state constraints optimal control problems (as in Hartl et al. (1995)) should be considered.