# Robust Rendezvous Planning Under Maneuver Execution Errors

Christophe Louembet, Denis Arzelier, Georgia Deaconu

► To cite this version:

Christophe Louembet, Denis Arzelier, Georgia Deaconu. Robust Rendezvous Planning Under Maneuver Execution Errors. Journal of Guidance, Control, and Dynamics, 2015, 38 (1), pp.76-93. 10.2514/1.G000391. hal-01078508

HAL Id: hal-01078508

https://hal.science/hal-01078508

Submitted on 29 Oct 2014

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Robust rendezvous planning under maneuver

execution errors

Christophe Louembet<fnref n="1" />, Denis Arzelier<fnref n="2" />, Georgia Deaconu<fnref n="3" />

CNRS; LAAS; 7 avenue du colonel Roche, F-31400 Toulouse; France

The problem of designing rendezvous guidance maneuver plan robust to thrusting errors is addressed in this paper. The aim of this paper is to develop tractable and robust guidance algorithms. Solving the rendezvous guidance problem via a direct approach leads to uncertain optimization problems while accounting for the Guidance, Navigation and Control (GNC) systems uncertainties and errors. A worst-case approach is considered in order to obtain tractable robust counterparts. The robustness certificates derived from these guidance algorithm provide the means to analyze the effects of the considered errors on the rendezvous mission. Several types of missions tested in a linear environment are used to illustrate the methodology.

Nomenclature

$a$ = semi-major axis ;

$e$ = eccentricity ;

$\mu$ = Earth gravitational constant;

$\nu$ = true anomaly ;

$\mathbf{X}$ = relative motion state;

<footnote n="1">¹ Univ de Toulouse; UPS, LAAS; F-31400 Toulouse; France.
Associate Professor, Université de Toulouse; louembet@laas.fr</footnote>

<footnote n="2">² Univ de Toulouse; LAAS; F-31400 Toulouse; France
Directeur de Recherche, LAAS-CNRS, Methods and Algorithms in control</footnote>

<footnote n="3">³ Univ de Toulouse; UPS, LAAS; F-31400 Toulouse; France.
Ph.D. graduate, LAAS-CNRS, Methods and Algorithms in control</footnote>
$\Phi(\nu_i, \nu_j)$ = transition matrix of relative motion from instant $\nu_j$ to $\nu_i$;

$N$ = number of velocity increments;

$\nu_i$ = impulses application times $i = 1, \dots, N$;

$\Delta V_i$ = velocity increment vector at $\nu_i$;

$\Delta V_{max}$ = actuators saturation in meter per second;

$u$ = disturbance variable;

$\nu$ = uncertainty set;

$\mathcal{T}$ = tolerance set;

$\Gamma$ = tolerance vector;

$0_l$ = null matrix of dimension $l$;

$\mathbb{I}_l$ = identity matrix of dimension $l$;

$1_l$ = vector of length $l$ where every entry is 1;

$\| \cdot \|_l$ = the $l$-norm of a vector;

$| \cdot |$ = the absolute value function;

# I. Introduction

Any orbital operation involving at least two vehicles (satellite or International Space Station servicing, debris removal and formation flight) include a rendezvous process which consists in bringing the actuated vehicle (chaser) into the proximity of the passive vehicle (target) via a series of elementary phases (phasing, far range rendezvous and close range rendezvous). Each phase is composed of orbital maneuvers that must be carefully planned to reach their goal. This paper is particularly dedicated to the far range rendezvous operations, which begins when the separation between the chaser and the target is sufficiently small to allow the relative navigation between the two spacecraft.

From a mathematical point of view, the rendezvous guidance problem consists in an impulsive fuel-optimal control problem when the gas ejection thrust is assumed to be impulsive. This particular
optimal control problem can be solved using the so-called primer vector theory that originates from
the seminal work of Lawden [1]. This theoretical corpus relies on optimality conditions derived
from Pontryagin's Maximum principle. This indirect approach has been the base for numerous
studies addressing the orbital rendezvous problem [1-6]. However, the algorithms from previously
cited works are, in general, tailored for off-line computations and are not intended to be part of the
spaceborne guidance software. In this context, numerical solutions based on direct methods [7] are
particularly appealing. In fact, direct methods consists in parametrizing the original rendezvous
guidance problem and converting it into a programming problem for which efficient solvers are
available. These approaches are known to be particularly effective in presence of state and input
constraints even if they generally do not provide the optimal solutions but only an approximation.
The application of direct methods to the rendezvous guidance problem was initiated by Waespy [8]
and later developed by Robertson et al. in the late nineties [9]. In these works, the constraints
are usually discretized in order to obtain a linear program. However, paper [10] recently proved ,
by converting the rendezvous problem in an semi-definite programming problem, that it is possible
to satisfy the continuously in time without discretisation. Numerous paper addressing different
rendezvous control problems, ranging from coordination between several spacecraft [11, 12], to
collision avoidance [13-15] often in conjunction with the Model Predictive Control (MPC) [11, 12, 16]
have been published since then. If impulsive thrust and linearized dynamics are assumed in the
previous cited papers[8-16], the direct methodology also permits to assume different type of thrust
model or non linear dynamics for rendezvous problem. In the case of continuous thrust, the solution
of the constrained and non linear optimal control problem has been obtained in [17] by performing a
sequence of convex and linear constrained optimal control problems solved by either a transcription
or a shooting method. The advantage of such a method is that the model can be augmented easily
with orbital disturbances without increasing its complexity. In [18], the thrust is modelled by pulse
and the optimal control is obtained by modulating the amplitude of the pulse addressing the fact
that the thrust time can be of few minutes width. In [19], the pulse width modulation is considered
and the optimal control problem is transformed in a 2-step linear program.

Among the challenges that arise from the needs of more spacecraft autonomy, the robustness
of the control method to the uncertainties inherent to GNC systems is a key issue [20]. These
uncertainties may originate from errors in the navigation system and/or the propulsion devices. For
example, in the far-range rendezvous phase, the chaser relative position and velocity provided by
the navigation system may be corrupted by measurement noise. In addition, the thrusts delivered
by the propulsion devices may be inaccurately executed. These errors usually lead, in real imple-
mentation, to the fact that open-loop maneuvers might steer the system far from the rendezvous
objectives.

In the case of infinite control horizon, closed-loop schemes may help to mitigate GNC system's errors and ensure system stability. In the literature, numerous works propose solutions to reject perturbations and disturbances of various types. The case of impulsive rendezvous has been addressed using the discrete linear quadratic regulator (LQR) techniques [21] (see the review paper [22] and references therein). Adaptive control techniques compose a family of control laws that enables the rejection of propulsion errors in the context of formation flight by exploiting Lyapunov stability criteria as in [23]. The receding horizon MPC techniques (see, for instance, [24, 26, 27]) have also been successful to tackle the problem of thrust errors mitigation. However, when it comes to the finite horizon rendezvous problem, it is much more complicated to get equivalent closed-loop certificates of stability and performances (in terms of precision) using these previously mentioned techniques.

In [25], a tube-MPC scheme has been tailored to compute feedback control laws robust to navigation uncertainties while guaranteeing a given precision at the end of the rendezvous mission. Note that the errors addressed in [25] are different in nature from the propulsion errors tackled in this paper.

Indeed, navigation uncertainties are usually treated, in the literature, as additive disturbances while thrust errors are addressed in the sequel as multiplicative disturbances (see [23]). At the best of our knowledge, the generalization of tube-MPC schemes to cope with multiplicative uncertainties like the ones affecting propulsion system is a very difficult challenge.

In this paper, the rendezvous planning problem is solved using linear programming problems
arising from the application of a direct methodology. In order to compute a robust maneuver plan
that provides a given level of rendezvous precision, the idea is to define an a priori model of the
uncertainties and errors affecting some data of the program. It is therefore necessary to resort to
specific techniques for solving robust optimization problems. The usual way to address an uncertain
program is to introduce in the program the probabilistic description of the uncertain data, consid-
ering them as random values. Then, the techniques from the stochastic programming theory can be
applied [29]. This has been done in numerous papers in the case of continuous thrust rendezvous
(see [30] and references therein), but few are those that treat the impulsive case. For instance, Luo
et al. modelled the impulse execution error as a zero-mean white noise process, leading to a multi-
objective program that is solved using genetic algorithms with no polynomial time convergence [31].
Other than the numerical complexity, this type of probabilistic approach can only provide proba-
bilistic guarantees for the constraints satisfaction.

The other way to deal with uncertain program consists in a deterministic worst-case approach. Sev-
eral works in the literature have addressed the robustness of the maneuver plan in the context of
navigation errors. How and Villerson described in [28] a "multiple-model" robust approach, that
consists in designing the input sequence to simultaneously satisfy the constraints for several initial
conditions. It is stated in [11] that testing between 1 and 10 initial conditions can lead to satisfying
results. However, the authors do not provide any certificate of robustness for this approach, they
only assert that such a choice is acceptable in practice thanks to properties of the linear propagation
of a convex hull. A worst case approach for dealing with navigation uncertainties which provides
a guaranteed certificate is described in [15]. Mueller and Larsson modelled the errors affecting the
initial state as belonging to an ellipsoid and then, presented a algorithm based on linear program-
ming that ensure the collision avoidance on a discrete time horizon.

If stochastic programming take advantage of the fact that probability distributions governing the be-
havior of GNC systems are usually available, the resolution technique have a high level of numerical
complexity. In addition, stochastic programming approach will provide probabilistic guarantees for
the constraints satisfaction while most of constraints cannot be softened in the rendezvous guidance
problem. For these reasons, a deterministic worst-case approach guaranteeing a priori robustness
and with good tractability is preferred in this paper.

The contribution of the present paper is to provide robust rendezvous planning algorithms
that guarantee the final rendezvous conditions in the presence of GNC uncertainties and errors.
Technically, the aim of the proposed guidance algorithm is to prevent the spread of these GNC
systems errors by working out the direct method algorithms through robust convex optimization
[33]. The main advantage is to provide the mission designer with constructive numerical algorithms
built on direct optimization schemes that also give a certificate of robust optimality as a by-product.
This approach makes it possible to analyze the worst possible effects of the GNC errors on a specific
mission. Note that uncertainties must be described through their bounds to match the requirement
of the robust optimization contrary to the stochastic programming.

The proposed guidance algorithms are built based on a three step methodology. First, a deterministic modelling of uncertainties model is proposed. Then, the effects of these errors on the optimal rendezvous guidance algorithm are exposed. Finally, a rendezvous guidance algorithm robust to above uncertainties that computes the best guaranteed precision performance is proposed. Please note that the scope of this paper is limited to propulsion devices errors (i.e. firing time and impulse execution errors). In fact, the proposed approach is not relevant to preclude the effects of the measurement noise. This issue is tackled separately by the authors using a different methods in the article [25].

This paper is organized as follows. Sections II provides a description of the rendezvous guidance problem as it is usually presented in the literature and proposes a polytopic relaxation of the rendezvous condition. Section A presents the general guidance problem for relaxed final conditions when uncertainties are considered along with the paradigms of the robust convex optimization. Propulsion devices errors are treated in sections III, IV and robust guidance programs are presented. Finally, section V presents the results obtained for several mission examples in order to validate the methodology, together with nonlinear simulations results.

## II. Direct resolution methods for the rendezvous guidance problem

After presenting the relative motion model and the guidance problem, the rendezvous guidance
problem is formulated as a linear program using a direct approach. If the rendezvous condition is
first defined as an exact condition, it is then shown that rendezvous condition must be relaxed to
handle errors and uncertainties in the rendezvous guidance algorithm.
A. the rendezvous guidance problem

This paper concentrates on the homing phase of the spacecraft rendezvous mission, which begins when the separation between spacecraft is sufficiently small to allow the guidance of the chaser spacecraft using relative navigation. In this framework, the spacecraft relative motion can be characterized expressed in the Local-Vertical-Local-Horizontal (LVLH) frame centered at the target spacecraft position (see figure 1).

<figure><img src="image_2.png" /><figcaption>Fig. 1: Illustration of the LVLH frame</figcaption></figure>

It is assumed that the chaser is moved using several ergol thrusters rigidly mounted on the vehicle frame. Under this assumption, the control can be modelled as impulsive signals that instantaneously affect the velocity component of the chaser relative state, $X \in \mathbb{R}^6$, in the LVLH frame. Under the impulsive control assumption, the linearized spacecraft relative dynamics of the relative state (position and velocity) $X$ expressed with the true anomaly of the leader spacecraft $\nu$ as the independent variable are given by

$$ \frac{d\mathbf{X}(\nu)}{d\nu} = A\mathbf{X}(\nu) + \mathbf{B}\sum_{i}\Delta\mathbf{V}(\nu_i) \quad (1) $$

where $A$ and $B$ are the matrices corresponding to the Tshauner-Hempel linearized and simplified dynamic equations [35]. The control vectors $\Delta\mathbf{V}(\nu_i)$ are defined in the leader LVLH frame such that $\Delta\mathbf{V}(\nu_i) = \Delta\mathbf{V}_i\delta(\nu - \nu_i)$ where the vector $\Delta\mathbf{V}_i$ belongs to $\mathbb{R}^3$ and $\delta(\nu - \nu_i)$ denotes the Dirac impulse at time $\nu_i$.
The rendezvous guidance problem requires the computation of a plan of N impulsive maneuvers {$\Delta \mathbf{V}_i$}$_{i=1,...,N}$ such that the finite-time trajectory $\nu \mapsto \mathbf{X}(\nu)$ satisfies the exact initial condition $\mathbf{X}(\nu_1) = \mathbf{X}_1$, the final condition $\mathbf{X}(\nu_f) = \mathbf{X}_f$ (see Fig. 2) and the actuators constraints, while minimizing a given cost (generally the ergol consumption).

<figure><img src="image_3.png" /><figcaption>Fig. 2: Deterministic rendezvous illustration</figcaption></figure>

The spacecraft rendezvous guidance problem can be stated as follows:

$$
\begin{array}{ll}
\min_{\Delta \mathbf{V}, \nu_i} & J(\cdot) \\
\text{under} &
\left\{
\begin{array}{l}
\displaystyle \frac{d\mathbf{X}(\nu)}{d\nu} = A\mathbf{X}(\nu) + \mathbf{B} \sum_{i=1}^{N} \Delta \mathbf{V}_i \delta(\nu - \nu_i) \\
\mathbf{X}(\nu_1) = \mathbf{X}_1, \quad \mathbf{X}(\nu_f) = \mathbf{X}_f \\
-\overline{\Delta \nu_i I_3} \le \Delta \mathbf{V}_i \le \overline{\Delta \mathbf{v}_i I_3}, \quad \forall \mathbf{i} = \mathbf{1}, \dots, \mathbf{N}
\end{array}
\right.
\end{array}
\qquad (2)
$$

where $\overline{\Delta v_i}$ denote the saturation bounds for the spacecraft thrusters. Note that $\overline{\Delta v_i}$ depend on the impulse time, $v_i$ and on the actuators limitation $\Delta V_{max}$

$$
\overline{\Delta v_i} = \sqrt{\frac{a^3(1-e^2)^3}{\mu} \frac{1}{1+e \cos v_i}} \Delta V_{max} \quad (3)
$$

B. Direct approach: Linear Programming

Direct shooting methods have shown their relevance for solving the rendezvous guidance problem
in the last decade [11, 13, 15, 37, 38]. Their usage are handy since closed form solutions of the
linearized dynamics (1) exist. For instance, Yamanaka and Ankersen developed in [36] a transition
matrix $\Phi_{ya}$ that enables the propagation of the chaser spacecraft relative state starting from an
initial state $X_1(\nu_1)$ under the action of $N$ impulsive controls:

$$X(\nu) = \Phi_{ya}(\nu, \nu_1)X(\nu_1) + \sum_{i=1}^{N} \Phi_{ya}(\nu, \nu_i)B\Delta V_i = \Phi_{ya}(\nu, \nu_1)X(\nu_1) + B\Delta V \quad (4)$$

where $\nu_1 < \nu_2 < \dots < \nu_N \le \nu$. $\Delta V$ denotes the stacked control vector

$$\Delta V = [\Delta V_1^T, \dots, \Delta V_N^T]^T \quad (5)$$

and

$$B = [\Phi_{ya}(\nu, \nu_1)B \quad \dots \quad \Phi_{ya}(\nu, \nu_N)B], \quad \text{with } B = \begin{bmatrix} O_3 \\ I_3 \end{bmatrix} \quad (6)$$

where $O_3 \in \mathbb{R}^{3 \times 3}$ is the null matrix and $I_3 \in \mathbb{R}^{3 \times 3}$ is the identity matrix.

By using the closed form solution provided by the transition matrix, the dynamic program (2) can be transformed into the following static one

$$\begin{array}{ll}
\min_{\Delta V_i, \nu_i} & J(\cdot) \\
\text{under} & \left\{
\begin{array}{l}
\mathbf{X}_f = \Phi_{ya}(\nu_f, \nu_1)\mathbf{X}(\nu_1) + \mathbf{B}\Delta\mathbf{V} \\
\mathbf{X}(\nu_1) = \mathbf{X}_1, \quad \mathbf{X}(\nu_f) = \mathbf{X}_f \\
-\overline{\Delta\nu_i}I_3 \le \Delta\mathbf{V}_i \le \overline{\Delta\mathbf{v}_i}I_3, \quad \forall \mathbf{i} = \mathbf{1}, \dots, \mathbf{N}
\end{array}
\right.
\end{array}
\quad (7)$$

The previous optimal control problem can be recast into a linear program by setting the impulse times $\{\nu_i\}_{i=1,...,N}$ a priori and electing an optimization criteria, $J(\cdot)$, that is linear in the decision variables $\Delta V = \{\Delta V_i\}_{i=1,...,N}^T$. Usually in the rendezvous problem, the criteria represents the fuel consumption. However, it will be showed next that, in the presence of GNC system errors, the choice of the optimization criteria is of importance and needs to be made carefully in order to include the robustness requirements. The nature of the objective and the description of the cost function are discussed further in a dedicated section after presenting the GNC system errors and their impacts on the rendezvous guidance problem.

### C. Effects of the GNC errors on the rendezvous guidance problem

As part of the GNC system, the performances of the guidance algorithms are tightly related to the performances of the navigation devices and propulsion equipment. Errors from these devices
and equipment can greatly deteriorate the guidance precision performances.

The navigation errors affect the initial state $\mathbf{X}_1$ used for the calculation of the guidance plan:

$$ \mathbf{X}_1 \in \mathcal{U}_{\text{nav}} \qquad (8) $$

The uncertainty set $\mathcal{U}_{\text{nav}}$ is directly related to the accuracy of the sensors and to the performances the navigation filter.

The impulses execution errors are twofold: errors on the firing times, $\nu_i$, and errors on the impulses execution. The spacecraft thrusts have been modelled as pure impulsive control, but practical execution of the maneuver make the trust profile looks like a continuous functions. To account for the modelling errors of the impulsive assumption, the firing times, $\nu_i$ can be considered as uncertain:

$$ \nu_i \in [\underline{\nu}_i, \overline{\nu}_i], \quad i = 1, \dots, N \qquad (9) $$

The impulses $\Delta\mathbf{V}_i$, are computed in the leader's LVLH frame, but the global effort is allocated on the different available gas thrusters. The follower's thrusters also need to be aligned before each impulse execution by means of attitude slew maneuvers. Hence, the precision of the control execution depends on the precision of this alignment. Thus, the thrust allocation and the alignment process introduce uncertainties in the impulses execution:

$$ \Delta\mathbf{V}_i \in \mathcal{U}_{\Delta\mathbf{V}_i} \qquad (10) $$

Previously, the rendezvous guidance problem has been defined with strict initial and final objectives [see Fig. 2 and Eq. (7)]. However, when the GNC uncertainties and errors (8) (9), (10) are integrated in the trajectory propagation (4), the final relative state $\mathbf{X}(\nu_f)$ can no longer be exactly computed;

$$ \mathbf{X}(\nu_f) = \Phi_{\text{ya}}(\nu_f, \nu_1)\mathbf{X}(\nu_1) + \sum_{i=1}^{N} \Phi_{\text{ya}}(\nu_f, \nu_i)\mathbf{B}\Delta\mathbf{V}_i \in \mathcal{X}_f, \quad \begin{cases} \mathbf{X}_1 \in \mathcal{U}_{\text{nav}} \\ \nu_i \in [\underline{\nu}_i, \overline{\nu}_i] \\ \Delta\mathbf{V}_i \in \mathcal{U}_{\Delta\mathbf{V}_i} \end{cases} \qquad (11) $$

In this case, the final rendezvous condition must be relaxed to enable for some tolerances on the final objective. The rendezvous terminal condition is replaced by a set membership constraint on
the spacecraft final relative state for all the possible instances of the uncertainties:

$$ \mathbf{X}_f - \mathbf{X}(\nu_f) = \mathbf{X}_f - \mathbf{\Phi}_{\text{ya}}(\nu_f, \nu_1)\mathbf{X}_1 - \mathbf{B}\mathbf{\Delta}\mathbf{V} \in \mathcal{T}, \quad \left\{ \begin{array}{l} \forall \mathbf{X}_1 \in \mathcal{U}_{\text{nav}} \\ \forall \nu_i \in [\underline{\nu}_i, \bar{\nu}_i] \\ \forall \mathbf{\Delta}\mathbf{V}_i \in \mathcal{U}_{\mathbf{\Delta}\mathbf{V}_i} \end{array} \right. \qquad (12) $$

The set $\mathcal{T}$ defines on admissible tolerance around the original rendezvous objective in presence of GNC errors. The rendezvous guidance problem (7) can be amended as the following uncertain program:

$$ \begin{array}{ll} \min_{\mathbf{\Delta}\mathbf{V}} & J(\cdot) \\ \text{under} & \left\{ \begin{array}{l} \mathbf{X}_f - \mathbf{\Phi}_{\text{ya}}(\nu_f, \nu_1)\mathbf{X}(\nu_1) - \mathbf{B}\mathbf{\Delta}\mathbf{V} \in \mathcal{T} \\ (\mathbf{X}_1, \nu, \mathbf{\Delta}\mathbf{V}) \in \mathcal{V}, \\ -\overline{\Delta\nu_i}I_3 \le \mathbf{\Delta}\mathbf{V}_i \le \overline{\mathbf{\Delta}\mathbf{v}_i}I_3, \forall \mathbf{i} = \mathbf{1}, \dots, \mathbf{N} \end{array} \right. \end{array} \qquad (13) $$

where the set of uncertainties, $\mathcal{V}$, as the direct product of the errors from (8) (9) and (10). It should be outlined that the tolerance set, $\mathcal{T}$, can interpreted as an inclusion set of all possible instance of the final state $\mathbf{X}(\nu_f)$ under GNC device errors.

In the sequel, the proposed work is focused on the impulses execution errors i.e. uncertainties on impulses time, $\nu$ and the mis-realization of the impulses, $\mathbf{\Delta}\mathbf{V}$.

### D. Relaxed polytopic rendezvous condition

In the proposed approach, $\mathcal{T}$ is chosen to be a polytopic set defined by the matrices $H_{\mathcal{T}} \in \mathbb{R}^{n_{\mathcal{T}} \times 6}$ and $K_{\mathcal{T}} \in \mathbb{R}^{n_{\mathcal{T}}}$ . The spacecraft relative state at the end of the rendezvous maneuver plan, $\mathbf{X}(\nu_f)$, is required to satisfy the following matrix inequality:

$$ H_{\mathcal{T}}\mathbf{X}(\nu_f) \le \mathbf{K}_{\mathcal{T}}(\mathbf{X}_f) \Leftrightarrow H_{\mathcal{T}}(\mathbf{B}\mathbf{\Delta}\mathbf{V} + \mathbf{\Phi}_{\text{ya}}(\nu_f, \nu_1)\mathbf{X}_1) \le \mathbf{K}_{\mathcal{T}}(\mathbf{X}_f) \qquad (14) $$

The $\mathbf{K}_{\mathcal{T}}$ vector is chosen such that the final set $\mathcal{T}$ is centered around the original rendezvous objective $\mathbf{X}_f$ see figure 3. Thus, the polytopic version of the uncertain optimal control problem (13)
is described by the following program:

$$
\begin{array}{ll}
\min_{\Delta \mathbf{V}} & J(\cdot) \\
\text{under} & H_{\mathcal{T}}(\mathbf{B}\Delta V + \Phi_{ya}(\nu_f, \nu_1)X_1) \le K_{\mathcal{T}}(X_f) \\
& (\nu, \Delta V) \in \mathcal{V}, \\
& -\overline{\Delta v_i} I_3 \le \Delta \mathbf{V}_i \le \overline{\Delta v_i} I_3, \forall \mathbf{i} = \mathbf{1}, \dots, \mathbf{N}
\end{array}
\qquad (15)
$$

<figure><img src="image_4.png" /><figcaption>Fig. 3: Relaxed rendezvous illustration: Polytopic tolerance arrival set</figcaption></figure>

Without loss of generality, the tolerance polytope $\mathcal{T}$ can be specified as a parallelotope with the
defining matrices $H_{\mathcal{T}}$ and $K_{\mathcal{T}}$ given by:

$$
H_{\mathcal{T}} = \begin{bmatrix} \mathbb{I}_6 \\ -\mathbb{I}_6 \end{bmatrix}, \quad \mathbf{K}_{\mathcal{T}} = \begin{bmatrix} \Gamma + \mathbf{X}_f \\ \Gamma - \mathbf{X}_f \end{bmatrix} \quad \text{with } \mathbb{R}^6 \ni \mathbf{\Gamma} > \mathbf{0} \qquad (16)
$$

Specifying a final tolerance in the form (16) presents computational advantages that will be ex-
ploited in what follows. Among these advantages, it enables to reduces the numbers of the defining
parameters (the vector Γ for instance). Moreover, the definition (16) of H_T and K_T defines the set
T as the inclusion set, in the interval analysis sense, for the spread zone of the different errors.

E. The objective function $J(\cdot)$

The typical objective for a rendezvous guidance algorithm is the minimization of the fuel cost of
the computed plan [20]. Assuming that the chaser spacecraft is equipped with six identical thrusters
rigidly mounted in its axes, the total fuel consumption cost is given by [32]:

$$
J = \| \Delta \mathbf{V} \|_{1} \tag{17}
$$
This piecewise linear cost function can be linearized by introducing slack variables $Z \in \mathbb{R}^{3N}$ in the guidance problems (7) and (15) such that

$$ \begin{array}{l} \Delta \mathbf{V}_i \le \mathbf{Z}_i, \quad \forall i = 1, \dots, 3N \\ -\Delta \mathbf{V}_i \le \mathbf{Z}_i \end{array} \qquad (18) $$

and setting $J = \sum_{i=1}^{3N} \mathbf{Z}_i$. The reader will note that:

$$ \min_{\Delta \mathbf{V}} \|\Delta \mathbf{V}\|_1 \Leftrightarrow \min_{\Delta \mathbf{V}, \mathbf{Z}} \sum_{i=1}^{3N} \mathbf{Z}_i \qquad (19) $$

On the other hand, in presence of GNC errors, besides minimizing the fuel cost, certifying the best rendezvous precision is of great interest. This precision objective can be described in terms of the dimensions of the tolerance set $\mathcal{T}$ which contains all the possible final states for the considered uncertainties. By considering the tolerance vector $\Gamma$ as decision variable, the dimensions of the final inclusion set can be minimized while maintaining the linear structure of the cost function. In this case the cost function can be chosen as:

$$ J = \sum_{i}^{6} \Gamma_{i} \qquad (20) $$

By choosing the cost function (20), the objective becomes able to compute the maneuver plan, $\Delta \mathbf{V}$, that is the least sensitive to the considered errors. Once computed, the tolerance set $\mathcal{T}$ gives a certified prediction of the spread zone for all the possible trajectories. It is also interesting to notice that the tolerance set has not to be set a priori but its size is minimized during the process.

Fuel is a critical resource for the spacecraft applications and even when it is not explicitly included in the objective function, it must be ensured that the total consumption does not exceed the allocated budget. For a given fuel budget, $M_{\Delta \mathbf{V}}$, the robust plan must be such that:

$$ \|\Delta \mathbf{V}\|_1 \le \sum_{i=1}^{3N} \mathbf{Z}_i \le M_{\Delta \mathbf{V}} \qquad (21) $$

The cost function (20) will be used in the robust algorithm presented in the next section along with the budget constraint (21) in order to achieve the best rendezvous precision while respecting the required limitations even in presence of errors.

The aim of this work is to obtain a tractable robust counterpart of the uncertain program (13) by using robust optimization concepts (cf appendix A). The uncertainties on impulse time and the
errors affecting the impulses execution are considered separately and for each case a linear program is obtained.

### III. Rendezvous guidance algorithm robust to uncertainties on impulses application time
#### A. Modeling impulses firing time uncertainties using interval analysis

In this section, the maneuver plan will be affected by some uncertainties on the maneuver firing time. The uncertainties are defined as:

$$v_i = v_i^0 + u_i \delta v_i, \quad i = 1, \dots, N \qquad (22)$$

where $v_i^0$ is the nominal firing time, $\delta v_i$ is the maximal disturbance and the perturbation vector $u \in \mathbb{R}^N$ belongs to the normalized interval $[-1, 1]$. The entries of $u$ are supposed to be independent from each other.

Uncertainties on the impulses times have two consequences on the linearized relative motion. First, since the transition matrix $\Phi_{ya}$ depends on the true anomaly, the propagation of the relative trajectory will be affected. Second, the bounds on the control amplitudes, $\Delta v$, also become uncertain since they depend on the position of the target spacecraft on its orbit (cf equation (3)). The present study does not address the uncertainties on the impulses bounds as their effect can be neglected since $\delta v_i$ is relatively small with respect to the orbital period.

In order to build the robust counterpart that accounts for the impulses time uncertainties, it is necessary to evaluate the maximum variation of the elements of the Yamanaka-Ankersen transition matrix on a given impulse time interval. Such information can be obtained through the computation of the interval inclusion function. Let an interval, denoted $[x]$, be a connected and finite set of $\mathbb{R}$ defined by:

$$[x] = [x^{-}, x^{+}] = \{x \in \mathbb{R} | x^{-} \le x \le x^{+}\} \qquad (23)$$

Let mid $([x])$ denote the middle of the interval $[x]$ and $\mathbb{R}$ denote the set of all intervals of $\mathbb{R}$

**Definition 1.** Consider a function $f$ from $\mathbb{R}^n$ to $\mathbb{R}^m$. The interval function $[f]$ from $\mathbb{R}^n$ to $\mathbb{R}^m$ is an inclusion function for $f$ if:

$$\forall [x] \in \mathbb{R}^n, \quad f([x]) \subset [f]([x]) \qquad (24)$$
The inclusion function can be computed using several methods, such as the natural inclusion function, the centered inclusion function or the Taylor inclusion function [39]. To reduce the conservation of the inclusion function, the Centred Inclusion Function (CIF) method is chosen to estimate the transition matrix inclusion denoted $[\Phi_{ya}](\nu_f, [\nu_i])$.

## B. Polytopic rendezvous robust counterpart

The final rendezvous condition is relaxed to polytopic final set to account for the presence of the impulse firing time uncertainties (see Fig. 4).

Under impulse firing time uncertainties, the polytopic rendezvous constraint (14) becomes uncertain since the matrix **B** is unknown and belongs to the interval matrix [**B**] such that

$$[\mathbf{B}] = [ [\Phi_{ya}](\nu_f, [\nu_1])B \dots [\Phi_{ya}](\nu_f, [\nu_N])B ] \quad (25)$$

<figure><img src="image_5.png" /><figcaption>Fig. 4: polytopic rendezvous condition under impulse firing time uncertainties</figcaption></figure>

Thus, the uncertain polytopic rendezvous condition is given by

$$H_{\mathcal{T}} \left( \Phi_{ya}(\nu_f, \nu_1) X(\nu_1) + \sum_{i=1}^{N} \Phi_{ya}(\nu_f, \nu_i) B \Delta \mathbf{V}_i \right) \leq \mathbf{K}_{\mathcal{T}} \quad (26)$$

where

$$[\Phi_{ya}](\nu_f, [\nu_i]) \ni \Phi_{ya}(\nu_f, \nu_i) = \hat{\Phi}_{ya}^0(\nu_f, \nu_i) + \delta\Phi_{ya}(\nu_f, \nu_i) \quad (27)$$

The matrix $\Phi_{ya}^0(\nu_f, \nu_i)$ is the center matrix of the inclusion matrix $[\Phi_{ya}](\nu_f, [\nu_i])$. The matrix $\Phi_{ya}^0(\nu_f, \nu_i)$ is generally different from the nominal transition matrix, $\Phi_{ya}(\nu_f, \nu_i)$. The perturbation matrix $\delta\Phi_{ya}(\nu_f, \nu_i)$ belongs to the interval matrix $[-\overline{\delta\Phi_i}, \overline{\delta\Phi_i}]$, $\overline{\delta\Phi_i}$ being the radius matrix of the inclusion matrix $[\Phi_{ya}](\nu_f, [\nu_i])$. The uncertain set $\mathcal{V}$ can be then expressed as:

$$\mathcal{V} = \{ \hat{\Phi}_{ya}(\nu_f, \nu_i) | \hat{\Phi}_{ya}(\nu_f, \nu_i) = \Phi_{ya}^0(\nu_f, \nu_i) + u_i \overline{\delta\Phi_i}, | u_i | \le 1, i = 1, \dots, N \} \quad (28)$$
According to [34], the robust counterpart of the uncertain rendezvous condition (26) can be written as:

$$H_{\tau} \left( \Phi_{ya}(\nu_f, \nu_1) \mathbf{X}(\nu_1) + \sum_{i=1}^{N} \Phi_{ya}^{0}(\nu_f, \nu_i) \mathbf{B} \Delta \mathbf{V}_i + \max_{\mathbf{u}_i} \{ \mathbf{u}_i \delta \overline{\Phi}_i \mathbf{B} \Delta \mathbf{V}_i \} \right) \leq K_{\tau} \quad (29)$$

which is equivalent to

$$H_{\tau} \Phi_{ya}(\nu_f, \nu_1) \mathbf{X}(\nu_1) + \sum_{i=1}^{N} \mathbf{H}_{\tau} \Phi_{ya}^{0}(\nu_f, \nu_i) \mathbf{B} \Delta \mathbf{V}_i + \mathbf{H}_{\tau} |\delta \overline{\Phi}_i \mathbf{B} \Delta \mathbf{V}_i| \leq \mathbf{K}_{\tau} \quad (30)$$

By observing that the slack variables $Z_i$, such that $|\delta \overline{\Phi}_i \mathbf{B} \Delta \mathbf{V}_i| \leq |\delta \overline{\Phi}_i \mathbf{B}||\Delta \mathbf{V}_i| \leq |\delta \overline{\Phi}_i|Z_i$, the previous inequality can be linearized to obtain the robust program to uncertain impulse times:

$$\begin{array}{l}
\displaystyle \min_{\Delta \mathbf{V}, \Gamma, Z} \sum_{i=1}^{6} \Gamma_i \\
\text{under} \quad \left\{
\begin{array}{l}
H_{\tau} \Phi_{ya}(\nu_f, \nu_1) \mathbf{X}(\nu_1) + \mathbf{H}_{\tau} \sum_{i=1}^{N} \Phi_{ya}^{0}(\nu_f, \nu_i) \mathbf{B} \Delta \mathbf{V}_i + |\delta \overline{\Phi}_i \mathbf{B}| Z_i \leq \mathbf{K}_{\tau} \\
\Delta \mathbf{V}_i \leq \mathbf{Z}_i \\
-\Delta \mathbf{V}_i \leq \mathbf{Z}_i \\
\sum_{i=1}^{3N} Z_i \leq M_{\Delta V} \\
[Z_{3i+1}, Z_{3i+2}, Z_{3i+3}]^T \leq \overline{\Delta v}_{i+1} I_3, \forall i = 0, \dots, N-1
\end{array}
\right.
\end{array}
\quad (31)$$

Note that the problem (31) has a linear structure.

# IV. Rendezvous guidance algorithm robust to errors in impulses execution

## A. Modeling the uncertainties

The performances of the thrust devices are described here in terms of amplitude and orientation precision (see figure 5).

$$\Delta \mathbf{V}_i = (\mathbf{1} + \lambda_i) \mathbf{M}_{ci} \Delta \mathbf{V}_i^0 \quad (32)$$

where the coefficient $\lambda_i$ represents the amplitude error and is such that

$$\lambda_i \in [-\varepsilon, \varepsilon] \quad (33)$$

and the Cardan rotation matrix $M_{ci}(\psi_i, \theta_i, \phi_i)$ represent the orientation error (see figure 6).
<figure><img src="image_6.png" /><figcaption>Fig. 5: Impulse execution errors description</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6: Cardan angles: rotation frome frame $(\vec{x}_0, \vec{y}_0, \vec{z}_0)$ to $(\vec{x}_3, \vec{y}_3, \vec{z}_3)$</figcaption></figure>

As long as the Cardan angles $\psi, \phi, \theta$ remain small, the Cardan rotation matrix can be approximated by

$$ M_{ci} = \begin{bmatrix} 1 & -\psi_i & \theta_i \\ \psi_i & 1 & -\phi_i \\ -\theta_i & \phi_i & 1 \end{bmatrix} \qquad (34) $$

where the Cardan angles are bounded by $\beta$, the maximal error angle

$$ |\psi_i| \le \beta, \quad |\theta_i| \le \beta, \quad |\phi_i| \le \beta \qquad (35) $$

Let $M_i$ denote the perturbation matrix such that

$$ M_i = (1 + \lambda_i) M_{ci} = \begin{bmatrix} 1 + \lambda_i & -(1 + \lambda_i)\psi_i & (1 + \lambda_i)\theta_i \\ (1 + \lambda_i)\psi_i & 1 + \lambda_i & -(1 + \lambda_i)\phi_i \\ -(1 + \lambda_i)\theta_i & (1 + \lambda_i)\phi_i & 1 + \lambda_i \end{bmatrix} \qquad (36) $$

Note that $M_i$ is composed by four independent element, thus, four perturbation variables $u_{ij}$ are
required to describe the affine uncertainty set $\mathcal{V}$.

$$ \mathcal{V} = \{ M_i \mid M_i = \mathbb{I}_3 + \sum_{j=1}^{4} u_{ij} M^j, \ |u_{ij}| \le 1, i = 1, \dots, N \} \quad (37) $$

where

$$ M^1 = \begin{bmatrix} \epsilon & 0 & 0 \\ 0 & \epsilon & 0 \\ 0 & 0 & \epsilon \end{bmatrix}, \quad M^2 = \begin{bmatrix} 0 & -(1+\epsilon)\beta & 0 \\ (1+\epsilon)\beta & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}, \quad M^3 = \begin{bmatrix} 0 & 0 & (1+\epsilon)\beta \\ 0 & 0 & 0 \\ -(1+\epsilon)\beta & 0 & 0 \end{bmatrix}, \quad M^4 = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & -(1+\epsilon)\beta \\ 0 & (1+\epsilon)\beta & 0 \end{bmatrix} \quad (38) $$

## B. Polytopic rendezvous robust counterpart

The rendezvous mission under errors on the impulse execution with a relaxed polytopic arrival set is illustrated in Fig 7. The uncertainties on the impulses execution can be translated into uncertainties on the transition matrix by using the description of the uncertainty set (37). In this case, the uncertain polytopic rendezvous condition can written as:

$$ H_{\tau} \Phi_{ya}(\nu_f, \nu_1) \mathbf{X}(\nu_1) + \mathbf{H}_{\tau} \mathbf{B} \mathbf{M}(\mathbf{u}) \mathbf{\Delta} \mathbf{V} \leq \mathbf{K}_{\tau} \quad (39) $$

where

$$ \mathbf{M}(\mathbf{u}) = \mathrm{diag}(\{M_1, \dots, M_N\}) = \mathbb{I}_{3N} + \underbrace{\begin{bmatrix} \sum_{j=1}^{4} u_{1j} M^j \\ \vdots \\ \sum_{j=1}^{4} u_{Nj} M^j \end{bmatrix}}_{\tilde{M}(\mathbf{u})}, \quad \|\mathbf{u}\|_{\infty} \le 1 \quad (40) $$

and the perturbation variable $\mathbf{u} \in \mathbb{R}^{4N}$ is norm-bounded such that $\|\mathbf{u}\|_{\infty} \le 1$.

<figure><img src="image_8.png" /><figcaption>Fig. 7: Rendezvous affected by impulses execution errors with polytopic arrival set</figcaption></figure>
A robust counterpart can be formulated by applying the results from [40]:

$$H_{\tau}\Phi_{ya}(\nu_f, \nu_i)\mathbf{X}(\nu_1) + \mathbf{H}_{\tau}\mathbf{B}\mathbf{\Delta}\mathbf{V} + \max_{\|\mathbf{u}\|_{\infty} \le 1} \{\mathbf{H}_{\tau}\mathbf{B}\mathbf{\Delta}\mathbf{M}(\mathbf{u})\mathbf{\Delta}\mathbf{V}\} \le \mathbf{K}_{\tau} \quad (41)$$

Then, it comes that

$$\max_{\|\mathbf{u}\|_{\infty} \le 1} \{H_{\tau}\mathbf{B}\mathbf{\Delta}\mathbf{M}(\mathbf{u})\mathbf{\Delta}\mathbf{V}\} = \max_{\|\mathbf{u}\|_{\infty} \le 1} \left\{ \sum_{i=1}^{N} H_{\tau}\Phi_{ya}(\nu_f, \nu_i)\mathbf{B}\left(\sum_{j}^{4} u_{ij}M^{j}\right)\mathbf{\Delta}\mathbf{V}_{i} \right\} \quad (42)$$

$$= \sum_{i=1}^{N} \sum_{j=1}^{4} \max_{\|\mathbf{u}\|_{\infty} \le 1} \{H_{\tau}\Phi_{ya}(\nu_f, \nu_i)BM^{j}\mathbf{\Delta}\mathbf{V}_{i}\mathbf{u}_{ij}\} \quad (43)$$

$$= \sum_{i=1}^{N} \sum_{j=1}^{4} |H_{\tau}\Phi_{ya}(\nu_f, \nu_i)BM^{j}\mathbf{\Delta}\mathbf{V}_{i}| \quad (44)$$

The robust counterpart to the relaxed guidance problem (13) under uncertainties on impulses execution (32) is finally obtained:

$$\min_{\mathbf{\Delta}\mathbf{V}, \mathbf{\Gamma}} \sum_{i=1}^{6} \Gamma_i \left\{ \begin{array}{l} \Phi_{ya}(\nu_f, \nu_i)\mathbf{X}(\nu_1) + \mathbf{H}_{\tau}\mathbf{B}\mathbf{\Delta}\mathbf{V} + \sum_{i=1}^{N} \sum_{j=1}^{4} |\mathbf{H}_{\tau}\Phi_{ya}(\nu_f, \nu_i)BM^j\mathbf{\Delta}\mathbf{V}_i| \le \mathbf{K}_{\tau} \\ \mathbf{X}(\nu_1) = \mathbf{X}_1 \\ \|\mathbf{\Delta}\mathbf{V}\|_1 \le \mathbf{M}\mathbf{\Delta}\mathbf{V} \\ -\overline{\Delta v_i}I_3I_3 \le \mathbf{\Delta}\mathbf{V}_i \le \overline{\mathbf{\Delta}\mathbf{v}_i}\mathbf{I}_3, \forall \mathbf{i} = \mathbf{1}, \dots, \mathbf{N} \end{array} \right. \quad (45)$$

By using the slack variables Z defined by (18) and the variables $Y \in \mathbb{R}^{4N}$ such that

$$-Y_{ij} \le H_{\tau}\Phi_{ya}(\nu_f, \nu_i)BM^j\mathbf{\Delta}\mathbf{V}_i \le \mathbf{Y}_{ij}, \quad \mathbf{i} = \mathbf{1}, \dots, \mathbf{N}, \quad \mathbf{j} = \mathbf{1}, \dots, \mathbf{4} \quad (46)$$

the problem (45) can be reformulated as the following linear program

$$\min_{\mathbf{\Delta}\mathbf{V}, \mathbf{\Gamma}, \mathbf{Z}, \mathbf{Y}} \sum_{i=1}^{6} \Gamma_i \left\{ \begin{array}{l} \Phi_{ya}(\nu_f, \nu_i)\mathbf{X}(\nu_1) + \mathbf{H}_{\tau}\mathbf{B}\mathbf{\Delta}\mathbf{V} + \sum_{i=1}^{N} \sum_{j=1}^{4} \mathbf{Y}_{ij} \le \mathbf{K}_{\tau} \\ -Y_{ij} \le H_{\tau}\Phi_{ya}(\nu_f, \nu_i)BM^j\mathbf{\Delta}\mathbf{V}_i \le \mathbf{Y}_{ij} \\ -\mathbf{Z}_i \le \mathbf{\Delta}\mathbf{V}_i \le \mathbf{Z}_i, \quad \text{, } i=1, \dots, N, j=1, \dots, 4 \\ \sum_{i=1}^{3N} \mathbf{Z}_i \le \mathbf{M}\mathbf{\Delta}\mathbf{V} \\ [Z_{3i+1}, Z_{3i+2}, Z_{3i+3}]^T \le \overline{\Delta v_{i+1}}I_3, \forall i=0, \dots, N-1 \end{array} \right. \quad (47)$$
V. Numerical examples

<table><caption>Table 1: Rendez-vous missions repertory</caption><thead><tr><th>Missions parameters</th><th>ATV 1</th><th>ATV 2</th><th>PROBA 3</th></tr></thead><tbody><tr><td>Semi-major axis, a [km]</td><td>6763</td><td>6763</td><td>37039.887</td></tr><tr><td>Eccentricity</td><td>0.0052</td><td>0.0052</td><td>0.80621</td></tr><tr><td>Inclination [degree]</td><td>52</td><td>52</td><td>60.15</td></tr><tr><td>Argument of Perigee [degree]</td><td>0</td><td>0</td><td>180.6</td></tr><tr><td>Right Ascension of the Ascending Node [degree]</td><td>0</td><td>0</td><td>173</td></tr><tr><td>Saturation [m/s]</td><td>5</td><td>5</td><td>0.8</td></tr><tr><td>Initial anomaly [degree]</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Mission duration [s]</td><td>2767</td><td>7200</td><td>141888</td></tr><tr><td>Initial relative state [m,m/s]</td><td>-30000<br>0<br>5000<br>8.154<br>0<br>0</td><td>-30000<br>0<br>5000<br>8.154<br>0<br>0</td><td>-5000<br>0<br>0<br>0<br>0<br>0</td></tr><tr><td>Final relative state [m,m/s]</td><td>-1000<br>0<br>0<br>0<br>0<br>0</td><td>-1000<br>0<br>0<br>0<br>0<br>0</td><td>-20<br>0<br>0<br>0<br>0<br>0</td></tr></tbody></table>

The optimal deterministic guidance algorithm (7) minimizing the ΔV-consumption objective (17) on one hand and the robust algorithms (31) and (45) on the other hand are compared on the benchmark missions from Table 1. The proposed rendezvous missions are chosen to cover low and high eccentricity target orbits and a large range of mission durations. The missions from Table 1 will be simulated in a Keplerian and linear environment first. Then, simulations are run in a nonlinear environment for a specific Prisma type mission to validate the approach described in this paper.
Most of the results are obtained considering a number of impulses, $N$, set to 50, equally distributed on the mission time horizon. Note that the influence of the number of impulses is also studied as well. It should be noted that all the simulations presented here are run in an open-loop fashion.

For the robust algorithms (31) and (45), the allocated $\Delta V$ budget, $M_{\Delta V}$, is calculated with respect to $J_{\Delta V}^*$, the optimal $\Delta V$-consumption obtained by solving (7) along with the cost function (17):

$$M_{\Delta V} = \alpha J_{\delta V}^* \qquad (48)$$

The aim of comparison is twofold. First, the numerical simulation must validate that the computed inclusion set $\mathcal{T}$ contains the spread of errors included in the defined uncertainty sets. The second aim of this numerical study is to ascertain the gain of the robust open-loop approach with respect to the deterministic algorithm (7) in terms of errors propagation by deterministic and robust maneuver plans and in terms of $\Delta V$ consumption. To make this study, for each case of uncertainties for tested mission, several graphics are proposed. On one hand figures 8, 10, 13, 15, 18, 19 propose four plots. Two plots present the fuel-optimal and the robust maneuver plans along $X_{lvlh}$ and $Z_{lvlh}$ axis. Two plots depict the results of Monte-Carlo simulations: the spread of the final relative states obtained by the propagation of the different errors by the fuel-optimal plan and the unrestricted-budget robust plan is exposed. On figures 9, 11, 12, 14, 16, 17, bar diagrams are presented. On this bar diagram, the level of guaranteed tolerance and spread are exposed for different fuel budget. For the tolerance this level corresponds to the robust cost (20). For the "Errors spread" bars, it is the sum of the widths of the spread along each axis. "Errors spread" bars can only be estimated from the Monte-Carlo runs contrary to the computed minimal tolerances. In sake of comparison, two information are added to the bar graphics: the level of errors spread for the fuel-optimal plan is indicated through the light grey line and the consumed $\Delta V$ is given by the black line.

## A. Uncertainties affecting the impulse firing time: numerical examples

In order to illustrate the efficiency of the algorithm (31), simulations are run considering random disturbances on each nominal firing location. Those disturbances are bounded within $\delta\nu_i$ such that it corresponds to 1 second errors.
In the case of the ATV 2 (Automated Transport Vehicle) mission (figures 8), the fuel-optimal plan consists of 4 impulses out of the possible 50, the robust plans are also of 4 impulses. The research of the robustness with respect to firing location errors is achieved by shifting the two first impulses to the left. In the robust plan the second and third impulses consume most part the fuel while the first and third ones are the most demanding for the fuel-optimal plan. One can also note that the impulses are only on $X_{lvlh}$-axis thrusts for the three plans.

On figures 8, one can observe that the precision of the mission is improved by the robust algorithm as expected after the computation of the optimal tolerance. The improvements are mainly on $Z_{lvlh}$ axis on position and $X_{lvlh}$ axis on velocity. The bar graphic 9 exposed the level of guaranteed tolerance and the errors spread for different budgets. Different facts can be observed from 9. First, a budget of 100% of the fuel-optimal consumption is needed to guarantee a better precision than the spread by the fuel optimal plan of the firing time errors. In fact, the level of spread of errors by the non robust plan is 10% larger than the guaranteed tolerance while the actual level of spread by the robust plan is 10% smaller. Moreover, when the budget is not restricted, it can be seen that the consumption is very similar to the minimal fuel cost $J_{\Delta V}^*$. This fact infers that the robustness do not necessarily imply an increase of the $\Delta V$ consumption but it is also a matter of impulsive thrusts location.

the advantages of using robust algorithm are more obvious in the case of PROBA 3 mission. Figures 10 show that the final states obtained by the propagating of the optimal plan are largely spread over 1800 m on Z axis and 5 m/s on $V_x$ axis while the optimal certified precision is about of 1 m along $X_{lvlh}$ and 7 m along $Z_{lvlh}$ in position and 20 mm/s and 1 mm/s in velocity. These results are obtained at the price of an 200% increase of the $\Delta V$ consumption.

However, for this difficult mission, improvement, with respect to the spread of the fuel-optimal plan, can be obtained with less than the fuel-optimal cost. In fact, it can be seen on figure 11 that only 90% of the nominal optimal consumption, $J_{\Delta V}^*$, is needed to certify a much better inclusion set and to greatly reduce the spread of the final disturbed states. Note that the level line of the non robust spread does not appear on figure 11 since the level of the spread is much larger than the worst computed tolerance. This fact is symbolized by the triangle. With a budget of 200%
of the optimal consumption, the guaranteed tolerance dimension are improved of 35%. Beyond
a $\Delta V$ budget of 200% of the optimal consumption, the guaranteed tolerance is not improved
significantly even if more $\Delta V$ consumption permits to tighten more the spread of errors. For both
missions, one can note that the gap between the guaranteed tolerance and the actual spread of the
errors is about 10 points of percentage for budget greater than 100% of the optimal-fuel consumption.

For the previous numerical experiments, the number of possible impulses is fixed to 50. In order
to evaluate the influence of $N$, we reproduce the experiments for different numbers of impulses,
ranging from 10 to 200 with an unrestricted budget (see figures 12). In the case of the ATV
2 mission, increasing the number of impulses only brings small improvements of the guaranteed
tolerance and errors spread while the consumption keeps unchanged. On the contrary, in the case
of PROBA 3 mission, the number of impulses has a stronger influence. The guaranteed tolerance
and the spread improve as $N$ increases but, at the price of more consumption.
<figure><img src="image_9.png" /><figcaption>Fig. 8: Mission ATV 2: Rendezvous under uncertain impulses time</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 9: Mission ATV 2: Rendezvous under uncertain impulses time</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 10: Mission PROBA 3: Rendezvous under uncertain impulses time</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 11: Mission PROBA 3: Rendezvous under uncertain impulses time</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 12: Rendezvous under uncertain impulses time with varying number of impulses $N$</figcaption></figure>

## B. Rendezvous maneuver plan robust to the impulses execution errors: numerical examples

To illustrate the properties of the robust guidance algorithm (45), amplitude errors $\epsilon = 0.001$ i.e. 0.1% and orientation errors $\beta = 1^\circ$ are considered on each executed impulsive maneuver. To simulate execution errors on each impulses, the angles $\psi_i$, $\theta_i$, $\phi_i$ and amplitude coefficient $\lambda_i$ are chosen randomly within their bound, then the computed $\Delta V_i$ is perturbed accordingly to (32). The results only present the in-plane propagation. However, one must keep in mind that mis-execution errors are spread over the three directions.

In the case of ATV 1 mission, figures 13 show that the minimal certified tolerances ensure much better precision in position than the one can be expected with the $\Delta V$ optimal plan. Moreover, the
spread of the errors is much more contained with the robust plan.

In Fig. 13, the optimal 4-impulses plan takes advantage of a long coasting period between the second and third impulse to reduce consumption of the ATV 1 mission. The unrestricted-budget robust maneuver plan is composed of 4 impulses. The most consuming impulsion of this robust plan is the second one that is placed halfway to the end of the mission. When limiting the budget, the robust strategy is slightly different: the halfway impulse is replaced by two impulses on $X_{lvlh}$ axis. However, it can be noticed that, when the robustness is sought, it is preferable to to thrust at midway to the end of the mission than having a coasting period.

Some information can be deduced from figure 14. First, the level of guaranteed tolerance is larger than the non robust plan spread of the errors, some improvement can be certified on given directions (for instance for the position). However, the spread is significantly smaller for robust plan than for optimal-fuel plan, even if it can not be certified for this particular mission. This fact is explained by the conservatism observed on figures 13, can also be highlighted in bar graph 14 since for every budget, there is a large gap between the guaranteed tolerance and the errors spread. This gap is between 65 and 35 points of percentage for budget greater than 100% of the fuel-optimal consumption.

Figures 15 and 16 present more spectacular results for the PROBA 3 mission since the fuel-optimal plan spread the impulses execution errors over nearly 30km along $X_{lvlh}$ axis while the robust plan certify tolerance width of 750m and limit the actual spread of the errors to about 320m on the same axis. The same phenomenon can be observed on the velocity where the fuel-optimal plan spread over 17m/s along the Z axis while the certified tolerance is about 0.5m/s and the errors are spread on 20 mm/s. Observing figures 15, one can note that specific thrust instants are to be preferred when robustness is sought for this particular mission. In addition, these robust thrust locations are independent from the allocated budget.

Reminding that the triangle in figure 16 means that the non robust spread of much larger than the certified tolerances, the robust plans permits a significant reduction of the maneuver plan sensitivity with respect to impulse execution errors. However, the dimension of the minimal
tolerance dimension shows large conservatism. The gap between the level guaranteed tolerance and
the actual spread is of between 90 and 40 % for the presented budgets.

Finally, the influence of the impulse number is analyzed. the previous experiments are repro-
duced assuming different numbers of impulses and with an unrestricted budget. For each Monte-
Carlo simulations run, the guaranteed tolerance, the spread of the errors and the consumption are
represented on Figures 17. On the other hand, in the case of errors on impulse firing time, increas-
ing the number of impulses, for the ATV 1 mission case, permits the improvement of the tolerance
and the errors spread. Moreover, it also slightly improves the robust consumption through a better
selection of impulses locations.

C. Results analysis

The both tested algorithms show the advantage to significantly reduce the spread of the errors
compared to the fuel-optimal plan. If this improvement can be certified a priori in the case of errors
on firing time, this is not always the case for impulses mis-execution errors. In fact, for firing time
error case, the computed tolerance sets are a pretty good inclusion set of the spread of the firing
time errors by the robust plans. This fact is due to the fact that the computation of the transition
inclusion matrix [Φ<sub>ya</sub>(ν<sub>f</sub>, [ν<sub>i</sub>])] is pretty accurate using the centered inclusion function [39] for errors
corresponding to 1 second (δν is about few milliradians). On the contrary, the propagation of the
mis-execution errors, represented by matrix M through the transition matrices Φ<sub>ya</sub>(ν<sub>f</sub>, ν<sub>i</sub>) imply
wrapping effects. The wrapping effect corresponds here to the artificial enlargement of the inclusion
that includes the image of an interval vector through matrix multiplication [39]. This effect is
represented in the equation (46), where the slack variables Y represent the inclusion set of the
spread of errors. These effect imply that the robust algorithm will seek to certify the rendezvous
precision for larger errors than needed.
<figure><img src="image_14.png" /><figcaption>Fig. 13: Mission ATV 1: Rendezvous under impulses execution errors</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 14: Mission ATV 1: Rendezvous under impulses execution errors</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 15: Mission PROBA 3: Rendezvous under impulses execution errors</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 16: Mission PROBA 3: Rendezvous under impulses execution errors</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 17: Rendezvous under impulses execution errors with varying number of impulses *N*</figcaption></figure>

D. Nonlinear simulations

The previously presented results have been obtained for the linearized spacecraft relative dy-
namics propagated using the Yamanaka and Ankersen transition matrix. This section illustrates
the performances of the robust algorithms when nonlinear dynamics are simulated instead.

The nonlinear simulator relies on the nonlinear propagation through Gauss equation integration
of each spacecraft and on the evaluation of the relative motion. This is illustrated for a PRISMA mis-
sion that is detailed in table 2.

A short mission duration has been chosen in order to be able to simulate the nonlinear dynamics
and orbital disturbances (J₂ term and atmospheric drag mainly) while maintaining the difference
<table><caption>Table 2: Mission PRISMA type</caption><tr><td>Semi-axis a, [km]</td><td>7011</td></tr><tr><td>Eccentricity e</td><td>0.004</td></tr><tr><td>Inclination [degree]</td><td>98</td></tr><tr><td>Right Ascension of the Ascending node [degree]</td><td>190</td></tr><tr><td>Argument of Perigee [degree]</td><td>0</td></tr><tr><td>Initial target anomaly [degree]</td><td>0</td></tr><tr><td>Saturation [N]</td><td>0.26</td></tr><tr><td>Mission duration [s]</td><td>1500</td></tr><tr><td>Initial relative state [m,m/s]</td><td>[300, 0, 0, 0, 0]</td></tr><tr><td>Final relative state [m,m/s]</td><td>[20, 0, 0, 0, 0]</td></tr></table>

between linear and nonlinear propagation within acceptable bounds. If the computed plans are implemented in an open-loop structure for this study, it has to be noted that for longer missions, a closed-loop control law has to be implemented to alleviate the modelling errors and external disturbances. For each type of uncertainties, the robust and the deterministic maneuver plans are tested by means of Monte-Carlo nonlinear simulations runs. The impulse time errors is equivalent to 1 second and the control execution errors is given by $\epsilon = 0.1\%$, $\beta = 1^{\circ}$.

The results obtained for the guidance algorithm robust to impulse trigger time uncertainties when using the nonlinear dynamics are presented in Fig. 18. The robust certificates provided by the rendezvous algorithm are validated also in a nonlinear environment while the non robust plan transgress the certified tolerance bound. The robust plan has also the advantages of only consuming 0.25% more than the optimal plan.

It can be observed in Fig. 19 that the maneuver plan robust to execution errors respects the arrival tolerance set even in the presence of environmental perturbations, while the nonrobust plan violates the tolerance box in position. The price of the certified robustness here is a increase of 8% of the $\Delta V$ consumption.
<figure><img src="image_19.png" /><figcaption>Fig. 18: Mission PRISMA : Rendezvous under uncertain impulsion times</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig. 19: Mission PRISMA: Rendezvous under impulsion execution errors</figcaption></figure>

VI. Conclusion

In this paper, the design of robust algorithms for rendezvous guidance problem under uncer-
tainties due to propulsion system errors is addressed. Tractable algorithms are developed on the
theoretical foundation of the convex robust optimization. In fact, polytopic tolerance arrival set
and interval descriptions of these errors leads to a linear description of the robust counterpart for
the optimal guidance problem. Numerical examples demonstrate the ability of the robust algorithm
to preclude the spread of the errors from propulsion systems even under fuel budget constraint.
Numerical simulations of several missions have been conducted in linear and non linear propagator.
The computed tolerance, the actual spread have been analysed with respect to ΔV consumption budget. It has been observed that minimizing the arrival tolerance set clearly changes the structure of the maneuver plans. Optimizing the guaranteed precision generally implies an increase of the ΔV cost. Most results show that, for a same level of consumption, a robust structure permits to mitigate the spread of the errors in comparison with the non robust optimal plan. Moreover, in some cases, the robust rendezvous algorithms are able to certify a better precision with a ΔV budget equivalent to the one obtained with the nominal plan.

A limitation of the proposed method comes from the fact that the affine description of the uncertainty sets, which the algorithms rely on, are obtained using a rather crude interval analysis, generating a possible high level of conservatism of the description. Thus, applying the proposed methodology to the case where the thrust errors are combined, leads to poor results that are too conservative. The computed tolerance boxes are unsound with respect to the rendezvous objectives. Future works will focus on developing tighter affine descriptions of the uncertainty set, thus minimizing the possible conservatism of the planning algorithm.

## APPENDIX A: ROBUST CONVEX OPTIMIZATION

In optimization framework, an uncertain program (A1) can be defined as follows [40]:

$$ \begin{array}{ll} \min_{\mathbf{x}} & f_0(\mathbf{x}) \\ \text{under} & f_i(\mathbf{x}, \mathbf{\Theta}_i(\mathbf{u}_i)) \le 0, \quad \forall \mathbf{\Theta}_i \in \mathcal{V}_i, \quad \mathbf{i} = 1, \dots, m \end{array} \qquad (A1) $$

where $\mathbf{x} \in \mathbb{R}^n$ denotes the vector of decision variables. The structure of the program is supposed to be fixed and described by the cost function $f_0$ and the constraints applications $f_i$. The data of the problem, $\mathbf{\Theta} \in \mathbb{R}^m$ is uncertain and depends on a disturbance variable $\mathbf{u}$ such that each entry $\mathbf{\Theta}_i$ evolve in uncertainty set $\mathcal{V}_i$:

$$ \mathcal{V}_i = \{\mathbf{\Theta}_i | \mathbf{\Theta}_i = \Pi(u_i), u_i \in \mathcal{U}_i \subset \mathbb{R}^{m_i}\} \qquad (A2) $$

The set $\mathcal{V}$ is the image of the disturbance set $\mathcal{U}$ through application $\Pi(\cdot)$. The paradigm of the robust convex optimization is based on the following statements [33]. First, $x$ must be obtained by solving (A1) without the exact knowledge of the data $\mathbf{\Theta}_i$. Second, the results are valid for any $\mathbf{\Theta}_i \in \mathcal{V}_i$. Third, for a given robust feasible solution $x$, the constraints $f_i(x, \mathbf{\Theta}_i(u_i)) \le 0$ are not
violated for all $\Theta_i \in \mathcal{V}_i$. This leads to the definition of a robust feasible solution to the uncertain program (A1).

**Definition 2.** The decision variable $x$ is said to be robustly feasible if and only if $f_i(x, \Theta_i) \le 0$ for all $\Theta$'s instances inside $\mathcal{V}$.

**Definition 3.** Let $x$ be a given robustly feasible solution. The guaranteed cost for the robust feasible solution $x$ in a worst-case cost sense is obtained by solving the maximization problem

$$ \max_{\Theta_i} \{f_0(\mathbf{x}) : \Theta_i \in \mathcal{U}_i, \forall \mathbf{i}\} \qquad (A3) $$

Thus, the optimal solution of the uncertain problem (A1) is obtained by solving a min-max problem, the so-called robust counterpart:

$$ \min_{\mathbf{x}} \max_{\Theta_i \in \mathcal{V}_i} f_0(x) \qquad (A4) \\ \text{under } f_i(x, \Theta_i(u_i)) \le 0, \quad \forall \Theta_i \in \mathcal{U}_i, i = 1, \dots, m $$

Robust convex programming aims at achieving a tractable description of the robust counterpart (A4) [33]. In the specific case of linear programming,

$$ \min_{\mathbf{x}} \max_{(\mathbf{A}, \mathbf{b}) \in \mathcal{U}} \gamma^T \mathbf{x} \\ \text{under } \mathbf{A}\mathbf{x} \le \mathbf{b}, \forall (\mathbf{A}, \mathbf{b}) \in \mathcal{U} $$

the convexity along with the tractability of the robust counterparts (A4) is ensured by considering only uncertainty sets affine in the disturbance variable

$$ \mathcal{V} = \left\{ [\mathbf{A}; \mathbf{b}] = [\mathbf{A}^0, \mathbf{b}^0] + \sum_{j=1}^{k} \mathbf{u}_j [\mathbf{A}^j, \mathbf{b}^j], \mathbf{u} \in \mathcal{V} \subset \mathbb{R}^k \right\} \qquad (A5) $$

Thus, the nature of $\mathcal{V}$ is completely described by the geometry of the $\mathcal{U}$ set. Convex robust counterparts have been developed for specific disturbance sets $\mathcal{U}$. One case of particular interest for this work is the interval case, $\|\mathbf{u}\|_{\infty} \le 1$ for which the robust counterpart is also a linear program [34].

## APPENDIX B: INCLUSION FUNCTIONS

Consider a function $\mathbf{f}$ from $\mathbb{R}^n$ to $\mathbb{R}^m$. The interval function, $[\mathbf{f}]$, from $\mathbb{R}^n$ to $\mathbb{R}^m$ is an inclusion function for $\mathbf{f}$ if

$$ \forall [\mathbf{x}] \in \mathbb{R}^n, \quad \mathbf{f}([\mathbf{x}]) \subset [\mathbf{f}]([\mathbf{x}]) $$
To illustrate the notion of inclusion function, consider a function $f$ from $\mathbb{R}^2$ to $\mathbb{R}^2$, with variables $x_1$ and $x_2$ that vary within intervals $[x_1]$ and $[x_2]$. As illustrated on figure 20, the interval vector $[f]([\mathbf{x}])$ is an inclusion function $[f]$ of $f$ since it guarantees to contain $f([\mathbf{x}])$. A inclusion function is not unique and its minimality depends on the computation technique.

<figure><img src="image_21.png" /><figcaption>Fig. 20: Image of a box by a vector function $f$ and two of its inclusion functions $[f]$ and $[f]^*$; $[f]^*$ is minimal</figcaption></figure>

## 1. Natural inclusion functions

To build an inclusion function, consider a function

$$f: \mathbb{R}^n \to \mathbb{R}, \quad (x_1, \dots, x_n) \mapsto f(x_1, \dots, x_n),$$

expressed as a finite composition of the operators $+$, $-$, $*$, $/$ and elementary functions (sin, cos, exp, sqrt, ...). An inclusion monotonic and thin inclusion function $[f] : \mathbb{R}^n \to \mathbb{R}$ for $f$ is obtained by replacing each real variable $x_i$ by an interval variable $[x_i]$ and each operator or function by its interval counterpart. This function is called the *natural inclusion function* of $f$. If $f$ involves only continuous operators and continuous elementary functions, then $[f]$ is convergent. If, moreover, each of the variables $(x_1, \dots, x_i)$ occurs at most once in the formal expression of $f$ then $[f]$ is *minimal*.

Natural inclusion functions are not minimal in general, because of the dependency and wrapping effects. The accuracy of the resulting interval strongly depends on the expression of $f$.

## 2. Centred inclusion functions

Let $f: \mathbb{R}^n \to \mathbb{R}$ be a scalar function of a vector $\mathbf{x} = (x_1, \dots, x_n)^T$. Assume that $f$ is differentiable over the domain given by the interval vector $[\mathbf{x}]$, and denote mid $([\mathbf{x}])$ by $\mathbf{m}$. The mean-value
theorem then implies that

$$ \forall \mathbf{x} \in [\mathbf{x}], \exists \mathbf{z} \in [\mathbf{x}] \text{ such that } f(\mathbf{x}) = f(\mathbf{m}) + \mathbf{g}^T(\mathbf{z})(\mathbf{x} - \mathbf{m}), $$

where $\mathbf{g}$ is the gradient of $f$, Thus

$$ \forall \mathbf{x} \in [\mathbf{x}], f(\mathbf{x}) = f(\mathbf{m}) + [\mathbf{g}^T]([\mathbf{x}])(\mathbf{x} - \mathbf{m}), $$

where $[\mathbf{g}^T]$ is an inclusion function for $\mathbf{g}^T$, so

$$ f([\mathbf{x}]) \subseteq f(\mathbf{m}) + [\mathbf{g}^T]([\mathbf{x}])([\mathbf{x}] - \mathbf{m}). $$

Therefore, the interval function

$$ [f_c](\mathbf{x}) \triangleq \mathbf{f}(\mathbf{m}) + [\mathbf{g}^T]([\mathbf{x}])([\mathbf{x}] - \mathbf{m}) $$

is an inclusion function for $f$, which shall call the *centred inclusion function*. To illustrate the interest of this function in the one-dimensional case, consider the function $[f_c](\mathbf{x})$ from $\mathbb{R}$ to $\mathbb{IR}$ defined by

$$ [f_c](\mathbf{x}) \triangleq \mathbf{f}(\mathbf{m}) + [\mathbf{f}']([\mathbf{x}])(\mathbf{x} - \mathbf{m}) $$

for any given $[\mathbf{x}]$. This function can be viewed as affine in $\mathbf{x}$ with an uncertain slope belonging to $[f']([\mathbf{x}])$. The graph of $[f_c](\mathbf{x})$ can be represented by a cone with centre $(\mathbf{m}, \mathbf{f}(\mathbf{m}))$ as illustrated on figure 21. It can be noticed that the smaller the width $w([\mathbf{x}])$ is, the better the cone approximates the function.

<figure><img src="image_22.png" /><figcaption>Fig. 21: Interpretation of the centred inclusion function</figcaption></figure>
When the width of $[x]$ is small, the effect of the pessimism possibly resulting from the interval evaluation of $[g]([x])$ is reduced by the scalar product with $[x] - m$, which is a small interval centred on zero.

ACKNOWLEDGMENTS

The presented works is part of project on optimal rendezvous guidance problems funded by the CNES Research and Technology Project grant R-S07/VF-0001-065. This project involves engineers from CNES (french space agency), EADS Astrium and academic researchers from LAAS-CNRS. The authors would like to thank J.C. Berges from the CNES who provided us with the missions scenarios studied in the paper.

REFERENCES

[1] Lawden, D. F., "General theory of optimal rocket trajectories," *Optimal trajectories for space navigation*, Butterworth, London, England, 1963, pp. 54-78. doi: 10.1016/0032-0633(64)90149-7

[2] Prussing, J. E., "Optimal multiple-impulse orbital rendezvous," Ph. D. Dissertation, Dept. of Aeronautics and Astronautics, Massachusetts Institute of Technology, Cambridge, MA, 1967.

[3] Jezewski, D. J., "Primer vector theory applied to the linear relative-motion equations," *Optimal control applications and methods*, Vol. 1, No. 4, 1980, pp. 387-401. doi: 10.1002/oca.4660010408

[4] Carter, T. E., "Optimal impulsive space trajectories based on linear equations," *Journal on Optimization Theory and Applications*, Vol. 70, No. 2, pp. 277-297. doi: 10.1007/BF00940627

[5] Carter, T. E., "Necessary and sufficient conditions for optimal impulsive rendezvous with linear equations of motion," *Journal of Guidance, Control and Dynamics*, Vol. 10, No. 3, 2000, pp. 219-227. doi: 10.1023/A:1008376427023

[6] Arzelier, D., Kara-Zaitri, M., Louembet, C. and Delibasi, A. "Using polynomial optimization to solve the fuel-optimal linear impulsive rendezvous problem," *Journal of Guidance, Control and Dynamics*, Vol. 34, No. 5, 2011, pp 1567-1572. doi: 10.2514/6.2010-7595

[7] Hull, D. G., "Conversion of optimal control problems into parameter optimization problems," *Journal of Guidance, Control and Dynamics*, Vol. 20, No. 1, Jan-Feb 1997, pp. 57-60. doi:10.2514/2.4033

[8] Waespy, C. M., "Linear-Programming Solutions for Orbital-Transfer Trajectories," *Operations Research*, Vol. 18, No. 4, 1970, pp. 635-653.
[9] Robertson, A., Inalhan, G., How, J., "Formation control strategies for a separated spacecraft inter-
ferometer" *Proceedings of the 1999 American Control Conference*, San Diego, CA, Vol 6, 1999, pp
4142-4147. DOI:10.1109/ACC.1999.786331

[10] Deaconu, G., Louembet, C. and Théron, A. "Designing Continuously Constrained Spacecraft Relative Trajectories for Proximity Operations". *Journal of Guidance, Control, and Dynamics*, (2014). DOI: 10.2514/1.G000283, published online on 4 Jun 2014.

[11] Villerson, M., "Coordination and control of multiple spacecraft using convex optimization techniques,"
Master's Thesis, Massachusetts Institute of Technology, Cambridge, MA, 2002.

[12] Kuwata, Y., Richards, A., Schouwenaars, T., How, J.P., "Distributed Robust Receding Horizon Control for Multivehicle Guidance," *Control Systems Technology, IEEE Transactions on*, Vol. 15, No. 4, July 2007, pp. 627-641. doi: 10.1109/TCST.2007.899152,

[13] Richards, A., Schouwenaars, T., How, J., and Feron, E., "Spacecraft trajectory planning with avoidance constraints using mixed-integer linear programming," *Journal of Guidance, Control, and Dynamics*, Vol. 25, No. 4, 2002, pp 755-764. doi:10.2514/2.4943

[14] Kim, Y., Mesbahi, M., and Hadaegh, F. Y., "Multiple-spacecraft reconfiguration through collision avoidance, bouncing, and stalemate," *Journal of Optimization Theory and Applications*, Vol. 122, No. 2, 2004, pp 323-343. doi:10.1023/B:JOTA.0000042524.57088.8b

[15] Mueller, J., and Larsson, R., "Collision avoidance maneuver planning with robust optimization,"
*Proceedings of the 7th International Conference on Guidance, Navigation and Control Systems*, European
Space Agency, Tralee, Count Kerry, Ireland, June 2008.

[16] Richards, A., and How, J., "Robust variable horizon model predictive control for vehicle maneuvering,"
*International Journal of Robust and Nonlinear Control*, Vol 16, No. 7, 2006, pp 333-351. DOI:
10.1002/rnc.1059

[17] Lu, P. and Liu, X., "Autonomous Trajectory Planning for Rendezvous and Proximity Operations by Conic Optimization," *Journal of Guidance, Control, and Dynamics*, Vol. 36, No. 2, 2013, pp. 375-389, DOI: 10.2514/1.58436

[18] Jifuku, R. and Ichikawa, A., "Optimal pulse strategies for relative orbit transfer along a circular orbit,"
*Journal of Guidance, Control, and Dynamics*, Vol. 34, No. 5, 2011, pp. 1329-1341, DOI: 10.2514/1.51230

[19] Vazquez, R., Gavilan, F. and Camacho E., "Trajectory Planning for Spacecraft Rendezvous with On/Off Thrusters," *Proceedings of the World Congress*, Vol.18, No. 1, pp 8473-8478, 2011, DOI: 10.3182/20110828-6-IT-1002.02445

[20] Fehse, W. *Automated Rendezvous and Docking of Spacecraft*. Cambridge Aerospace Series, Cambridge
University Press, Cambridge, UK., 2003, pp 76-111.

[21] Redding, D., Adams, N. and Kubiak, E. "Linear-quadratic stationkeeping for the STS orbiter", *Journal of Guidance, Control and Dynamics*, Vol.12, No. 2, 1989, pp 248-255.

[22] Kristiansen, R. and Nicklasson, P. "Spacecraft formation flying: a review and new results on state feedback control", *Acta Astronautica*, Vol 65, No.11, 2009, pp 1537-1552, DOI: 10.1016/j.actaastro.2009.04.014

[23] Lim, H.C., Bang, H., "Adaptive control for satellite formation flying under thrust misalignment", *Acta Astronautica*, Volume 65, Issues 1-2, July-August 2009, pp 112-122. DOI: 10.1016/j.actaastro.2009.01.022.

[24] Mazal, L., Mingotti, G., and Gurfil, P., "Optimal On-Off Cooperative Maneuvers for Long-term Satellite Cluster Flight", *Journal of Guidance, Control, and Dynamics*, Vol. 37, No. 2, March-April 2014, pp. 391-402. DOI: 10.2514/1.61431

[25] Deaconu, G., Louembet, C. and Theron, A. "Minimizing the effects of the navigation uncertainties on the spacecraft rendezvous precision". *Journal of Guidance, Control and Dynamics*, Vol. 37, No. 2 (2014), pp. 695-700. DOI: 10.2514/1.62219

[26] Breger, L., Inalhan, G., Villerson, M. and How, J. "Cooperative spacecraft formation flying: model predictive control with open and closed loop robustness". In P. Gurfil, edito, *Modern Astrodynamics*, chapter 8, Elvesier, New York, 2006, pp 237-277.

[27] Hartley, E., Trodden, P., Richards, A., Maciejowski, J. "Model predictive control system design and implementation for spacecraft rendezvous", *Control Engineering Practice*, Vol. 20, No. 7, 2012, pp 695-713, DOI: 10.1016/j.conengprac.2012.03.009.

[28] How, J. and Villerson, M., "Analysis of the impact of sensor noise on formation flying control," *Proceedings of the 2001 American Control Conference*, Arlington, VA, Vol 5, 2001, pp 3986-3991. DOI:10.1109/ACC.2001.946298

[29] Kall, P. and Wallace, S. "Stochastic programming," Wiley, 1994, pp 110-158

[30] Gavilan, F., Vazquez R. and Camacho, E.F. "Chance-constrained Model Predictive Control for Spacecraft Rendezvous with Disturbance Estimation," *Control Engineering Practice*, Vol.20, No. 2, 2012, pp 111-122. doi: 10.1016/j.conengprac.2011.09.006

[31] Luo, Y., Tang, G., Li, Y., Li, H. "Optimization of Multiple-Impulse, Multiple-Revolution, Rendezvous-Phasing Maneuvers," *Journal of Guidance, Control and Dynamics*, Vol.30, No. 4, 2007, pp 946-952. DOI: 10.2514/1.25620

[32] Ross, I.M. "Space trajectory Optimization and L1-Optimal Control Problems". In P. Gurfil, edito,
*Modern Astrodynamics*, chapter 6, page 155-188. Elvesier Academic press, 2007.

[33] Ben-Tal, A., and Nemirovski, A., "Robust convex optimization," *Mathematics of Operations Research*, Vol 23, No. 4, 1998, pp 769-805. DOI:10.1287/moor.23.4.769

[34] Ben-Tal, A. and Nemirovski A. *Robust solutions of linear programming problems contaminated with uncertain data*. *Math. Programming*, Vol. 88, No. 3, 2000, pp 411-424 DOI: 10.1007/PL00011380

[35] Tschauner, J., "The elliptic orbit rendezvous," *AIAA Journal*, Vol. 5, No. 6, 1967, pp. 1110-1113. doi: 10.2514/3.4145.

[36] Yamanaka, K., and Ankersen, F., "New state transition matrix for relative motion on an arbitrary elliptical orbit," *Journal of Guidance, Control and Dynamics*, Vol. 25, No. 1, 2002, pp 60-66. doi: 10.2514/2.4875

[37] Larsson, R., Berge, S., Bodin, P., and Jönsson, U., "Fuel Efficient Relative Orbit Control Strategies for Formation Flying and Rendezvous within PRISMA" *Advances in the Astronautical Sciences*, 2006, vol. 125, paper 25.

[38] Breger, L. and How, J., "Safe trajectories for autonomous rendezvous of spacecraft," *Journal of Guidance, Control and Dynamics*, Vol. 31, No. 5, 2008, pp 1478-1489. DOI: 10.2514/1.29590

[39] Jaulin, L., Kieffer, M., Didrit, O., and Walter, E., *Applied Interval Analysis with Examples in Parameter and State Estimation, Robust Control and Robotics*. Springer-Verlag, Berlin, Germany 2001, pp 11-44.

[40] Ben-Tal, A., EL Ghaoui, and L., Nemirovski A. *Robust Optimization* Princeton Univ. Press, Princeton, NJ, 2009, pp 3-26