Clément Gazzino

Denis Arzelier

Luca Cerri

Damiana Losa

Christophe Louembet

Christelle Pittet

A Three-step Decomposition Method for Solving the Minimum-Fuel Geostationary Station Keeping of Satellites Equipped with Electric Propulsion

Keywords: GEO Station Keeping, electric propulsion, fuel optimal control problem, Pontryagin Maximum Principle, Switched Systems Theory

In this paper, a control scheme is elaborated in order to perform the station keeping of a geostationary satellite equipped with electric propulsion while minimizing the fuel consumption. The use of electric thrusters imposes to take into account some additional non linear and operational constraints that make the overall station keeping optimal control problem difficult to solve directly. That is why the station keeping problem is decomposed in three successive control problems. The first one consists in solving a classical optimal control problem with an indirect method initialized by a direct method without enforcing the thrusters operational constraints. Starting from this non feasible solution for the genuine problem, the thrusters operating constraints are incorporated in the second problem, whose solution produces a feasible but non optimal control profile via two different ways. Finally, the third optimizes the commutation times thanks to a method borrowed to the switched systems theory. Simulation results on a realistic example validate the benefit of this particular control scheme in the reduction of the fuel consumption for the geostationary station keeping problem.

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 to this end GEO spacecraft are equipped with electric and/or chemical thrusters.

Chemical propulsion systems have been and are still widely used. For these propulsion systems with high thrust capabilities, SK control laws are usually designed assuming an impulsive idealization of the thrust, as described for example in [START_REF] Campan | Les techniques de maintien à poste de satellites géostationaires[END_REF][START_REF] Soop | Handbook of Geostationary Orbits[END_REF][START_REF] Sidi | Spacecraft Dynamics and Control[END_REF]. This type of propulsion is still used for the derivation of autonomous SK laws (see for instance the references [START_REF] Chao | Onboard stationkeeping of geosynchronous satellite using a global positioning system receiver[END_REF][START_REF] No | Simple Approach to East-West Station Keeping of Geosynchronos Spacecraft[END_REF][START_REF] Emma | Algorithm for Autonomous Longitude and Eccentricity Control for Geostationary Spacecraft[END_REF][START_REF] Park | A new approach to on-board stationkeeping of GEOsatellites[END_REF][START_REF] Gomes | Low-thrust out-of-plane orbital station-keeping maneuvers for satellites[END_REF]). The idea of using electric propulsion for station keeping dates back to the sixties. The references [START_REF] Boucher | Electric Propulsion for Control of Stationary Satellites[END_REF][START_REF] Molitor | Ion Propulsion System for Stationary-Satellite Control[END_REF][START_REF] Barrett | On the Application of electric Propulsion to Satellite Orbit Adjustement and Station Keeping[END_REF] describe generalities about the SK of a GEO satellite, and the reference [START_REF] Hunziker | Low-Thrust Station Keeping Guidance for a 24-Hour Satellite[END_REF] designs a control strategy in a simplified case. Some theoretical developments have then been conducted in the eighties. The reference [START_REF] Eckstein | Optimal Station Keeping by Electric Propulsion With Thrust Operation Constraints[END_REF] analyses the effect of electric North/South correction with constraints on the firing times while [START_REF] Eckstein | Station acquisition and Station-Keeping with Low Thrust Systems[END_REF] model the secular evolution of the satellite and perform station keeping with a fixed number of thrusts. [START_REF] Eckstein | Optimal autonomous station keeping of geostationary satellites[END_REF] derive a control strategy distinguishing short and long horizons. The patent [START_REF] Anzel | Method and apparatus for a Satellite Station Keeping[END_REF] proposes a new configuration with four thrusters mounted on the anti-nadir face of the satellite and derives a correction strategy while analyzing the effect of the orbital disturbances. This reference also describes the advantages (complete control of the three orbit vectors, remaining diagonal pair of thruster for a complete control in case of failure of one thruster) of this specific propulsion system based on HS 702 satellite with a xenon ion propulsion system. It has been widely used as a typical configuration for electric propulsion in the literature [START_REF] Gopinath | Optimal low thrust orbit transfer from gto to geosynchronous orbit and stationkeeping using electric propulsion system[END_REF], [START_REF] Losa | Electric Station Keeping of Geostationary Satellites: a Differential Inclusion Approach[END_REF], [START_REF] Losa | Station Keeping of Geostationary Satellites with On-Off Electric Thrusters[END_REF] and more recently [START_REF] Walsh | MPC for Coupled Station Keeping, Attitude Control, and Momentum Management of Low-Thrust Geostationary Satellites[END_REF]. Nowadays the electric propulsion is a viable alternative to the chemical one, in particular in the case of SK of GEO satellite ( [START_REF] Peukert | OHB-System's View on Electric Propulsion Needs[END_REF]), despite limiting thrust operations constraints (large on board power needs, mission requirements restricting the duration of use of the electric power system, impossibility to perform SK maneuvers at eclipse epochs, minimum elapsed time between two consecutive firing, on-off profile of the thrusters, thrust allocation). Indeed, the bigger specific impulse of electric thrusters leads to consequent savings in fuel consumption and a reduction of the satellite mass, enabling increased payload capacity and improved satellite lifetime.

Considering these technological and operational features, optimal control strategies for electric SK taking various constraints into account have to be carefully designed. The problem of fuel-optimal station keeping is in general expressed as an optimal control problem (OCP). Several types of techniques for solving the resulting optimal sta-tion keeping control problem exist. When the thrust is considered as impulsive or when simple models are used to describe the disturbing forces, analytical solutions provide control laws, as in [START_REF] Sukhanov | On one approach to the optimization of low-thrust station keeping manoeuvres[END_REF]. Otherwise, it is necessary to resort to numerical methods, such as direct collocation based methods as described in [START_REF] Elnagar | The Pseudospectral Legendre Method for Discretizing Optimal Control Problem[END_REF][START_REF] Hargraves | Direct Trajectory Optimization using nonlinear programming and collocation.pdf[END_REF][START_REF] Hull | Conversion of optimal control problems into parameter optimization problems[END_REF][START_REF] Betts | Survey of Numerical Methods for Trajectory Optimization[END_REF][START_REF] Betts | Mesh refinement in direct transcription methods for optimal control[END_REF][START_REF] Betts | Very Low Thrust Trajectory Optimization Using a Direct SQP Method[END_REF]. For this family of approaches, the state and the control variables are discretized to transcript the original optimal control problem into a non linear programming problem. Conversely, indirect methods rely on the derivation of necessary optimality conditions expressed via the socalled Pontryagin Maximum Principle (PMP) and on the numerical solution of the two-points boundary value problem obtained from these optimality conditions. In order to counteract uncertainties affecting the control law, the use of Model Predictive Control algorithms is proposed in [START_REF] Weiss | Opportunities and Potential of Model Predictive Control for Low-Thrust Spacecraft Station-Keeping and Momentum-Management[END_REF][START_REF] Weiss | Model Predictive Control for Simultaneous Station Keeping and Momentum Management of Low-Thrust Satellites[END_REF]. To deal with on-off models of thrusts, the references [START_REF] Vazquez | Trajectory Planning for Spacecraft Rendezvous in elliptical orbits with On / Off Thrusters[END_REF][START_REF] Vazquez | Pulse-Width Predictive Control for LTV Systems with Application to Spacecraft Rendezvous[END_REF] use the Pulse Width Modulation technique to generate rectangular profiles from a continuous one. [START_REF] Losa | Electric Station Keeping of Geostationary Satellites: a Differential Inclusion Approach[END_REF] 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 [START_REF] Losa | Station Keeping of Geostationary Satellites with On-Off Electric Thrusters[END_REF].

In the reference [START_REF] Gazzino | Optimal Control for Minimum-Fuel Geostationary Station Keeping of Satellites Equipped with Electric Propulsion[END_REF], the two first steps of a decomposition technique for solving the station keeping control problem is presented. For the first step, an indirect method based on the application of the PMP with mixed controlstate constraints is applied to solve a simplified optimal SK control problem, without considering the hard constraints on the control law (thrust constraints such as latency between two bursts of the same thruster and no simultaneous thrusting for instance). In a second step, a numerical approach is used to enforce all the thrust constraints left apart at the first step, thanks to dedicated equivalence schemes. If the operational constraints are respected, the spacecraft trajectory is composed of thrusting arcs separated by coasting arcs. Thus, the control profile switches from a time interval where all thrusters are off to a time interval for which one thruster is on, and vice-versa. Therefore, the system can be viewed as a switched system composed by one subsystem per thruster and one subsystem describing the dynamics during coasting arcs. The reference [START_REF] Gazzino | Solving the Minimum-Fuel Low-Thrust Geostationary Station Keeping Problem via the Switching Systems Theory[END_REF] takes advantage of the method proposed by [START_REF] Xu | Optimal control of switched systems based on parameterization of the switching instants[END_REF] consisting in computing the optimal switching times of switched systems thanks to a time change of coordinates. The idea of the proposed paper is hence to demonstrate the benefit of solving the station keeping problem with a three-step method. The two first steps are the ones of [START_REF] Gazzino | Optimal Control for Minimum-Fuel Geostationary Station Keeping of Satellites Equipped with Electric Propulsion[END_REF] and the third step optimizing the switching times between coasting arcs and actuated arcs highly improves the fuel consumption. Therefore, the proposed contribution describes the three-step decomposition method in a unified framework, allowing to illustrate the benefit of this three-step method on a real example given by our industrial partner and aerospace manufacturer Thales Alenia Space.

Problem statement

Geostationary Spacecraft Dynamics

The motion of a spacecraft orbiting the Earth on a Geostationary Earth Orbit (GEO) can be described with the equinoctial orbital elements as defined in [START_REF] Battin | An Introduction to the Mathematics and Methods of Astrodynamics[END_REF]:

x eoe = a e x e y i x i y M Θ T ∈ R 6 , ( 1
)
where a is the semi-major axis, (e x , e y ) the eccentricity vector, (i x , i y ) the inclination vector, M Θ = ω+Ω+M -Θ the mean longitude with Ω the right ascension of the ascending node, ω the perigee's argument, M the mean anomaly and Θ(t) the right ascension of the Greenwich meridian. The dynamics of the spacecraft may be represented by the following non linear state-space model:

dx eoe dt = f L (x eoe , t) + f G (x eoe , t)u. ( 2
)
f L ∈ R 6 is the Lagrange contribution part of the external disturbing forces. For a GEO spacecraft, the main disturbing forces are the non-spherical part of the Earth gravitational potential that mainly affects the mean longitude M Θ , the solar radiation pressure (SRP) that mainly affects the eccentricity vector (e x , e y ) and the gravitational attraction of the sun and the moon that mainly affects the inclination vector (i x , i y ) (see [START_REF] Shrivastava | Orbital Perturbations and Stationkeeping of Communication Satellites[END_REF]). These forces are described in detail by the CNES ORANGE model (cf. [START_REF] Campan | ORANGE: Orbital analytical model for geosynchronous satellite[END_REF]) or in the reference [START_REF] Soop | Handbook of Geostationary Orbits[END_REF]. f G ∈ R 6×3 is the Gauss contribution part for the disturbing forces that do not derive from a potential. In the case of a GEO spacecraft, this contribution consists in the acceleration produced by the thrusters that can be viewed as a disturbing one. u = u R u T u N T ∈ R 3 is the control vector expressed in the local orbital RT N frame (also written RSW ) defined in [START_REF] Soop | Handbook of Geostationary Orbits[END_REF] by (see Figure 1):

• N is the unit vector along the kinetic momentum;

• R is the unit vector along the direction Earth's center -satellite;

• T completes the right-handed orthogonal direct basis.

In order to deal with the station keeping problem, the station keeping point is defined as:

x sk = a sk 0 0 0 0 M Θ sk T (3)
where a sk is the geosynchronous semi-major axis and M Θ sk is the station mean longitude.

The relative dynamics equations are developed by linearization of Equation (2) at the station keeping point (3) evolving on an ideal circular geostationary keplerian orbit. By denoting x = x eoe -x sk , the relative state model for the SK problem reads: where the matrices A ∈ R 6×6 , B ∈ R 6×3 , C ∈ R 3×6 and D ∈ R 6 are defined as follows:

dx dt = A(t)x + D(t) + B(t)u, ( 4
)
X G Y G Z G Equatorial Orbital Plane Ω ω ν i G Plane R T N S r
A(t) = ∂ ∂x eoe f L (x eoe (t), t) xeoe=x sk , ( 5
)
B(t) = f G (x sk , t), (6)
D(t) = f L (x sk , t). ( 7
)
The geographical coordinates of the satellite:

y eoe = T (x eoe , t)x eoe , ( 8
)
are of interest because the station keeping problem consists in constraining them in the vicinity of the station keeping geographical position y sk = r sk 0 λ sk t where r sk is the geosynchronous radius and λ sk is the station keeping geographical longitude. The relative geographical position with respect to the station-keeping position is denoted by:

y = y eoe -y sk = T (x sk , t)x = C(t)x, (9)
and is obtained by linearizing Equation ( 8) with respect to the station keeping geographical position. After the linearization, the station keeping requirements on the latitude and the longitude of the spacecraft are expressed by constraining the relative geographical position with respect to the centre of the SK window: ∀t ∈ [t 0 , t f ],

[0 1 0]C(t)x(t)
δ and [0 0 1]C(t)x(t) δ, [START_REF] Molitor | Ion Propulsion System for Stationary-Satellite Control[END_REF] with δ being the half width of the SK window in the latitude and longitude directions and [t 0 , t f ] the optimisation horizon. A trajectory satisfying the SK constraints ( 10) is called a SK-feasible trajectory (see Figure 2).

Electric Propulsion System

The considered satellite is equipped with four electric thrusters mounted on the anti-nadir face, each of them having an orientation defined by a cant angle θ and a slew angle α. These angles define the North-East, the South-East, the North-West and the South-West directions of thrust. The satellite dynamics can be stated considering the four thrusts provided by the four engines as control variables. The control u(t) expressed in the local orbital frame is a linear combination of the four thrusts such that u = ΓF , where Γ

X G Y G Z G δ δ λ sk Greenwich Meridian Θ(t) y sk δϕ δλ
= [Γ 1 | Γ 2 | Γ 3 | Γ 4 ] ∈ R 3×4 and F = [F 1 F 2 F 3 F 4 ] T ∈ [0, F max ] 4 .
The thrust direction vectors Γ j ∈ R 3 are given by: Γ j = 1 m -sin θ j cos α j -sin θ j sin α j -cos θ j T , (11) where the cant angles θ j and slew angles α j are defined exactly as in [START_REF] Anzel | Method and apparatus for a Satellite Station Keeping[END_REF] (see Figure 3). The mass m of the spacecraft is supposed to be constant. For the sake of simplicity, the thrust vector is normalized by the maximum level of thrust F max . It is thus possible to write F = F max F with F ∈ [0, 1] 4 . As the thrust profile is made of a series of on-off thrusts, the rectangular signal of the j th thrust of the thruster i 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 (see Figure 4). The control function is thus rewritten as:

F (t) = F t, {t i j , ∆t i j } i=1,...,4 j=1,...,Pi , =              P1 j=1 1 [t 1 j -∆t 1 j ,t 1 j +∆t 1 j ] (t) P2 j=1 1 [t 2 j -∆t 2 j ,t 2 j +∆t 2 j ] (t) P3 j=1 1 [t 3 j -∆t 3 j ,t 3 j +∆t 3 j ] (t) P4 j=1 1 [t 4 j -∆t 4 j ,t 4 j +∆t 4 j ] (t)              . ( 12
)
where P i is the number of thrusts of thruster i. Technological operational constraints inherent to the proposed actuation system are defined by: (i) thrusters cannot have simultaneous thrusts;

∆t i,j t i,j t Fi 1
(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 by 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 .

The restricted power available on-board prevents two thrusters from being active simultaneously and the minimum time latency between two thruster firings is imposed in order to allow an efficient battery recharge (see [START_REF] Losa | Station Keeping of Geostationary Satellites with On-Off Electric Thrusters[END_REF]).

The constraint for the time latency between the thrust k of thruster i and the thrust l of thruster j is mathematically expressed as:

|t i,k -t j,l | -(∆t i,k + ∆t j,l ) K i,j , (13)
for k = 1 . . . P i and l = 1 . . . P j , where K i,j = T s if i = j (constraint (iii)) and K i,j = T d otherwise (constraint (iv)). In addition, some other convenient constraints force the thrusters to be active only during the resolution time interval:

t i,j -∆t i,j t 0 and t i,j + ∆t i,j t f . (

A trajectory satisfying the operational constraints for the considered electric propulsion system is called operationally feasible.

Fuel-optimal Station Keeping Control Problem

The main goal of the station keeping system is to maintain the longitude and the latitude of the satellite in a box defined by its size δ by acting on the orbital parameters via the 4 thrusters while reducing the fuel consumption to extend the operational lifetime of the satellite. The associated optimal control problem is in general defined over a fixed horizon for the computation of optimal control laws. The performance index to minimize in a Station Keeping operation is naturally defined as:

J = t f t 0 4 thruster i=1 P i j=1 |uR ij (t)| + |uT ij (t)| + |uN ij (t)| dt (15) = Fmax t f t 0 4 thruster i=1 ||Γi||1 P i j=1 | Fi,j(t)|dt (16)
= 2Fmax

4 thruster i=1 ||Γi||1 P i j=1 ∆ti,j, . (17)
Removing the constant multiplicative terms that does not affect the optimal solution, the objective function may be reduced to:

J = 4 i=1 Pi j=1 ∆t i,j . ( 18
)
Considering all the constraints described above, the minimum-fuel SK problem to solve may be summarized as the following optimal control problem: Problem 1. Find the sequence of dates {t i,j } and durations {∆t i,j }, for i = 1 . . . 4, j = 1 . . . P i solutions of the minimization problem:

min ti,j ,∆ti,j J = 4 i=1 Pi j=1 ∆t i,j , s.t.                        ẋ(t) = A(t) x(t) + D(t) + B(t) F t, {t i j , ∆t i j } i=1,...,4 j=1,...,Pi , x(t 0 ) = 0, |[0 1 0]C(t)x(t)| δ, |[0 0 1]C(t)x(t)| δ, 2∆t i,j T l , t i,j -∆t i,j t 0 , t i,j + ∆t i,j t f , |t i,k -t j,l | -(∆t i,k + ∆t j,l ) K i,j , ( 19
) with B(t) = F max ΓB(t).
•

The Problem 1 raises some difficult mathematical issues. Firstly, due to the parametrization of the rectangular functions by t i,j and ∆t i,j , the dynamics define transcendental equations to be satisfied for the problem to solve. Secondly, the state constraints and the logical actuation constraints make the previous optimal control problem difficult to handle with existing methods. In addition, if the optimal number of thrusts per thruster is not known a priori, the problem to solve is a very challenging nonlinear mixed integer programming problem. For all these reasons, traditional solvers often fail to find a solution to the Problem 1.

The proposed paper aims thus at finding an efficient way to solve the fuel optimal station keeping problem by splitting the overall problem in a sequence of subproblems. The proposed resolution method consists in three different steps:

(A) The thrusters constraints ( 13) and ( 14) 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 initialized by a direct method solution dedicated to the search of adjoint variables as described in [START_REF] Grimm | Adjoint estimation from a direct multiple shooting method[END_REF] and [START_REF] Stryk | Direct and indirect methods for trajectory optimization[END_REF]. The PMP is applied to derive the first order necessary optimality conditions and the associated Two Points Boundary Value Problem (TP-BVP). This TPBVP is then solved via a collocation method. In order to recover an on-off profile from the continuous solution obtained, a threshold parameter ς = ς 1 ς 2 ς 3 ς 4 ∈ R 4 has to be chosen. At the end of this step, the state trajectory fulfills the SK requirements but may not be operationally feasible.

(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 solutions compatible with actuation constraints (i)-(iv) and ( 14). The resulting state trajectory is operationally feasible but not necessarily SK-feasible.

(C) The solution of step (B) being very sensitive to the threshold parameter ς chosen at step (A) the proposed third part optimizes the switching times for each thruster with the objective to reduce the fuel consumption. This step is based on an optimization technique stemming from the switched systems theory (see [START_REF] Xu | Optimal control of switched systems based on parameterization of the switching instants[END_REF]) and takes advantage of the decomposition of the overall systems into several subsystems: one for each thruster being on and one when all the thrusters are off. The ensuing trajectory is both operationally and SK-feasible.

To sum up, instead of solving the overall SK problem at once, this problem is split into three smaller and simpler subproblems that are easier to solved in general. Section 3 will be devoted to the first part described above: solution of the simplified OCP via the hybrid method. Then in Section 4, the thrusters constraints are enforced via the resolution of a second minimization problem called equivalence scheme. In Section 5, the switching times for each thrust are optimized using the switched systems theory, noticing that a system with a bang-bang control profile can be considered as being made up of two subsystems switching from one with control on to one with control off, and vice-versa. A commutation between a coasting arc and a fired arc is therefore interpreted as a commutation between the two subsystems.

Solution of the simplified OCP (step (A))

Removing the operational actuation constraints ( 13) and ( 14) allows to simplify the original optimal control problem. In particular, the thrust functions Fi are not modeled a priori as rectangular functions parametrized by t i,j and ∆t i,j and the simplified OCP to be solved reads as : Problem 2. Find the functions t → Fi (t), i = 1, . . . , 4, solutions of the minimisation problem: min

F (t)∈[0;1] 4 J = t f t0 4 i=1 Fi (t)dt, (20)
with the constraints:

     ẋ(t) = A(t) x(t) + D(t) + B(t) F (t), |[0 1 0]C(t)x(t)| δ, |[0 0 1]C(t)x(t)| δ, x(t 0 ) = 0, x(t f ) = 0, (21)
• Note that the constraint on the final state vector is not strictly necessary here from a theoretical point of view. However, this additional constraint may be heuristically justified by the numerous numerical experiments for which it proved to be helpful in reducing the sensitivity of the solution of Step A with respect to tuning parameters by preventing the longitude to end up on the boundary of the SK window.

Problem 2 is a minimum-fuel linear OCP defined on a fixed horizon with constraints both on the state and the control vectors. This problem is solved with the PMP.

The state geographical constraints are handled through the following additional state variable whose dynamics is given by (see the reference [START_REF] Naidu | Optimal Control Systems[END_REF] for instance): [START_REF] Sukhanov | On one approach to the optimization of low-thrust station keeping manoeuvres[END_REF] with :

ẋ7 (t) = ψ(t, x(t)) = 4 i=1 1 2 ψ 2 i 1 + sign(ψ i (t, x(t))) ,
ψ 1 (t, x(t)) = C 2 (t)x(t) -δ, ( 23a
)
ψ 2 (t, x(t)) = -C 2 (t)x(t) -δ, ( 23b
)
ψ 3 (t, x(t)) = C 3 (t)x(t) -δ, ( 23c
)
ψ 4 (t, x(t)) = -C 3 (t)x(t) -δ, ( 23d
)
where C l (t) is the l th line of the matrix C(t).

The Hamiltonian of the system is:

H(x(t), F (t), λ(t)) = 4 i=1 Fi (t) + µ 1 ψ(x(t)) + λ(t) t A(t) x(t) + D(t) + B(t) F (t) , ( 24
)
where λ(t) ∈ R 6 is the adjoint vector, and µ 1 is a constant parameter. Choosing the value of this parameter allows to give more importance for minimizing the fuel consumption or enforcing of the SK constraints.

Applying the PMP described in [START_REF] Bonnard | Mécanique céleste et contrôle des véhicules spatiaux[END_REF] and [START_REF] Naidu | Optimal Control Systems[END_REF], the minimization condition:

H x * (t), F * (t), λ * (t) = min u∈[0,1] 4 H (x * (t), u(t), λ * (t)) (25
) leads to the following first order necessary switching condition for the optimal control:

∀i = 1, . . . , 4, F * i (t) = 1 if 1 + [λ * T (t) B(t)] i < 0, 0 if 1 + [λ * T (t) B(t)] i > 0, (26) where [
h] i stands for the i th component of a generic vector h. The previous necessary condition can be rewritten using the switching function:

v → f s (v) = 1 2 (1 -sign v 1 ) . . . (1 -sign v 4 ) T , ( 27
)
leading to the optimal control:

F * (t) = f s 1 + λ T (t) B(t) , ( 28
)
where 1 is a vector of ones.

The switching function:

f s (1 + λ T B) = 1 2 1 -sign 1 + λ T (t) B(t) (29)
defines a non continuous piecewise constant control profile. As the TPBV Problem 3 is solved with the bvp4c function of Matlab that requires continuous functions, the method described by [START_REF] Mazal | Optimal On-Off Cooperative Maneuvers for Long-term Satellite Cluster Flight[END_REF] is used in order to approximate the discontinuous control profile into a continuous one. f s is thus approximated by:

f s (v) ≈ fs (v, M ) = 1 - 1 π arctan (M v) + π 2 , ( 30
)
where M is a parameter controlling the slope of the approximating switching function.

The two-point boundary value problem is built from the canonical equations of Hamilton and the transversality conditions and is defined as Problem 3 Problem 3. Find the functions x and λ solutions of the minimization problem:

ẋ(t) = A(t)x(t) + D(t) + B(t) fs (1 + λ T B, M ), λ(t) = -A(t) T λ(t) -µ 1 dψ dx (x(t)),
x(t 0 ) = 0, x(t f ) = 0, λ(t 0 ) and λ(t f ) free.

(31)

•

This TPBVP is initialized by a direct collocation method, as described in the references [START_REF] Bulirsch | Combining direct and indirect methods in optimal control: Range maximization of a hang glider[END_REF] and [START_REF] Bonnard | Optimal control with state constraints and the space shuttle re-entry problem[END_REF] for instance, and then solved with an homotopy method in order to gradually increase the slope of the switching function by increasing accordingly the value of the parameter M of the approximated switching function.

A threshold parameter ς = ς 1 ς 2 ς 3 ς 4 T ∈ R 4

transforms the continuous profile into a discontinuous one as the Figure 5 shows. The decision variables, the middle times of the thrusts and their half-width durations are recovered. This transformation permits to extract from the on-off profile the number of thrusts for each thrusters P i , i = 1, . . . , 4. According to Figure 5, it is clear that the value of the parameter ς has a direct influence on the associated half width thrust duration. As the solution of Problem 3 does not satisfy in general the operational actuation constraints described in Section 2.2, the next section proposes an additional step for which an auxiliary problem enforces the actuation constraints on a new and equivalent control law, while preserving the structure and the overall effect of the thrust.

t F i ς b ς a ∆t i,j;a ∆t i,j;b t i,j

Enforcing the operational thrusters constraints (step (B))

The solution of the TPBV Problem 2 consists in a series of rectangle signals naturally resulting from the application of the PMP. However, this control law does not respect the actuation constraints left apart in the step (A). For instance, simultaneous activation of two different thrusters may occur. The aim of this second part is to find an equivalent control law ensuring that the resulting trajectory will be operationally feasible. 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. Let FBV P be the control obtained by solving Problem 3. Finding a Consumption Based Equivalent (CBE) control for the satellite is then equivalent to solve Problem 4 defined as follows: Problem 4. Find the t i,j and the ∆t i,j , i = 1, . . . , 4, j = 1, . . . , P i minimising: min ti,j ,∆ti,j

4 i=1 FBVP,i (t) 1 -2 Pi j=1 ∆t i,j , (32)
subject to the constraints:

     2∆t i,j T l , |t i,k -t j,l | -(∆t i,k + ∆t j,l ) K i,j , t i,j -∆t i,j 0, t i,j + ∆t i,j T. ( 33
)
• This problem is a non linear optimization problem where FBVP,i (t) 1 is the L 1 norm of the i th component of the solution of Problem 3.

The second way to obtain an equivalent control respecting the actuation constraints is to define an "effect-based" equivalent (EBE) control. As the dynamics equation is given by Equation (4), the state vector at time t is given by: [START_REF] Gazzino | Solving the Minimum-Fuel Low-Thrust Geostationary Station Keeping Problem via the Switching Systems Theory[END_REF] where Φ(t, t 0 ) is the state transition matrix, implicitly defined by the homogeneous differential equation Φ(t, t 0 ) = A(t)Φ(t, t 0 ), Φ(t 0 , t 0 ) = Id (see [START_REF] Antsaklis | Linear Systems[END_REF]).

x(t) = Φ(t, t 0 )x(t 0 ) + t t0 Φ(t, τ ) [D(τ ) + B(τ )Γu(τ )] dτ,
To get a control profile that has the same effect at time t f as the solution of Problem 3 and that respects the actuation constraints, the Problem 5 defined as follows has to be solved: Problem 5.

min t i,j ,∆t i,j 6 i=1 t f t 0 Φ(T, τ ) B(τ ) FBVP (τ ) i dτ - 4 k=1 P k j=1 ∆t k,j t k,j +∆t k,j t k,j -∆t k,j [Φ(T, τ ) B(τ )]idτ 2 , ( 35
)
such that the constraints (33) are satisfied.

•

Both Problems 4 and 5 can be solved by classical non linear optimization solvers.

The state trajectories obtained from the control profiles solutions of Problems 4 and 5 are very sensitive to the threshold parameter ς chosen at the end of the step (A). For some values of ς, the computing trajectory may be operationally feasible but not restricted to the SK window. The benefit of the step (B) is to provide a sequence of thrusting/coasting arcs admissible for the operational constraints (i)-(iv). This sequence is described by {U } = {U 1 , . . . , U K }, with:

U k ∈            0 0 0 0     ,     1 0 0 0     ,     0 1 0 0     ,     0 0 1 0     ,     0 0 0 1            , ( 36
)
assuming that there are K intervals over which the thrust is constant. The middle time of the thrust t i,j and the half width durations of the thrusts ∆t i,j are ranked in the increasing thrust middle times order t 1 , . . . , t P , assuming that there are P thrusts, and transformed to the commutation times as:

s 2k-1 = t k -∆t k , s 2k = t k + ∆t k , , k = 1, . . . , P. ( 37
)
In this case, the following equality holds : K = 2P + 1.

In order to guarantee both the operationally and the SK-feasibility of the state trajectory, a step (C) relying on the switched systems theory with the results of step (B) as initialization is used to enforce all the constraints while reducing the consumption.

Optimization of the Commutation Times with the Switched Systems Theory (step (C))

The reference [START_REF] Xu | Optimal control of switched systems based on parameterization of the switching instants[END_REF] presents a technique to optimize the commutation times between the different subsystems. This method resorts to a time change of coordinates and a parametrization of the commutation times. The control profile computed at the end of the step (B) is an on-off control profile. The idea of this step is to notice that a bang-bang profile is a sequence of commutations between a system for which the control is on and another system for which the control is off. Hence, the date when the thruster switches on or switches off can be seen as a commutation between the subsystem for which the control is on to the subsystem for which the control is off, or viceversa. This remark applies to the four thrusters of the considered satellite. Hence the method proposed by [START_REF] Xu | Optimal control of switched systems based on parameterization of the switching instants[END_REF] will be applied in order to optimize the commutations between the coasting and the firing arcs. To this end, it is mandatory to assume that the order of the sequence of the active subsystems is an operationally feasible one. As the control profile computed at the end of step (B) is operationally feasible, this firing sequence can be used. Seeking these on and off times is equivalent to find the optimal commutation times between the subsystems.

Noting that each firing arc is separated by a coasting arc, the odd commutation times thus correspond to the beginning of a thrust and the even ones to the ending of a thrust. The notation of the commutation times is extended to s 0 = t 0 and s 2P +1 = t f . Therefore, minimizing the fuel consumption is equivalent to minimize the length of the intervals on which the thrusters are on. Hence the objective function is written as:

J C ({s k }) = F max P i=1 (s 2i -s 2i-1 ). ( 38
)
The state constraints are handled by mean of the penalization term described by Equation ( 22) added to the objective function. The new objective function for the step (C) of the decomposition method reads thus: [START_REF] Grimm | Adjoint estimation from a direct multiple shooting method[END_REF] with µ 3 a parameter that balance the objective function between the fuel consumption minimization and the stay in the SK window. Note that µ 3 may be different from µ 1 defined in Equation [START_REF] Hargraves | Direct Trajectory Optimization using nonlinear programming and collocation.pdf[END_REF]. The tuning of this parameter may be a difficult issue to settle since selecting µ 3 too large results in an overconsumption while a too small µ 3 leads to an infeasible trajectory.

JC (x, {s k }) = F max P i=1 (s 2i -s 2i-1 ) + µ 3 2P +1 k=1 (s k -s k-1 ) k k-1 ψ(t(τ ), x)dτ,
As in the reference [START_REF] Xu | Optimal control of switched systems based on parameterization of the switching instants[END_REF], the time change of coordinates is used to parametrize the switching times:

t = s k + ∆ k (τ -k) if t ∈ [s k , s k+1 ], (40)
with ∆ k = s k+1 -s k . The system dynamics is rewritten with the new time variable:

dx(τ ) dτ =                                ∆ 2k-2 A(t(τ ))x(τ ) + D(t(τ )) if τ ∈ [2k -2, 2k -1], ∆ 2k-1 A(t(τ ))x(τ ) + D(t(τ )) + B(t(τ ))ΓF max U k if τ ∈ [2k -1, 2k], . . . ∆ 2P A(t(τ ))x(τ ) + D(t(τ )) if τ ∈ [2P, 2P + 1], (41)
for k = 1, . . . , P . The state vector can now be considered as a function of the new time variable τ and of the switching times s k : x = x(τ, {s k }).

The station keeping OCP reads thus in the switching systems framework: Problem 6. Find the optimal, switching sequence {s k }, k ∈ {1, . . . , 2P }, with P fixed, of the min-imisation problem:

min x,{s k } J C (x, {s k }) = F max P i=1 (s 2i -s 2i-1 ) +µ 3 2P +1 k=1 (s k -s k-1 ) k k-1 ψ(t(τ ), x)dτ, s. t.                            ẋ(τ ) =              ∆ 2k-2 A(t(τ ))x(τ ) + D(t(τ )) if τ ∈ [2k -2, 2k -1], ∆ 2k-1 A(t(τ ))x(τ ) + D(t(τ )) +B(t(τ ))ΓF max U k if τ ∈ [2k -1, 2k] x(0) = 0, s 2k -s 2k-1 T l , s 2k+1 -s 2k T g , ( 42
) with T g = T s for the firing arcs of the same thruster and T g = T d for the firing arcs of two different thrusters.

•

As the switching times determine the structure of the control profile, the state trajectory is recovered by propagation of the system dynamics from the initial condition. Problem 6 can be solved with a descent algorithm. To this end, it is necessary to compute the derivative of the performance index in order to obtain the descent direction. The computation of this derivative requires to distinguish the odd and the even commutation times:

dJ C ds 2l-1 = -F max + µ 3 2l-1 2l-2 ψ(t(τ ), x(τ )dτ -µ 3 2l 2l-1 ψ(t(τ ), x(τ ))dτ + µ 3 (s 2l-1 -s 2l-2 ) 2l-1 2l-2 (τ -2l + 2) ∂ψ ∂t (t(τ ), x(τ )) + ∂ψ ∂x (t(τ ), x(τ )) ∂x ∂s 2l-1 (τ ) dτ + µ 3 (s 2l -s 2l-1 ) 2l 2l-1 (-τ + 2l + 2) ∂ψ ∂t (t(τ ), x(τ )) + ∂ψ ∂x (t(τ ), x(τ )) ∂x ∂s 2l-1 (τ ) dτ, (43)
and :

dJ C ds 2l = F max -µ 3 2l+1 2l ψ(t(τ ), x(τ )dτ + µ 3 2l 2l-1 ψ(t(τ ), x(τ ))dτ + µ 3 (s 2l+1 -s 2l ) 2l+1 2l (-τ + 2l + 1) ∂ψ ∂t (t(τ ), x(τ )) + ∂ψ ∂x (t(τ ), x(τ )) ∂x ∂s 2l-1 (τ ) dτ + µ 3 (s 2l -s 2l-1 ) 2l 2l-1 (τ -2l + 1) ∂ψ ∂t (t(τ ), x(τ )) + ∂ψ ∂x (t(τ ), x(τ )) ∂x ∂s 2l-1 (τ ) dτ. (44)
As the switching times {s k } are independent of the time variable τ , the derivatives of the state vector with respect to the commutation times verify the following relation:

d dτ ∂x(τ ) ∂s k = ∂ ∂s k dx(τ ) dτ . ( 45
)
It is thus possible to derive a dynamic equation for ∂x(τ ) ∂t k by differentiation of the dynamics equation of the state vector x(τ ) with respect to the commutation time s k .

Hence the function τ → ∂x(τ ) ∂s k is obtained by integrating the following dynamics equations:

d dτ ∂x(τ ) ∂s 2l-1 =                                                                          •(s 2k-1 -s 2k-2 )A(t(τ )) ∂x(τ ) ∂s 2l-1 if τ ∈ [2k -2, 2k -1], k = l, •(s 2k -s 2k-1 )A(t(τ )) ∂x(τ ) ∂s 2l-1 if τ ∈ [2k -1, 2k], k = l, •A(t(τ ))x(τ ) + D(t(τ )) + (s 2l-1 -s 2l-2 ) (τ -2l + 2) ∂A(t) ∂t x(τ )+ ∂D(t) ∂t + A(t(τ )) ∂x(τ ) ∂s 2l-1 if τ ∈ [2l -2, 2l -1], • -A(t(τ ))x(τ ) + D(t(τ )) + B(t(τ ))U l + (s 2l -s 2l-1 ) (-τ + 2l + 2) ∂A(t) ∂t x(τ ) + ∂D(t) ∂t + ∂B(t) ∂t U l + A(t(τ )) ∂x(τ ) ∂s 2l-1 if τ ∈ [2l -1, 2l], (46) and d
dτ ∂x(τ ) ∂s 2l =                                                                        •(s 2k-1 -s 2k-2 )A(t(τ )) ∂x(τ ) ∂s 2l if τ ∈ [2k -2, 2k -1], k = l, •(s 2k -s 2k-1 )A(t(τ )) ∂x(τ ) ∂s 2l if τ ∈ [2k -1, 2k], k = l, • -[A(t(τ ))x(τ ) + D(t(τ ))] + (s 2l+1 -s 2l ) (-τ + 2l + 1) ∂A(t) ∂t x(τ )+ ∂D(t) ∂t + A(t(τ )) ∂x(τ ) ∂s 2l if τ ∈ [2l, 2l + 1], •A(t(τ ))x(τ ) + D(t(τ )) + B(t(τ ))U l + (s 2l -s 2l-1 ) (τ -2l + 1) ∂A(t) ∂t x(τ ) + ∂D(t) ∂t + ∂B(t) ∂t U l + A(t(τ )) ∂x(τ ) ∂s 2l if τ ∈ [2l -1, 2l]. ( 47
)
The derivative of the penalization function with respect to the state vector is:

∂ψ ∂x (t, x(τ )) = 2C 2 (t) [C 2 (t)x(t) -δ] sign C 2 (t)x(t) -δ + 1 + 2C 2 (t) [C 2 (t)x(t) + δ] sign -C 2 (t)x(t) -δ + 1 + 2C 3 (t) [C 3 (t)x(t) -δ] sign C 3 (t)x(t) -δ + 1 + 2C 3 (t) [C 3 (t)x(t) + δ] sign -C 3 (t)x(t) -δ + 1 , (48)
and its derivative with respect to time is:

∂ψ ∂t (t, x(τ )) = 2 dC 2 (t) dt x [C 2 (t)x(t) -δ] sign C 2 (t)x(t) -δ + 1 + 2 dC 2 (t) dt x [C 2 (t)x(t) + δ] sign -C 2 (t)x(t) -δ + 1 + 2 dC 3 (t) dt x [C 3 (t)x(t) -δ] sign C 3 (t)x(t) -δ + 1 + 2 dC 3 (t) dt x [C 3 (t)x(t) + δ] sign -C 3 (t)x(t) -δ + 1 .
(49) The ensuing nonlinear non convex optimization problem is solved by an interior point algorithm based on [START_REF] Byrd | An interior point algorithm for large-scale nonlinear programming[END_REF] and encapsulated in the function fmincon in MATLAB©.

Putting all the three steps (A), (B) and (C) together, the SK problem can therefore be solved according to Algorithm 1.

Numerical Results

The proposed three-step decomposition method is applied on a realistic example and its benefit on the reduction of the fuel consumption while enforcing both the operational and SK constraints is demonstrated.

The considered satellite weights 4850 kg and is equipped with the 4 electric thrusters propulsion system described in Section 2.2 with thrusters in the North-East, North-West, South-East and South-West directions. This satellite has to be controlled in order to remain close to its geostationary position at a fixed longitude λ sk and a fixed latitude ϕ sk = 0. The station keeping is performed over the time horizon defined by t 0 = 0 and t f = 1 week using the threestep technique described in this paper.

Solution of step (A)

The slope parameter M for the approximated switching function is arbitrarily chosen to be increased from 10 0 to 10 4 and from trial and error numerical simulations, it appears that the best choice for the tuning parameter µ 1 is µ 1 = 1.10 4 . On the Figure 6 depicting the continuous control profile, it is possible to notice that the operational constraints are not respected, as the thrusters 1 and 3 have simultaneous firings as well as the thrusters 2 and 4. Therefore, the step (A) is not sufficient to solve the whole Algorithm 1 Solve the minimum fuel SK problem Require: x(t 0 ) the initial position Ensure: the trajectory is operationally and SK-feasible {Step (A): SK problem without the operational con-straints} 1: Remove the operational constraints (i)-(iv) from Problem 1 2: Fi (t) i=1,...,4,t∈[t0,t f ] ← the optimal control for the SKfeasible trajectory computed by an indirect method initialized by a direct method (see the reference [START_REF] Gazzino | Optimal Control for Minimum-Fuel Geostationary Station Keeping of Satellites Equipped with Electric Propulsion[END_REF]) 1 shows that these two different values of the threshold parameter lead to different number of thrusts per thrusters. The two considered threshold parameters are ς a = 0.8 0.8 0.8 0.8 T and ς b = 0.9 0.9 0.9 0.9 T . The threshold values depend on the configuration of the satellite and on the mission characteristics. Finely tuning these parameters is pivotal as the number of thrusts will remain fixed for the steps (B) and (C). In this study, these proposed values have been evaluated by proceeding a linear search. It is possible to find in the literature some rules of thumb for geostationary SK strategies, in particular in [START_REF] Campan | Les techniques de maintien à poste de satellites géostationaires[END_REF], [START_REF] Soop | Handbook of Geostationary Orbits[END_REF] and [START_REF] Sidi | Spacecraft Dynamics and Control[END_REF]. In these references, the effect of the perturbing forces on the orbit is studied, and some general SK laws are derived. The North-South effect of the Sun and the Moon attractions are the most pregnant forces that must be corrected each half-orbit, once in the North direction and half an orbit later in the South direction. The East-West drift is meanwhile corrected by setting different thrust durations for each thrusters. This SK strategy was successfully used in an industrial context in [START_REF] Anzel | Method and apparatus for a Satellite Station Keeping[END_REF]. The control profile of Figure 6 shows that the two South thrusters (F 2 , F 4 ) have a thrust at the beginning of each day, and the two North thrusters (F 1 , F 3 ) have a thrust half a day later. On the control profile, the thrusts have always different durations in order to compensate for the East-West drift. The physical rule of thumb can thus be recovered from a systematic optimization process.

Solution of step (B)

The trajectory obtained at the end of step (A) is SKfeasible but not operationally feasible. Therefore, the step (B) is mandatory in order to enforce the operational constraints. As described in the Section 4, the recovery of the middle time of the thrust and their half duration are very sensitive to the threshold parameter ς. Examples are given and some of them exhibit an operationally feasible but not SK-feasible state trajectory, depending of the value of ς. For ς a , the trajectory computed after the CBE and EBE schemes are operationally feasible but not SK feasible (see Figure 7). Despite the fact that some of the trajectories are not SK-feasible, the aim of this step is just to enforce the operational constraints since an operationally feasible initial solution is mandatory for the step (C), whose aim is to enforce SK constraints. -: CBE scheme with ςc, --: EBE scheme with ς d , --: station keeping window.

Solution of step (C)

The step (C) optimizes the commutation times of the control profile solution of step (B). The firing sequence from the step (B) is supposed to be the optimal one, and its commutation times have to be optimized. For both equivalence schemes of the step (B), this step manages to minimize the fuel consumption while enforcing the operational as well as the SK constraints. For the two state trajectories of Figure 8 computed with ς c and ς d , the step (C) reduces the fuel consumption, as the Table 2 shows. Figures 9 and 10 display the control profile optimized thanks to the switching system framework. On these Figures, it appears that the length of the thrusts have been reduced, implying a decrease of the consumption. Figures 11 and12 illustrate the fact that even if the trajectory resulting from step (B) is not SK-feasible, the optimization of the commutation times performed in step (C) manages to enforce the SK constraints, The Figures 13 and14 show that this step reduces the consumption as well.

For these computations, µ 3 = 10 8 . As the integration of the dynamics equation ( 41) and its time derivatives with respect to the switching times ( 46) and ( 47) are numerically integrated with large steps in order to reduce the computation time, some margins have been taken and the size of the SK window has been reduced by 10%. In Table 2, several values of the threshold parameter ς used at the end of the step (A) are given, as well as the stemmed number of thrusts P . This Table 2 shows clearly the fuel consumption reduction between the steps (B) and (C).

In the reference [START_REF] Eckstein | Optimal Station Keeping by Electric Propulsion With Thrust Operation Constraints[END_REF], the SK strategy leads to a velocity increment budget of 3.204 m/s for 10 days, which leads to 2.24 m/s for a week. In [START_REF] Soop | Handbook of Geostationary Orbits[END_REF]Table 3,p.150], it is possible to find inclination correction requirements for one year. According to the authors strategy, the annual velocity correction to be performed is a least 78.67 m/s, i.e. 1.51 m/s for a week. Therefore, the proposed three-step correction algorithm creates SK strategies that require less fuel consumption than the ones existing in the litterature.

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. In a first step, a classical optimal control problem is solved with state constraints using a precise indirect method initialized by a collocation based direct method. In a second step, two ways of dealing with the thrusters operational constraints are proposed resulting in two different fuel consumption results. In a third step, an optimization of the commutation times coming from the theory of the switched systems is used in order to optimize the on and off times of the firing sequence obtained at the end of the second step. In fact, the system can be naturally decomposed into several subsystems, one per thruster whose control is on and one corresponding to a coasting arc. As demonstrated by the numerical examples, the threshold parameter used in the first step has a high impact on the control profile computed at the second step, so that the resulting trajectory may not be SK-feasible in some cases. However, as this trajectory is operationally feasible, the third step objective is to enforce both the operational and the SK constraints. Despite these positive results, some issues remain open. The state constraints are taken into account with a penalization function, but alternative formulations of the first order necessary conditions coming from the Pontryagin Maximum Principle exist and should be considered. Furthermore, in the last step, the dynamics and its time derivative have to be numerically integrated, once per commutation time. This leads to long resolution time, that could be lowered using semi-analytic integration tools for instance.

Figure 1 :

 1 Figure 1: Geocentric inertial reference frame (G, X G , Y G , Z G ) and local orbital frame (S, R, T, N ).

Figure 2 :

 2 Figure 2: Station Keeping window.

Figure 3 :

 3 Figure 3: Orientation of the four thrusters (picture from TAS).

Figure 4 :

 4 Figure 4: Parametrization of the j th thrust of the thruster i.

Figure 5 :

 5 Figure 5: Effect of the threshold parameter ς on the half-width duration of the thrusts.

F4Figure 6 :

 6 Figure 6: Continuous control profile at the end of step (A).

  Figures 8 displays the state trajectory computed with the CBE scheme and ς c = 0.89 0.89 0.93 0.93 T and the state trajectory computed with the EBE scheme and ς d = 0.95 0.9 0.95 0.95 T . In these two cases, a fine tuning of the threshold parameter leads to operationally and SK-feasible trajectory.

Figure 7 :Figure 8 :

 78 Figure 7: Trajectories in the (λ, ϕ) plane after the step (B).-: CBE scheme with ςa, --: EBE scheme with ςa, --: station keeping window.

F4Figure 9 :

 9 Figure 9: Control profiles : -: after the step (B) with ςc and the CBE scheme (total thrust duration: 24.35 hours), --: after step (C) (total thrust duration: 72.43 min).

F4Figure 10 :

 10 Figure 10: Control profiles :-: after the step (B) with ς d and the EBE scheme (total thrust duration: 7.17 hours), --: after step (C) (total thrust duration: 26.78 min).

Figure 11 :Figure 12 :

 1112 Figure 11: Trajectories in the (λ, ϕ) plane : -: after the step (B) with ςa and the CBE scheme, --: after step (C), --: station keeping window.

F4Figure 13 :

 13 Figure 13: Control profiles : -: after the step (B) with ςa and the CBE scheme, --: after step (C).

F4Figure 14 :

 14 Figure 14: Control profiles : -: after the step (B) with ςa and the EBE scheme, --: after step (C).

  3: ς ← threshold parameter to transform the continuous control profile in a rectangle profile 4: {t i,k , ∆t i,k , {U k }} ← the middle times t i,k and the half durations ∆t i,k of the thrusts, {Step (B): equivalence schemes} 5: {t i,k , ∆t i,k } ← the middle times and half durations computed by the CBE or the EBE. {Step (C): commutation times optimization with the switching systems technique} 6: {s k } ← the commutation times computed from {t i,k , ∆t i,k }

	7: {s (0) k } ← {s k }
	8: {s k } ← the solution of Problem 6 with an interior
	point algorithm
	station keeping problem. The Table

Table 1 :

 1 Number of thrusts per thrusters for the extracted on-off control profile at the end of step (A) with two values of the threshold parameter, ςa and ς b .

Table 2 :

 2 Consumptions at the end of steps (B) and (C) with respect to the chosen threshold parameter ς at the end of step (A).

	Threshold	T	P	Equivalence Consumption	SK	Consumption	SK
	ς 1 ς 2 ς 3 ς 4			scheme	for step (B)	performed	for step (C)	performed
		T
	0.89 0.89 0.93 0.93	13	CBE	3.51 m/s	yes	0.31 m/s	yes
		T
	0.89 0.88 0.93 0.93	14	CBE	4.41 m/s	yes	0.33 m/s	yes
		T
	0.9 0.9 0.95 0.95	10	EBE	1.65 m/s	yes	0.24 m/s	yes
		T
	0.95 0.9 0.95 0.95	7	EBE	2.05 m/s	yes	0.17 m/s	yes
		T
	0.8 0.8 0.8 0.8	22	CBE	21.45 m/s	no	0.52 m/s	yes
		T
	0.8 0.8 0.8 0.8	22	EBE	22.52 m/s	no	0.52 m/s	yes


Fuqi Song
email: fuqi.song@ims-bordeaux.fr

Gregory Zacharewicz
email: gregory.zacharewicz@ims-bordeaux.fr

David Chen
email: david.chen@ims-bordeaux.fr

AN ARCHITECTURE FOR INTEROPERABILITY OF ENTERPRISE INFORMATION SYSTEMS BASED ON SOA AND SEMANTIC WEB TECHNOLOGIES

Keywords: Enterprise Information Systems Interoperability, Semantic interoperability, Semantic Web, SOA, Ontology

Enterprises interoperability becomes a necessary and critical task for enterprises as the information and complexity of systems increase exponentially. To adapt existing information systems to new requirements and build bridges between them, a conceptual vision and understanding of interoperability can facilitate the technical development to reduce the conflicts and gaps between heterogeneous systems. Now the focus shifts from the syntactic concern to semantic issues. Semantic heterogeneity currently becomes the next barrier and challenge to face in enterprise interoperability, it promises to play a major role to enable interoperability of enterprise information systems. The approach can be transposed at enterprise decision level, where heterogeneous business processes also need to be interoperable. SOA helps to explore a loosely coupled architecture with many advantages, especially for enhancing interoperability among enterprises. To obtain business process interoperability, service orchestration and choreography could coordinate and compose the services and processes effectively. This paper, first gives a brief introduction of main issues in this domain, and then recalls a framework for enterprise interoperability. Based on these methodologies and the interoperability framework, a new architecture for enterprise interoperability of information systems using SOA ideology and semantic web technologies is proposed, whose final objective is to enable and enhance enterprise interoperability and then create business values.

INTRODUCTION

The development of information systems is a key issue which influences the progress of enterprises nowadays. Since the information systems have been applied for decades in enterprises, numerous information systems had been developped in various domains and industries, so they are coexisting at the same time independently. However isolated information systems cannot adapt the rapid business changes, such as the globalization of production, economic crisis, marketing changes, etc. Thereby new demands that came up requires multiple information systems work together to undertake the tasks that one single system could not complete. This is the problem of interoperability, whose goal is to fulfill the gap between new business requirements and existing information systems in enterprises.

The main basic barrier, which prevents to achieve interoperability, is the heterogeneity, both structural and semantic. Data integration is a solution to solve the problem of information system interoperability before. It mainly solves the structural heterogeneity, because these systems were relatively simple. However, as the amount and complexity of information systems increase, the structural heterogeneity becomes a less important issue, instead, the semantic heterogeneity becomes new challenge and focus.

Besides the semantics issues, the issue of business processs among different information systems is crutial for enterprises. Because the ultimate goal of enabling information system interoperability is to make business interoperable. The objective of business process interoperability is to coordinate and compose the services into process to fulfill new and more complex business requirements.

Thus semantic and business process issues are two significant points for constructing a flexible and effective architecture for information systems interoperability. This paper aims to propose an operational architecture to facilitate interoperability of enterprise information systems from IT and enterprise modeling perspectives. On on hand, with applying SOA ideology and methodology it provides a loosely coupled architecture and enhance the business process integration. On the other hand, the use of semantic web technogloies allows solving semantic heterogeneity.

Framework for Enterprise Interoperability (FEI)

Enterprise interoperability is the ability for two systems to understand one another and to use functionality of one another. In enterprises, solutions need to be implemented from multiple aspects to facilitate interoperability. In [START_REF] Chappell | Framework for enterprise interoperability[END_REF][START_REF] Chen | Architectures for enterprise integration and interoperability: Past, present and future[END_REF] a three dimensional framework (see figure 1) has been proposed to represent the ability of interactions between enterprises systems. First the framework identifies the barriers which are obstacles for achieving interoperability among enterprises: conceptual, technological and organizational. Then four concern points are defined: data, service, process and business, these points should be taken into account when architectures and solutions are built. Generic approaches used for enterprise interoperability are categorized into: integrated, unified and federated.

Figure 1: Framework for enterprise interoperability

Service-Oriented Architecture (SOA)

Service-oriented architecture (SOA) is a flexible set of design principles used during the phases of systems development and integration in computing.

A system based on a SOA will package functionality as a suite of interoperable services that can be used within multiple separate systems from several business domains. Traditional methods are highly coupled, this doesn't facilitate the enterprise interoperability, and is not flexible for scalability, extendibility and maintenance. Service-orientation requires loose coupling of services with operating systems, and other technologies that underlie applications. SOA separates functions into distinct units, or services, which developers make accessible over a network in order to allow users to combine and reuse them in the production of applications [START_REF] Bell | Service-Oriented Modeling: Service Analysis, Design, and Architecture[END_REF].

A service is usually defined as the realization of business functionality via software that anyone can use, anywhere, to compose new applications by using their services in the context of new or modified processes [START_REF] Izza | Exploiting semantic web services in achieving flexible application integration in the microelectronics field[END_REF]. The goal of the web services effort is to achieve interoperability between applications by using Web standards. Web services use a loosely coupled integration model to allow flexible interaction of heterogeneous systems in a variety of domains including business-to-consumer, business-to-business and enterprise application.

Service orchestration and choreography are two terms used to describe the composition and collaboration of services. A service is a unit of functionality, which could be regarded as atomic to provide users certain functions. Most of time, these services are independent with each other, this may cause the waste and redundancy, because one service may be reused for different business processes.

Service orchestration refers to an executable business process that may interact with both internal and external web services, whereas service choreograph means more collaborative in nature, where each party involved in the process describes the part they play in the interaction. Orchestration differs from choreography in that it describes a process flow between services, controlled by a single part. More collaborative in nature (see figure 2), choreography tracks the sequence of messages involving multiple parties, where no one party truly "owns" the conversation [START_REF] Peltz | Web service orchestration and choreography -a look at WSCI and BPEL4WS[END_REF]. The following of the paper is organized as; section 2 analyzes the semantic issues in information systems interoperability. Section 3 describes the proposed architecture and elaborates each layer. In section 4, we discuss the proposed architecture compared with some existing approaches. Last section presents some future work that will be done to this architecture. The methods for data interoperability have a trend shifting from structural to semantic [START_REF] Ziegler | Three decades of data integration -all problems solved?[END_REF].

Because the volume of information increase exponentially and more systems get involved, to find useful information and retrieve them from the huge amount of data becomes necessary. The traditional structure-oriented approaches have limitations for retrieving "meaningful" information, because of insufficient emphasis on semantic aspects. [START_REF] Buccella | Ontology-based data integration methods: A framework for comparison[END_REF] classify heterogeneity into four categories: (1) structural heterogeneity, involving different data models;

Semantic Heterogeneity

(2) syntactical heterogeneity presents different languages and data representations;

(3) systemic heterogeneity, involving hardware and operating systems; and (4) semantics heterogeneity involves different concepts and their interpretations.

In general, the semantic heterogeneity deals with three types of concepts: the semantically equivalent concepts, the semantically unrelated concepts, and the semantically related concepts. In the first case, different models use different terms to refer the same concept, for example, one system may use "teacher" whereas the other uses "professor". In the second case, the same term may be used by different systems to denote completely different concepts, such as "dear" may refer to "expensive" or address a person; and in the last case, different classifications may be performed, for example one system classifies "student" as "male" and "female" and other system as "bachelor" and "master".

Ontology

Ontology is regarded as the main approach to solve the semantic problems. A classical definition considers ontology as a specification of a conceptualization [START_REF] Gruber | A Translation Approach to Portable Ontology Specifications[END_REF]. Ontology can represent the meaning of vocabularies themselves and their relations. Ontology is frequently used in data integration for solving the content explication.

In [START_REF] Wache | Ontology-based integration of information -a survey of existing approaches[END_REF], many approaches for information integration based on ontology are surveyed, and three possible modes are concluded: single ontology, multiple ontologies and hybrid approaches. Each mode has its constraints and benefits, single ontology is easy to implement but is not adaptive when adding or removing a data source, whereas multiple ontology implement costly but has more flexibility with data source changes. Hybrid mode facilitates comparison among different sources due to the shared vocabulary but hard to maintain.

Semantic Web

The semantic web provides a common framework that allows data to be shared and reused across application, enterprise, and community boundaries.

The ultimate goal of the Web of data is to enable computers to do more useful work and to develop systems that can support trusted interactions over the network (W3C, 2010). Semantic web is a group of methods and technologies to enable the construction of a connected and reusable data network. (Bratt, 2007) presented a stack to illustrate the main semantic web technologies.

A SOA AND SEMANTIC WEB THCHNOLOGIES BASED APPROACH

The objective of our approach is to propose a loosely coupled architecture for information systems interoperability. To obtain high interoperability, we explore the techniques and methods that could facilitate and benefit it, thereby SOA and semantic web technologies are applied to design this architecture. The idea is not developing a new dependant and separate information system, but establishing a transparent layer among enterprises information systems so that they can share and exchange information without semantic barriers.

Adaption to the Framework

Based on the enterprise interoperability framework and the approaches used in [START_REF] Zacharewicz | A Multi Agent/HLA Platform for Enterprises Interoperability: Short-Lived Ontology Based[END_REF], we focus in the framework to a specific scope: information system interoperability. We consider all the barriers and three concern points: data, service and process. A federated approach is applied due to the setup of a semantic layer. The mapping is shown in figure 3 to situate the proposed architecture in the Framework.

Overview

The architecture can be applied in one enterprise or among several enterprises (ex: partner relation, collaborations) to enable business corporation among them via their own information systems, and keep their own applications independent. The architecture includes five layers as illustrated in figure 4. Each layer is independent to the others, which communicate through defined interfaces and protocol. The lower levels are more concerned with data, in the contrast, the higher level focuses more on business and users.

Figure 4: Overview of the architecture

In short, the data interoperability part is carried out in the two lower layers, and the service and business interoperability part is done in the two upper levels, the highest level deals with user interaction. The architecture covers three interoperability concerns: data, service and processes. A semantic layer upon data sources is designed to emphasize the semantic heterogeneity of information, not only the information semantics, the semantic web services are also built in the service layer to tackle conceptual interoperability barriers. In addition, the SOA-based methods are used for solving the technological barriers and providing a federated architecture. The business coordination and collaboration take place at process layer, and get end users are involved via the user interface layer, to remove the organizational barriers.

First, a semantic layer is built upon data sources, which is equivalent to the information in the data sources, through a rule-based processor to map between the semantic layer and data sources. It does not aim to transform the data from the data sources to the ontology server, but a mapping connection built between them. When retrieving or updating the data, the real action needs to be effected on the data sources through the semantic layer. And a monitor component is set up to keep semantic adaptive and evolutionary when the information and semantic change.

Then the upper level is service layer, in which the services are regarded as atomic. The services layer is core part of the architecture, because it provides functionalities and connects information and business processes. Namely, it retrieves information from lower level to implement the services, and then provides the atomic services to compose business process. Compared with the single services in this layer, the business process layer is linked to the real business process performed in enterprises. The processes are composed by the atomic services in services layer. In order to integrate the services into process, services orchestration and choreography are applied for coordinating or composing services. In particular, the description languages WS-BPEL (Web Services Business Process Execution Language) and WS-CDL (Web Services Choreography Description Language) are used for service orchestration and choreography respectively.

The highest level is the user interface level, which the user interacts with directly, in this layer the systems are the ones existed in enterprises or the future system to be developed. It connects to the enterprise service bus (ESB) with business process layer and service layer, in order to call the services to accomplish the business processes.

Elaborations of Each Layer

In this section, we define each layer of the architecture; the detailed illustration is shown in figure 5.

Data Sources Layer

This layer represents the data sources located in different enterprises and domains, They can be accessed directly or through the interfaces they provide. These sources have three characteristics: heterogeneous, distributed and autonomous. The details of each characteristic are as follows.

Heterogeneous: They have different structure and format, the most common storage method is rational database, which is a well-structured format. Another one is file. However, the formats of files vary a lot depending on its property. For example HTML/XML/XHTML is well-structured sources, whereas excel sheet, word document are semistructured sources. As well as non-structured sources, such as plain text document (.txt) and some other documents with special designed format.

Distributed: The sources locate on different domains either physically or logically. Such as in one enterprise, the database could be deployed to several cities, because there is one branch in each city. Even two databases are deployed in one same server, they belong to different database system, for example, Oracle and DB2 respectively, this is an instance of logically distributed. Autonomous: Each source is logically autonomous, each source can be regarded as one system, and they are dependent to each other. They are connected to the systems via interfaces or API (Application Programming Interface).

Semantic Layer

Upon the structural data sources, a semantic layer is established to represent the semantics of the information. Not all information will be taken into account when building the semantic layer, but only the information which is useful and valuable to the integrated business. Thus a rule-based processor is constructed to map and retrieve data from the lower data sources layer. With defined rules, the processor could filter and control the data so as to build a concise and effective semantic layer. The processor can also update the data sources when a change request is sent from the upper level.

Ontology is used for representing the semantics of information. An ontology server is deployed as the storage of the ontology. As we introduced in the introduction in the recall section, the ontology server is applied with the global schema approach due to its advantages. A meta-data storage is designed to abstract the global ontology vocabularies. Furthermore it could be used to extend the global ontology schema to multiple and hybrid schema when single mode has too much limitations and constraints for a specific future system.

There are many languages for representing ontology, each one has its highlight and constraints, the selection will depend on the specific focus and requirements. (Pulido et al, 2006) and [START_REF] Maniraj | Ontology Languages -A Review[END_REF] give a review of the ontology languages, (Corcho et al, 2002) presents the different ontology languages with detailed comparisons and introduces the tools for building ontologies.

User interface layer

Monitor

Another important issue when building semantic is to keep it adaptive and evolutionary. Because the information in data sources changes frequently, such as adding a new database or removing some files. The ontology needs to keep consistent with the real information layer. Even the data itself doesn't change, the meaning of information could change date by date. Because we are dealing with semantics, when the meaning of a word changes, the old ones become meaningless and incorrect, so keeping evolutionary is necessary and important. To do these tasks, a component called monitor is set up to keep the ontology update-to-date and correct.

Services Layer

Service layer lies above the information and semantics, these services are the existing and new developed ones located in different information systems and enterprises. Traditionally, these services are independent and not related; however, they could be reconstructed to compose into a process that fulfills certain tasks. The upper level will perform the service integration. In this layer, the services are gathered to establish a semantic web services layer, compared with the conventional web services, semantic web services has the ability to understand the services and the relations between them.

A component called query engine is designed for querying and answering between the semantic layer and the services layer. The query engine is ontology oriented, such as, SPARQL (W3C recommendation, 2008), so that it can query the ontology in the semantic layer. An optimizing mechanism is set up in this query engine, for optimizing the query request sent from application level and the answers.

Business Process Layer

Business process layer represents the business activities in enterprises, these business processes happen internally or among enterprises. These processes are composed with the atomic services in services layer. Service orchestration and choreography are applied to coordinate and compose services into business process, for example WS-BPEL and WS-CDL are specific description language for business process integration. WS-CDL is an XML-based language that describes peer-to-peer collaborations of parties by defining, from a global viewpoint, their common and complementary observable behavior; where ordered message exchanges result in accomplishing a common business goal (W3C, 2004). WS-BPEL is an OASIS standard executable language for specifying actions within business processes with web services, WS-BPEL export and import information by using web service interfaces exclusively (OASIS, 2007).

A business integrator is designed to perform the process integration task, and in charge of connecting to the ESB (Enterprise Service Bus) [START_REF] Chappell | Framework for enterprise interoperability[END_REF]. The dotted line between the process and the services in the figure 7 is logically connected. It means a certain process is composed with several services, but physically, the process and services are connected via ESB.

User Interface Layer

This layer represents the user level, namely, the interface the users interact with, such as the human resources system, orders management system, inventory management orders and etc. Via the defined interface, the applications connect to the ESB and explore the process and service to complete user requirements. New modules and functionalities could be added and developed within the current systems. So this could keep the dependability of the current systems, they don't need to be changed to adapt the new requirements, because the cost for updating the legacy or developed systems is very high, we try to minimize the unnecessary changes and maximize the value-added features.

DISCUSSION

In this paper, we addressed the importance of semantic and business process for enterprise information system interoperability. Based on the ideology of SOA, we proposed a new architecture to integrate enterprises information systems, which uses semantic web technologies to construct a semantic layer upon data sources, and service orchestration and choreography to integrate business process. The top level is the information systems, with which the end users interact directly. The applications connect to the ESB through the defined interface, and then invoke the processes and services to fulfill business requirements. [START_REF] Ziegler | Three decades of data integration -all problems solved?[END_REF] sums up the usual methods for data integration by six levels of integration. They are described from the architectural perspective of information systems. In particular, at the user level only "Manual" integration, which has a very low integration degree, is available. The most often used approach is "Uniform Data Access" to solve the data integration problems, which provide user a unified user interface of physically distributed data, such as portals and mediated query systems. "Common data storage" uses a method to transfer the data sources to a new data storage, warehouse and operational data stores are such kind of instances.

Compared with these integration approaches, interoperability has more the meaning of coexistence, autonomy and federate environment, whereas integration refers to the concepts of coordination, coherence and standardization [START_REF] Chen | Architectures for enterprise integration and interoperability: Past, present and future[END_REF]. The interoperable architecture will provide added values for enterprise information systems interoperability. First, due to the semantic technologies used in the architecture, the "understanding" problems can be facilitated among different information systems. Secondly, with SOA, the architecture is loosely coupled, where each layer is only connected through interface and protocols. The data sources, services and business processes could be added or removed without affecting the other parts of the architecture, gaining flexibility for the enterprise. Thirdly, it provides high scalability; this architecture could be applied internally in one enterprise for its sub branches or different departments, as well as among different enterprises and partners. Involved applications and services can be reused ones; also new developed ones can be added easily. We believe this approach is promising regarding web 2.0 orientation of all enterprise applications.

FUTURE WORK

The future work will include, first, more details to describe the running mechanism of the architecture, such as, the interfaces between each layer. Second, we will elaborate and refine the semantic layer to define how the ontology presents the semantics of information; third, the aspects of enterprises business level (e.g. decision level) will be taken into account to make the architecture adapted the complex situation in enterprises. The architecture will be applied to a project in the following research, to figure out the improvements and get feedbacks from industry.

Figure 2 :

 2 Figure 2: Web services orchestration and choreography

Figure 3 :

 3 Figure 3: Position of the proposed architecture in the framework

Figure 5 :

 5 Figure 5: Five layers of the architecture