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

HAL
open science

# Metric approximation of minimum time control systems

Jean-Baptiste Caillau, Jean-Baptiste Pomet, Jeremy Rouot

► To cite this version:

Jean-Baptiste Caillau, Jean-Baptiste Pomet, Jeremy Rouot. Metric approximation of minimum time control systems. 2017. hal-01672001

HAL Id: hal-01672001

https://inria.hal.science/hal-01672001

Submitted on 22 Dec 2017

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Metric approximation of minimum time control systems*

J.-B. Caillau†

J.-B. Pomet‡

J. Rouot§

November 2017

## Abstract

Slow-fast affine control systems with one fast angle are considered. An approximation based on standard averaging of the extremal is defined. When the drift of the original system is small enough, this approximation is metric, and minimum time trajectories of the original system converge towards geodesics of a Finsler metric. The asymmetry of the metric accounts for the presence of the drift on the slow part of the original dynamics. The example of the $J_2$ effect in the two-body case in space mechanics is examined. A critical ratio between the $J_2$ drift and the thrust level of the engine is defined in terms of the averaged metric. The qualitative behaviour of the minimum time for the real system is analyzed thanks to this ratio.

**Keywords.** Slow-fast control systems, minimum time, averaging, Finsler metric, $J_2$ potential of two-body problem

**MSC classification.** 49K15, 70Q05

## 1 Averaging of slow-fast minimum time control systems

We consider the following slow-fast control system on an $n$-dimensional manifold $M$:

$$ \dot{I} = \varepsilon F_0(I, \varphi, \varepsilon) + \varepsilon \sum_{i=1}^{m} u_i F_i(I, \varphi, \varepsilon), \quad |u| = \sqrt{u_1^2 + \cdots + u_m^2} \le 1, \quad (1) $$

$$ \dot{\varphi} = \omega(I) + \varepsilon G_0(I, \varphi, \varepsilon) + \varepsilon \sum_{i=1}^{m} u_i G_i(I, \varphi, \varepsilon), \quad \omega(I) > 0, \quad (2) $$

<footnote>*Work supported by the French Space Agency (CNES R&T contract no. R-S13/BS-005-012).</footnote>

<footnote>†LJAD, Univ. Côte d'Azur & CNRS/Inria, Parc Valrose, F-06108 Nice (caillau@unice.fr). Also supported by FMJH and EDF (PGMO grant no. 2016-1753H).</footnote>

<footnote>‡Inria Sophia Antipolis Méditerranée, Univ. Côte d'Azur, 2004 route des lucioles, F-06902 Valbonne (jean-baptiste.pomet@inria.fr).</footnote>

<footnote>§École Polytechnique Féminine, 2 rue Sastre, F-10430 Rosières-près-Troyes (jeremy.rouot@epf.fr).</footnote>
with $I \in M$, $\varphi \in \mathbf{S}^1$, $u \in \mathbf{R}^m$, and fixed extremities $I_0, I_f$, and free phases $\varphi_0, \varphi_f$. All the data is periodic with respect to the single fast angle $\varphi$, and $\omega$ is assumed to be positive on $M$. Extensions are possible to the case of several phases but resonances have then to be taken into account. According to Pontrjagin maximum principle, time minimizing curves are projections onto the base space $M \times \mathbf{S}^1$ of integral curves (extremals) of the maximized Hamiltonian below:

$$H(I, \varphi, p_I, p_\varphi, \varepsilon) := p_\varphi \omega(I) + \varepsilon K(I, \varphi, p_I, p_\varphi, \varepsilon),$$

$$K := H_0 + \sqrt{\sum_{i=1}^{m} H_i^2},$$

$$H_i(I, \varphi, p_I, p_\varphi, \varepsilon) := p_I F_i(I, \varphi, \varepsilon) + p_\varphi G_i(I, \varphi, \varepsilon), \quad i = 0, \dots, m.$$

There are two types of extremals: abnormal ones that live on the level set $\{H=0\}$, and normal ones that evolve on nonzero levels of the Hamiltonian. One defines the averaged Hamiltonian $\bar{K}$ as

$$\bar{K} := \bar{H}_0 + \bar{K}_0, \quad \bar{H}_0 := \langle p_I, \bar{F}_0 \rangle,$$

$$\begin{aligned} \bar{K}_0(I, p_I) &:= \frac{1}{2\pi} \int_0^{2\pi} \sqrt{\sum_{i=1}^{m} H_i^2(I, \varphi, p_I, p_\varphi = 0, \varepsilon = 0)} \, d\varphi \\ &= \frac{1}{2\pi} \int_0^{2\pi} \sqrt{\sum_{i=1}^{m} \langle p_I, F_i(I, \varphi, \varepsilon = 0) \rangle^2} \, d\varphi. \end{aligned}$$

It is smooth on the open set $\Omega := \mathcal{C}\bar{\Sigma}$ where

$$\Sigma := \{(I, p_I, \varphi) \in T^*M \times \mathbf{S}^1 \mid (\forall i = 1, m) : \langle p_I, F_i(I, \varphi, \varepsilon = 0) \rangle = 0\},$$

$$\bar{\Sigma} := \bar{\varpi}(\Sigma) \quad \bar{\varpi} : T^*M \times \mathbf{S}^1 \to T^*M.$$

Indeed, the canonical projection $\bar{\varpi}$ that forgets the fiber $\mathbf{S}^1$ is a closed mapping as the factor $\mathbf{S}^1$ is compact, so $\bar{\Sigma}$ is closed. On also defines the open submanifold $M_0 := \Pi(\Omega)$ of $M$. We assume that $M_0$ is connex.

Under the assumption

$$(A1) \quad \text{rank}\{\partial^j F_i(I, \varphi, \varepsilon = 0) / \partial \varphi^j, \ i = 1, \dots, m, \ j \ge 0\} = n, \ (I, \varphi) \in M \times \mathbf{S}^1,$$

one has

**Proposition 1.** The symmetric part $\bar{K}_0 : (\Omega \subset \mathcal{T}^*M \to \mathbf{R}$ of the tensor $\bar{K}$ is positive definite and 1-homogenous. It so defines a symmetric Finsler co-norm.

*Remark 1.* Condition (A1) is related to the controllability of the original system without drift ($F_0$). It actually amounts to checking the rank of the Lie algebra generated by $F_1, \dots, F_m$ avec and $\partial/\partial\varphi$.

Let us recall that a Finsler norm is a function $F : TM \to \mathbf{R}$ that is smooth on $TM\setminus 0$ and such that

(i) $F(x, \lambda v) = \lambda F(x, v)$, $\lambda > 0$ (the norm is said to be symmetric or absolute value if $F(x, -v) = F(x, v)$),
$$ (ii) \partial^2 F^2(x, v)/\partial v^2 > 0. $$

The fact that the tensor in (ii) depends on $v$ is the main difference with the Riemannian setting. Let now $x$ and $y$ belong to $M$, and let $d(x, y)$ be the infimum of final times $t_f$ over all $\mathcal{C}^1$ curves $\gamma$ connecting the two points with speed bounded by one:

$$ \gamma(0) = x, \quad \gamma(t_f) = y, \\ F(\gamma(t), \dot{\gamma}(t)) \le 1, \quad t \in [0, t_f]. $$

Having so defined the metric $d$ associated with $F$, one defines geodesics to be constant speed curves whose short segments minimize length. Finsler co-norms are defined in the same fashion on the cotangent bundle. Let $F^* : T^*M \to \mathbf{R}$ be smooth on $T^*M \setminus 0$ and such that

$$ (i) F^*(x, \lambda p) = \lambda F^*(x, p), \quad \lambda > 0, $$

$$ (ii) \partial^2(F^*)^2(x, p)/\partial p^2 > 0. $$

Then $F^*$ is a Finsler co-norm, dual to a Finsler norm as both are related through the Legendre transform. More precisely, set

$$ F(x, v) := \max_{p \text{ s.t. } F^*(x,p) \le 1} \langle p, v \rangle. $$

The mapping $F$ defines a Finsler norm whose geodesics are integral curves of the Hamiltonian $F^*$ restricted to the level set $\{F^* = 1\}$ (see, e.g., [12]). One actually has $F^*(x, p) = F(x, v)$, $v := \ell_x^*(p)$, where $\ell_x^* : T_x^*M \to (T_x^*M)^* \simeq T_x M$ is the Legendre transform

$$ \ell_x^* : p \mapsto \frac{1}{2} \frac{\partial^2(F^*)^2}{\partial p^2}(x, p)(p, \cdot) $$

We now assume

$$ (A2) \quad \bar{K}_0(I, \bar{F}_0^*(I)) < 1, \quad I \in M, $$

where $\bar{F}_0^*$ is the inverse Legendre transform of $\bar{F}_0$. Under this new assumption, one has

**Proposition 2.** The tensor $\bar{K} = \bar{H}_0 + \bar{K}_0$ is positive definite and defines an asymmetric Finsler co-norm.

*Remark 2.* For small enough $\varepsilon > 0$, this condition is related to the local controllability of the original system, drift $F_0$ included. It measures the ability of the controlled vector fields $F_1, \dots, F_m$ and their brackets with $\partial/\partial\varphi$ to compensate for the drift. (See next section.)

The geodesics are the integral curves of the Hamiltonian $\bar{K}$ restricted to the level set $\{\bar{K} = 1\}$,

$$ \frac{dI}{d\tau} = \frac{\partial \bar{K}}{\partial p_I}, \quad \frac{dp_I}{d\tau} = -\frac{\partial \bar{K}}{\partial I}, \\ I(0) = I_0, \quad I(\tau_f) = I_f, \quad \bar{K}(I_0, p_I(0)) = 1, $$
and $\tau_f = d(I_0, I_f)$ for minimizing ones. Up to a reparameterization of time ($ds = \omega(I)d\tau$), the geodesics also are integral curves of

$$ \bar{h}(I, p_I) := \frac{\bar{K}(I, p_I) - 1}{\omega(I)} $$

restricted to the level $\{\bar{h} = 0\}$. We study in the next section the convergence properties of the original system towards this metric when $\varepsilon \to 0$.

## 2 Approximation properties

On order to identify the slow and fast part on the extremal flow of the original system, we use the following ansatz. For $\varepsilon > 0$, we normalize min. time extremals according to $\bar{K}(I_0, p_I(0)) = 1$. (Under the previous assumptions, $\bar{K}$ indeed defines a Minkowski norm on the fiber $T_{I_0}^*M_0$.) Now, as

$$ H = H(0) = p_{\varphi}(0)\omega(I_0) + \varepsilon K(I_0, \varphi(0), p_I(0), p_{\varphi}(0), \varepsilon) = O(\varepsilon) $$

since $p_{\varphi}(0) = 0$, and since $\varphi(0) \in S^1$ and $p_I(0)$ are bounded. So, outside resonance ($\omega$ is assumed to be positive on $M_0$),

$$ p_{\varphi} = -\varepsilon \cdot \frac{K - H/\varepsilon}{\omega(I)} = O(\varepsilon). $$

This ensures the possibility of division par $\varepsilon$, that is the existence of a smooth function $h$ such that $p_{\varphi} = -\varepsilon h$. This is crucial since

$$ \dot{p}_I = -p_{\varphi}\omega'(I) + \varepsilon \frac{\partial K}{\partial I}(I, \varphi, p_I, p_{\varphi}, \varepsilon) $$

this allows to identify $p_I$ as a slow variable ($\dot{p}_I = O(\varepsilon)$). More precisely, the following holds.

**Lemma 1.** For small enough $\varepsilon > 0$, the fixed point equation

$$ h = \frac{K(I, \varphi, p_I, -\varepsilon h, \varepsilon) - k}{\omega(I)} $$

has a unique solution $h = h(I, \varphi, p_I, k, \varepsilon)$ smoothly depending on $(I, \varphi, p_I, k, \varepsilon)$.

This symplectic reduction eliminates $p_{\varphi}$ and one can rewrite the extremal flow in the standard form to perform averaging:

$$ \dot{I} = \varepsilon \frac{\partial H}{\partial p_{\varphi}} \frac{\partial h}{\partial p_I}, \quad \dot{p}_I = -\varepsilon \frac{\partial H}{\partial p_{\varphi}} \frac{\partial h}{\partial I} = O(\varepsilon), $$

$$ \dot{\varphi} = \omega(I) + O(\varepsilon). $$

Changing time $t$ to $s = \varepsilon\varphi$,

$$ \frac{dI}{ds} = \frac{\partial h}{\partial p_I}(I, s/\varepsilon + \varphi(0), p_I, k(\varepsilon), \varepsilon), $$

$$ \frac{dp_I}{ds} = -\frac{\partial h}{\partial I}(I, s/\varepsilon + \varphi(0), p_I, k(\varepsilon), \varepsilon), \quad s \in [0, s_f], $$
with

$$k(\varepsilon) := K(I_0, \varphi(0), p_I(0), p_{\varphi} = 0, \varepsilon) = H/\varepsilon \ge 0.$$

Using averaging, we are led to approximate the curves of the original system with the integral curves of the previously defined $\bar{h}$,

$$\begin{aligned} \bar{h}(I, p_I) &= \frac{1}{2\pi} \int_{0}^{2\pi} h(I, \varphi, p_I, k = 1, \varepsilon = 0) \, d\varphi, \\ &= \frac{\bar{K}(I, p_I) - 1}{\omega(I)}. \end{aligned}$$

At this stage, it is not clear whether the choice $k = 1$ is justified or not, and whether we must restrict to $\{\bar{h} = 0\}$ or not.

In time $s = \varepsilon\varphi$, we define the shooting function

$$\bar{S}(s_f, p_{I_0}) := (I(s_f, I_0, p_{I_0}) - I_f, \bar{h}(I_0, p_{I_0}))$$

associated with the two-point boundary value problem (denoting $z = (I, p_I)$)

$$\frac{dz}{ds}(s) = \vec{h}(z(s)), \quad s \in [0, s_f],$$

$$I(0) = I_0, \quad I(s_f) = I_f, \quad \bar{h}(I(0), p_I(0)) = 0.$$

For $\varepsilon > 0$, and for any fixed $\varphi_0 \in \mathbf{S}^1$, we use the same normalization of $p_I(0)$ by $\bar{h} = 0$ and also define the shooting function

$$S_{\varepsilon}(s_f, p_{I_0}) := (I(s_f, I_0, p_{I_0}, \varepsilon) - I_f, \bar{h}(I_0, p_{I_0})) \quad (3)$$

associated with the two-point boundary value problem

$$\frac{dI}{ds} = \frac{\partial h}{\partial p_I}(I, s/\varepsilon + \varphi_0, p_I, k = 1, \varepsilon),$$

$$\frac{dp_I}{ds} = -\frac{\partial h}{\partial I}(I, s/\varepsilon + \varphi_0, p_I, k = 1, \varepsilon), \quad s \in [0, s_f].$$

The approximation result below is key for the rest of the study.

**Proposition 3.** Let $I_0$ and $I_f$ in $M_0$, and let $(\bar{s}_f, \bar{p}_{I_0})$ be a regular zero of $\bar{S}$. For any $\varepsilon > 0$, and whatever $\varphi_0 \in \mathbf{S}^1$, there exists a zero $(s_f(\varepsilon), p_{I_0}(\varepsilon))$ of $S_{\varepsilon}$ such that

$$s_f(\varepsilon) \rightarrow \bar{s}_f, \quad p_{I_0}(\varepsilon) \rightarrow \bar{p}_{I_0} \quad \text{quand } \varepsilon \rightarrow 0.$$

The proof of this proposition relies on the following fixed point result applied to the family of shooting functions (3) for $\varepsilon > 0$.

**Lemma 2.** Let $f : \mathbf{R}^n \rightarrow \mathbf{R}^n$ be a continuously differentiable mapping having a regular zero at $x = 0$, and let $f_{\varepsilon} : \mathbf{R}^n \rightarrow \mathbf{R}^n$, $\varepsilon > 0$, be continuous mappings converging uniformly towards $f$ on a neighbourhood of the origin when $\varepsilon \rightarrow 0$. Then, there exist $\varepsilon_0 > 0$ together with a mapping $x : [0, \varepsilon_0] \rightarrow \mathbf{R}^n$ that is continuous at $\varepsilon = 0$ such that $x(0) = 0$ and

$$f_{\varepsilon}(x(\varepsilon)) = 0, \quad \varepsilon \in (0, \varepsilon_0].$$
<figure><img src="image_3.png" /><figcaption>Figure 1: Typical swallowtail singularities of the value function $\varepsilon \mapsto t_f(\varepsilon)$ obtained when following the characteristics with a continuation on $\varepsilon$. Numerical simulation from [10].</figcaption></figure>

**Remark 3.** One cannot expect $\mathcal{C}^1$ regularity for the value function $\varepsilon \mapsto t_f(\varepsilon)$. Indeed, if one obtains the value function using a continuation on $\varepsilon$ and computing min. time extremals, the continuation allows to follow the characteristics and to cross the singularities of the associated Hamilton-Jacobi-Bellman equation. Due to the existence of local minima (because of the free initial and final phases), swallowtail singularities are encountered, typically [10]. These singularities accumulate as $\varepsilon \to 0$. (See Figure 1.)

**Lemma 3 (Verification lemma).** Let $\varepsilon > 0$, and let $\varphi_0 \in \mathbf{S}^1$. To any zero of $S_\varepsilon$ correspond a extremal $(I, \varphi, p_I, p_\varphi)$ of H and a final time $t_f$ such that

$$I(0) = I_0, \quad I(t_f) = I_f, \quad p_\varphi(0) = O(\varepsilon), \quad p_\varphi(t_f) = O(\varepsilon).$$

In order to state our main convergence result, we make the following strong assumptions on the metric defined by $\bar{K}$.

(A3) The metric is geodesically convex on $M_0$.

(A4) Whatever $I_0$ and $I_f$ in $M_0$, $I_f \notin \text{Cut}(I_0)$, there exist $\varepsilon_0 > 0$, a compact neighbourhood $K(I_0, I_f)$ of the minimizing geodesic, and $\eta > 0$ such that, for any $\varepsilon \in (0, \varepsilon_0]$, any admissible trajectory (whose final time is $t_f$) not contained in $K(I_0, I_f)$ it holds that

$$\varepsilon t_f \ge d(I_0, I_f) + \eta.$$
A partial result in the direction of (A3) is proved in [6] for the Finsler metric associated with the two-body potential (geodesic convexity of the so-called meridian half-planes in the of 2D case; compare with [5]). A quantitative study of the original dynamical system is required for the estimation in (A4). In the two-body case, one has to analyze the effect of the singularities of the dynamics at $n = 0$ (parabolic resonance) and at $n = \infty$ (pure collision) ($n$ being the mean motion). The issue of loss of regularity due to $\pi$-singularities must also be addressed [4]. Note that, in the statement of (A4), we use the fact that provided the final point does not belong to the cut locus of the initial one, a unique (forward) geodesic of the Finsler metric connects them.

**Proposition 4.** Let $I_0$ and $I_f$ belong to $M_0$, $I_f \notin \text{Cut}(I_0)$. Then, for small enough $\varepsilon > 0$, existence holds for the original minimum time control problem (1)-(2).

**Theorem 1.** Let $I_0$ and $I_f$ belong to $M_0$, $I_f \notin \text{Cut}(I_0)$. Let $(I_\varepsilon, \varphi_\varepsilon, p_{I_\varepsilon}, p_{\varphi_\varepsilon})_\varepsilon$ be a family of minimizing extremals, and let $(\bar{t}_f(\varepsilon))_\varepsilon$ be the associated family of minimum times. Then, denoting $z_\varepsilon := (I_\varepsilon, p_{I_\varepsilon})$, one has

$$ \|z_{\varepsilon} - \bar{z}\|_{\infty} = O(\varepsilon) + O(k(\varepsilon) - 1), \quad \varepsilon \bar{t}_f(\varepsilon) \to d(I_0, I_f), \quad \varepsilon \to 0, $$

where $\bar{z}$ is the Hamiltonian lift of the minimizing geodesic connecting $I_0$ to $I_f$.

### 3 Application to space mechanics

We consider the the two-body potential case,

$$ \ddot{q} = -\mu \frac{q}{|q|^3} + \frac{u}{M}, \quad |u| \le T_{\max}. $$

Thanks to the super-integrability of the $-1/|q|$ potential, the minimum time control system is slow-fast with only angle (the longitude of the evolving body) if one restricts to the case of transfers between elliptic orbits ($\mu$ is the gravitational constant). In the non-coplanar situation, we have to analyze a dimension five symmetric Finsler metric. In order to account for the Earth non-oblateness, we add to the dynamics a small drift $F_0$ on the slow variables. In the standard equinoctial orbit elements, $I = (a, e, \omega, \Omega, i)$, the $J_2$ term of order $1/|q|^3$ of the Earth potential derives from the additional potential ($r_e$ being the equatorial radius)

$$ R_0 = \frac{\mu J_2 r_e^2 (1 - e^2)^{-3/2}}{|q|^3} \left( \frac{1}{2} - \frac{3}{4} \sin^2 i + \frac{3}{4} \sin^2 i \cos(2\omega + 2 - \varphi) \right). $$

As a result, the system now has to small parameters (depending on the initial condition). One is due to the $J_2$ effect, the other to the control:

$$ \varepsilon_0 = \frac{3J_2r_e^2}{2a_0^2}, \quad \varepsilon_1 = \frac{a_0^2 T_{\max}}{\mu M}. $$

Here, $a_0$ is the initial semi-major axis, $T_{\max}$ the maximum level of thrust, and $M$ the spacecraft mass. We make a reduction to a single small parameter as
follows: Defining $\varepsilon := \varepsilon_0 + \varepsilon_1$ and $\lambda := \varepsilon_0/(\varepsilon_0 + \varepsilon_1)$, one has

$$ 
\begin{aligned}
\dot{I} &= \varepsilon_0 F_0(I, \varphi) + \varepsilon_1 \sum_{i=1}^{m} u_i F_i(I, \varphi), \\
&= \varepsilon \left( \lambda F_0(I, \varphi) + (1 - \lambda) \sum_{i=1}^{m} u_i F_i(I, \varphi) \right).
\end{aligned}
$$

There are two regimes depending on whether the $J_2$ effect is small against the control ($\varepsilon_0 \ll \varepsilon_1$ and $\lambda \to 0$) or not ($\varepsilon_0 \gg \varepsilon_1$ and $\lambda \to 1$). The critical ratio on $\lambda$ can be explicitly computed in metric terms.

**Proposition 5.** In the average system of the two-body potential including the $J_2$ effect, $\bar{K} = \lambda \bar{H}_0 + (1 - \lambda) \bar{K}_0$ is a metric tensor if and only if $\lambda < \lambda_c(I)$ with

$$ \lambda_c(I) = \frac{1}{1 + \bar{K}_0(I, \bar{F}_0^*(I))}. $$

The relevance of this critical ratio for the qualitative analysis of the original system is illustrated by the numerical simulations displayed in Figures 2 to 6. For a given initial condition $I_0$ on the slow variables, we let the drift $F_0$ alone act: We integrate the flow of $F_0$ during a short positive duration $\tau_d$, then compute the trajectory of the averaged system to go from this point $I(\tau_d)$ back to $I_0$. For $\lambda < \lambda_c(I_0)$, the tensor $\bar{K}$ is a metric one, and this trajectory is a geodesic. As $\tau_d$ tends to zero, the time $\tau_f$ to come back from $I(\tau_d)$ tends to zero when $\lambda < \lambda_c(I_0)$. For $\lambda \ge \lambda_c(I_0)$, finitess of this time indicates that global properties of the system still allows to control it although the metric character of the approximation does not hold anymore. (See Figure 2.) The behaviour of $\tau_f$ measures the loss in performance as $\lambda$ approaches the critical ratio. This critical value depends on the initial condition and gives an asymptotic estimate of whether the thrust dominates the $J_2$ effect or not. Beyond the critical value, the system is still controllable, but there is a drastic change in performance. As the original system is approximated by the average one, this behaviour is very precisely reproduced on the value function of the original system for small enough $\varepsilon$. (See Figures 3 to 6.)

## References

[1] Agrachev, A. A.; Sachkov, Y. L. Control Theory from the Geometric Viewpoint. Springer, 2004.

[2] Beletsky, V. V. Essays on the motion of celestial bodies. Birkhäuser, 1999.

[3] BepiColombo mission: sci.esa.int/bepicolombo

[4] Bombrun, A.; Pomet, J.-B. The averaged control system of fast oscillating control systems. SIAM J. Control Optim. 51 (2013), no. 3, 2280-2305.

[5] Bonnard, B.; Caillau, J.-B. Geodesic flow of the averaged controlled Kepler equation. Forum Math. 21 (2009), no. 5, 797-814.
<figure><img src="image_4.png" /><figcaption>Figure 2: Value function $\lambda \mapsto \tau_f(\lambda)$, $\tau_d \to 0$ (averaged system). On this example, $a = 30$ Mm, $e = 0.5$, $\omega = \Omega = 0$, $i = 51$ degrees (strong inclination), and $\lambda_c \approx 0.4239$. The value function is portrayed for $\tau_d = 1e-2, 1e-3, 1e-4, 1e-5$.</figcaption></figure>

[6] Bonnard, B.; Henninger, H.; Nemcova, J.; Pomet, J.-B. Time versus energy in the averaged optimal coplanar Kepler transfer towards circular orbits. *Acta Appl. Math.* **135** (2015), no. 2, 47–80.

[7] Caillau, J.-B.; Cots, O.; Gergaud, J. Differential pathfollowing for regular optimal control problems. *Optim. Methods Softw.* **27** (2012), no. 2, 177–196.

[8] Caillau, J.-B.; Daoud, B. Minimum time control of the restricted three-body problem. *SIAM J. Control Optim.* **50** (2012), no. 6, 3178–3202.

[9] Caillau, J.-B.; Dargent, T.; Nicolau, F. Approximation by filtering in optimal control and applications. *IFAC PapersOnLine* **50** (2017), no. 1, 1649–1654. Proceedings of the 20th IFAC world congress, Toulouse, July 2017.

[10] Caillau, J.-B.; Farrés, A. On local optima in minimum time control of the restricted three-body problem. *Recent Advances in Celestial and Space Mechanics*, 209–302, Mathematics for industry **23**, Springer, 2016.

[11] Caillau, J.-B.; Noailles, J. Coplanar control of a satellite around the Earth. *ESAIM Control Optim. and Calc. Var.* **6** (2001), 239–258.

[12] A brief introduction to Finsler geometry. Course notes, 2006. math.aalto.fi/~fdahl/finsler/finsler.pdf

[13] Dargent, T.; Nicolau, F.; Pomet, J.-B. Periodic averaging with a second order integral error. *IFAC PapersOnLine* **50** (2017), no. 1, 2892–2897. Proceedings of the 20th IFAC world congress, Toulouse, July 2017.
<figure><img src="image_5.png" /><figcaption>Figure 3: Value function $\lambda \mapsto \tau_f(\lambda)$, $\tau_d \to 0$ (original system, $\varepsilon = 1e-3$). On this example, $a = 30$ Mm, $e = 0.5$, $\omega = \Omega = 0$, $i = 51$ degrees (strong inclination), and $\lambda_c \approx 0.4239$. The behaviour of the value function for the original system matches very precisely the behaviour of the averaged one. (See also Figure 4 for a even lower value of $\varepsilon$.)</figcaption></figure>

[14] Edelbaum, T. N. Optimal low-thrust rendez-vous and station keeping, *AIAA J. 2* (1964), no. 7, pp. 1196–1201.

[15] Geffroy, S.; Epenoy, R. Optimal low-thrust transfers with constraints. Generalization of averaging techniques. *Acta Astronaut. 41* (1997), no. 3, 133–149.

[16] Hampath software: *hampath.org*

[17] Lisa Pathfinder mission: *sci.esa.int/lisa-pathfinder*
<figure><img src="image_6.png" /><figcaption>Figure 4: Value function $\lambda \mapsto \tau_f(\lambda)$, $\tau_d \to 0$ (original system, $\varepsilon = 1e-4$). On this example, $a = 30$ Mm, $e = 0.5$, $\omega = \Omega = 0$, $i = 51$ degrees (strong inclination), and $\lambda_c \approx 0.4239$. The behaviour of the value function for the original system matches very precisely the behaviour of the averaged one.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5: Value function $\lambda \mapsto \tau_f(\lambda)$, $\tau_d \to 0$ (original system, $\varepsilon = 1e-3$). On this example, $a = 11.675$ Mm, $e = 0.75$, $\omega = \Omega = 0$, $i = 7$ degrees (weak inclination), and $\lambda_c \approx 0.2287$. The behaviour of the value function for the original system matches very precisely the behaviour of the averaged one. (See also Figure 6 for a even lower value of $\varepsilon$.)</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6: Value function $\lambda \mapsto \tau_f(\lambda)$, $\tau_d \to 0$ (original system, $\varepsilon = 1e-3$). On this example, $a = 11.675$ Mm, $e = 0.75$, $\omega = \Omega = 0$, $i = 7$ degrees (weak inclination), and $\lambda_c \approx 0.2287$. The behaviour of the value function for the original system matches very precisely the behaviour of the averaged one.</figcaption></figure>