Predictive control algorithm for spacecraft rendezvous hovering phases

The paper proposes a predictive control law for orbital rendezvous hovering phases. The proposed algorithm provides an impulsive control that steer the vehicle to a set of the periodic relative orbits enclosed in a hovering zone. The control law is computed by determining a point lying in the intersection between a semi-algebraic set and a hyperplane using an alternating projections algorithm. The present work proposes several improvement compare to previous published works based on a two-impulse strategy: the control consists of one impulse instead of two, the periodic and admissible orbit that the chaser is steered to is not defined a priori, and the budget and saturation conditions are taken into account. The efficiency of the proposed predictive control algorithm is bench tested on a nonlinear simulator.


I. INTRODUCTION
In the future missions of on-orbit inspection and on-orbit servicing, the ability to maintain a relative periodic motion is a crucial issue for any rendezvous spacecraft. The purpose of this research is to develop a predictive control law that enables the control of periodic relative orbits that are included in a given linear region for sake of safety.
Constrained relative trajectories have been studied under the term of "hovering". Hovering has been defined in [1] as follow To develop a control strategy to place a deputy satellite inside a specific lobe defined in the chief body-fixed frame and keep it there in a fuel-optimal manner.
In most of works on the subject, as [1]- [7], the hovering strategy is based on two different strategies. First, the "teardrop" strategy consists in computing an impulsive thrust such that the initial and final relative positions on the bounds of the tolerance subset are equal for a given time of flight [1], [8]. Second, the "pogo" strategy, where the impulsive control is computed each time the vehicle hit the bounds of the lobe to maintain it as longtime as possible while minimizing the fuel consumption [1], [4], [7].
However, some contributions have recently taken advantage of the naturally periodic relative orbits to design hovering control laws [9]- [11]. In [9], a methodology to compute the coordinates of periodic relative orbits included in a polytopic region is exposed. This method has been extended to the non periodic relative trajectories in [11].
A control law proposed in [10] consists of two impulsive maneuvers that stabilize the spacecraft with respect to a given constrained periodic relative orbit previously computed. Even though this control law has been proven stable for any eccentric target orbit and requires very few computational effort, its optimality in terms of fuel consumption is not guaranteed.
The present article proposes a different strategy to hover inside a given polytopic subset of the relative state space. The purpose is to steer the vehicle to the set of the periodic relative orbits enclosed in a hovering zone.
To achieve such a goal, a predictive methodology is applied to provide the impulsive control. At each call, the controller computes, if it exists, the coordinates of a relative periodic orbit included in the hovering zone that go through the current position. The impulsive control is then deduced from these coordinates. Moreover, during the control computation, spacecraft abilities and propellant resource are also taken into account.
Although the present work extends the results from the conference article [12], it differs in several points. In the present paper, the chaser is steered to a periodic and admissible orbit not defined a priori, contrary to what is done in [12]. Consequently, the control consists of one impulse instead of two for the techniques from [12]. Finally, the proposed methodology permits to consider in the problem saturation and propellant budget constraints.
As it is exposed thereafter, the set of the orbits coordinates admissible with respect to the above-mentioned constraints is described as the intersection between a semi-algebraic set and a hyperplane. Obtaining the control mainly consists in computing one point of this intersection set. This is achieved using an alternating projections algorithm. Note that this algorithm by means of its properties make the use of a complete optimization procedure unnecessary as explained in this paper.
The paper is structured as follows. Section II gives a brief description of the relative motion and recall the trajectories parametrization used in the sequel. Section III describes, in terms of the relative trajectories coordinates, the constraints on the relative orbits and the controls. Section IV presents the alternating projection algorithm that permits to compute the coordinates of an admissible relative trajectory and the control. suming Keplerian orbits, the relative dynamics can be given by the linearized and simplified Tschauner-Hempel equations [13] (time is replaced by true anomaly as independent variable): where e is the leader orbit eccentricity. Defining the vector T , a state-space representation for the non-autonomous system (1) is: Note that the state expressed in time and expressed in true anomaly are linked by means of the following similar transformation: Using the transition matrix, Φ(ν, ν 0 ), proposed by Yamanka and Ankersen in [14] to propagate a initial state X(ν 0 ) and considering an unforced motion: Expanding eq. (4), the following relative trajectories description can be deduced: where c ν = cos ν, s ν = sin ν, J ν ν0 is given by: and T , the parameter vector, is expressed as a linear function ofX(ν), Please refer to [11] for a more detailed equation (7). 1 Earth-Centered Inertial 2 Local Vertical / Local Horizontal

III. CONSTRAINED PERIODIC ORBITS DESCRIPTION
In this section, the admissible trajectories are described. The admissible orbits must have two main properties: the periodicity and the inclusion in a polytopic set. These two properties are described thereafter as constraints on the parameter vector D.

A. Periodicity
The relative trajectories between spacecraft are not generally periodic. The periodicity constraint guarantees a relative motion that does not require any correction in absence of disturbances. ImposingX(ν + 2π) =X(ν), ∀ν ≥ ν 0 one can deduce that this constraint can be expressed as: Analyzing the equations in (5), condition (8) becomes evident, since the only non-periodic and divergent term J ν ν0 is always multiplied by d 0 .

B. Polytope constraints using non-negative polynomials
During the rendezvous hovering phases, the follower is required to remain in the interior of a certain limited region of the space. This region is going to be modeled as a polytope (without any generality loss, a rectangular parallelepiped will be used): Defining X pos (t) = [x(t) y(t) z(t)] T , this constraint can be rewritten using the matrix H = [I 3×3 − I 3×3 ] T and the vector K = [x max y max z max x min y min z min ] T : Replacing the time t by the true anomaly of the leader ν as independent variable, we obtain: Assuming that the periodicity constraint is satisfied and introducing the following variable change: the trigonometric terms in (5) are eliminated and we obtain the following rational expressions: the coefficients p xi , p yi and p zi depending linearly on D: Using the variable change (12) and the rational expressions (13) in (11), we can obtain the following polynomial inequalities: that can be reformulated as: The non-negativity conditions in (16) result in an infinite number of constraints to be respected. Deaconu presents in [9] a way to transform this infinite number of constraints in a finite number of convex constraints using the following theorem presented by Nesterov in [15]: Theorem 1: Let P (w) be a polynomial of degree 2n on R and p ∈ R 2n+1 the vector containing its coefficients. Then: where S n+1 + is the cone of the symmetric semi-definite positive matrices of size n + 1 and H n,l are the Henkel matrices: This is equivalent to affirm that the cone of non-negative univariate polynomials is the image of the cone of positive semi-definite matrices by the linear operator Λ * (·).
Let γ i be the vectors containing the coefficients of the polynomials Γ i (w). Using (17), we can transform conditions in (16) into: In the end of this section, constraints on the control are given: the thrusters saturation and the propellant budget. Considering saturation constraints is one of the main advantages of the predictive control law. Moreover, the propellant is a rare resource for spacecraft systems since they can not be refueled. Thus, propellant budget constraints is also proposed to take into account the scarceness of the propellant.

C. Thrusters saturation
Let ∆V = [∆V x ∆V y ∆V z ] T be the vector representing the impulsive velocity correction executed by the follower. Considering that the spacecraft have a pair of equivalent propellers symmetrically and oppositely disposed on each axis and assuming that the saturation limit for each propeller is ∆V max > 0, the saturation constraint can be written as: These inequalities can be transformed in a matrix positive semi-definite constraint: Mi 0, s = x, y, z (21)

D. Propellant budget
An maximum budget of propellant σ is available for each impulsion executed. This constraint is expressed as: To check condition (22), we must verify if: IV. COMPUTING THE CONTROL LAW Accordingly to the model predictive control scheme, the control is computed periodically in order to bring the system on a periodic trajectory inside a given polytope respect to actuator saturation and budget constraints. Using the previous equations, links between the variables of the problem and its restrictions can be stated. The Fig. 2 presents these links, each ellipse representing a constraint.
Constrained variables Transformations H,K Given a set of constraints {H, K, ∆V max , σ} and an initial relative state X 0 , we want to determine an impulse ∆V respecting both saturation and budget constraints such that the state after the impulse, X + = X 0 +B∆Ṽ ,equivalent to the state D + (by means of (7)) whose first component d 0 is zero to copt the periodicity and that generates (by means of (14), (15) and (16)) six non-negative polynomials to ensure the polytopic inclusion of the orbits

A. Closed convex sets intersection problem
The constraints can be summarized as follows: Polytope constraints: Thrusters saturation: We remark that these constraints are convex, since they are mathematically represented by matrices positive semidefiniteness (on the left) and equalities (on the right).
Let Q be the symmetric block diagonal matrix given by: This matrix is sparse by construction. Table I shows that taking symmetry into account we only need 49 variables to completely describe it. We can use Q to reformulate the positive semi-definiteness conditions and the equations linking the variables in (24) as follows: where S 32 is the set of 32 × 32 symmetric matrices and m is the number of equations.
Equations tr(A i Q) = b i describe an affine subset of S 32 and Q 0 represents the symmetric semi-definite matrices cone S 32 + . Therefore, solving the problem is equivalent to find an element Q belonging to the intersection of these sets.

B. Alternating projections algorithm
For two closed and convex sets, we can apply the alternating projections algorithm proposed by Boyd and Dattorro in [16] to determine the existence of an intersection: Let C and D be two closed convex subsets of E, let || · || E be a norm on E and let P C and P D denote the projections on C and D, respectively. The algorithm consists in:

Algorithm 1 Alternating Projections
Choose an arbitrary Q 0 and compute Q C It has been proven that [17]: i.e. a pair of points in C and D that have minimum distance.
Thereafter, let C be the affine set described by the equations from system (26), let D be the symmetric semi-definite matrices cone S 32 + and let ||·|| F be the Frobenius norm. Since these sets are closed and convex, the alternating projections algorithm can be implemented to determine an element belonging to their intersection. In this case, projection P C is given by: where a i are solution of the following equation: where The projection P D is given by: by considering the spectral decomposition is written as

C. Initialization and stopping criteria
As any iterative algorithm, the projection algorithm's behavior depends on two items: the initial guess and the stopping criteria. The initial guess to provide to the alternating projection algorithm is the matrix Q 0 . In fact, the alternating projection algorithm provides to the initial guess the nearest element in the intersection of the subsets (if exists) [18]. Thus, the choice of Q 0 can have a great impact the solution Q * as can be seen of 3. In this paper, two strategies are tested: the initial guess Q 0 is either the null matrix or the solution Q * of the previous call of the projection algorithm. As it is exposed in section VI, different strategies for choosing the initial iterate will lead to different controlled trajectories. The stopping criteria is based on the distance in the Frobenius norm between Q D k and Q C k : The algorithm stops when where ǫ Q > 0 and ǫ δ > 0. The first condition means the iterate Q k is close enough to the intersection and the second condition reveals that the algorithm reaches a minimal and non zero distance between the convex subsets so that the intersection is an empty set.

V. STABILITY BY MEANS OF INVARIANCE
The stability around a periodic trajectory is understood in the sense given in [19], which requires the characterization of the evolution of the trajectories that start arbitrarily close to the periodic one. According to [20, theorem 3.5], the stability of periodic system is characterized by the eigenvalues of its monodromy matrix M = Φ(ν + 2π, ν). In the case of the linearized relative motion, M has at least two eigenvalues µ i such that |µ i | > 1. This leads to the conclusion that the autonomous systemX ′ (ν) =Ã(ν)X(ν) is naturally unstable. Thus, without control, the trajectories, described by (4), that started arbitrary close to any periodic trajectory will diverge from it.
The proposed predictive control law has for goal to steer the system in the vicinity of the set of admissible periodic orbits. Thus, the stability will be defined in the sequel as the invariance of given subsets described thereafter. Now, two sets are given in order to define the stability of the system controlled by means of the proposed predictive algorithm: • S 0 {H, K}, the set of vectors D(ν) ∈ R 6 that describes periodic trajectories satisfying the polytope constraints; • S 1 {H, K, ∆V max , σ}, the set of vectors D(ν) ∈ R 6 for which it is possible to find one impulse ∆V , respecting the saturation constraints and the propellant budget, such that D + (ν), the state right after the impulse, belongs to S 0 {H, K}.
The set S 0 {H, K} naturally invariant because of the periodicity respected by each of its elements.
The set S 1 {H, K, ∆V max , σ} is also a class of equivalence under the developed control law, since it is possible to compute an impulsive correction via the alternating projections algorithm that generates a new state belonging to S 0 {H, K}. Consequently, the set S 0 is invariant and attractive under the developed control law with S 1 as convergence set. If the current state does not belong to S 1 , an ad hoc control law must be used to bring the system in S 1 . For instance, an optimal impulsive plan computed using techniques from [11] can bring the chaser satellite to the vicinity of the hovering box from any relative state.
The Figure 4 represents the sets S 0 and S 1 and the stability and instability zones in R 6 .  This mission has been implemented on a simulator based on Gauss nonlinear equations for the relative motion, taking into account the atmospheric drag and Earth's oblateness effects (see [21] for details).
In order to assess the impact of the initialization techniques described in IV-C on the resultant trajectory and the consumption, a warm start initialization, taking the last iteration solution as initial guess, and a cold start initialization, using a constant Q 0 matrix for all the iterations, were performed.
Simulations with different eccentricities and time-intervals between control calls are conducted to test the algorithm with both initialization procedures. The aim of these simulations is first to verify the designed control law efficiency in respecting the different constraints while rejecting the unmodeled dynamics, and second to highlight the importance of the algorithm initialization procedure. The results obtained for these simulations are presented in Figure 5. We can remark that, on one   hand, the smaller eccentricity, the lower consumption is. On the other hand, the consumption increases with the period between controls. Another conclusion can be drawn from Figure 5: the control algorithm using a cold start technique generates greater consumptions compared to the one using a warm start technique. This can be explained by the fact that cold start leads to computation of parameters of admissible orbit that evolve much more than the one computed from warm start procedure.
On Figure 6, controlled trajectories obtained by means of the nonlinear simulator previously described are exposed. The control algorithm is called every 20 second for this specific simulation. First, the simulated trajectories remains admissible during the 5-orbit period of hovering. However, it can be observed that both initialization procedures lead to different geometries respecting the tolerance box constraints despite the orbital perturbation. The total ∆V consumption over the 5 orbital periods is of 4.4 mm/s and 18 mm/s for the warm start and cold start control respectively. Note that, the ∆V budget and saturation constraints are also satisfied. These facts highlight that the admissible orbit computed using a warm start initialization is closer to current position than the one computed from a cold start. In fact, after every call of the control algorithm, the chaser will evolve close to the computed orbit. With the warm start initialization, the next computed orbit will remains in the vicinity of the previously computed orbit. In the contrary, the cold start procedure leads to the computation of a completely different admissible orbit. The present paper develops a predictive control law for the orbital rendezvous. This control law aims to make the chaser satellite hover in a given tolerance set in order to maintain a way point situation. As the relative motion has a periodic dynamics, the predictive control algorithm take advantage of the existence of periodic orbits for this system to develop a law that maintain the hovering while comsuming propellant essentially to reject perturbation, contrary to pogo or tear drop strategies available in the literature. This control law consists in steering periodically the chaser spacecraft into the subset of the admissible orbits. Moreover, budget and saturation condition are also accounted in the control computation. One of the main features of the proposed control law is that it is based on the alternating projection algorithm. This iterative algorithm provides to any initial guess the nearest point in the intersection subset of convex sets when it is not empty. We studied the influence of the initial guess choice on the controlled trajectories and the consumption through nonlinear simulations. Further studies should also focus on the numerical issues like convergence time as function of the stopping criteria tolerance. Another problem of interest is the characterization of the invariant set S 1 , that could speed up the decision of whether or not to apply the projection algorithm, once it should be possible to determine the algorithm convergence by verifying the belonging of the current state to S 1 .