# On the convergence of time-optimal maneuvers of fast-oscillating control systems

Lamberto Dell'Elce, Jean-Baptiste Caillau, Jean-Baptiste Pomet

► To cite this version:

Lamberto Dell'Elce, Jean-Baptiste Caillau, Jean-Baptiste Pomet. On the convergence of time-optimal maneuvers of fast-oscillating control systems. ECC 2021 - European Control Conference, Jun 2021, Virtual, France. pp.2008-2013, 10.23919/ECC54610.2021.9655101 . hal-03391620

HAL Id: hal-03391620

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

Submitted on 21 Oct 2021

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.
# On the convergence of time-optimal maneuvers of fast-oscillating control systems

Lamberto Dell'Elce,¹ Jean-Baptiste Caillau² and Jean-Baptiste Pomet¹

**Abstract**— For a control system with one fast periodic variable, with a small parameter measuring the ratio between time derivatives of fast and slow variables, we consider the Hamiltonian equation resulting from applying Pontryagin maximum principle for the minimum time problem with fixed initial and final slow variables and free fast variable. One may perform averaging at least under normalization of the adjoint vectors and define a “limit” average system. The paper is devoted to the convergence properties of this problem as the small parameter tends to 0. We show that using the right transformations between boundary conditions of the "real" and average systems leads to a reconstruction of the fast variable on interval of times of order $1/\epsilon$ where $\epsilon$ is the small parameter. This is only evidenced numerically in this paper. Relying on this, we propose a procedure to efficiently reconstruct the solution of the two point boundary problem for nonzero $\epsilon$ using only the solution of the average optimal control problem.

## I. INTRODUCTION

Periodic averaging is extremely common. For a differential equation in a certain number of slow variables and one periodic fast variable, with a small parameter $\epsilon > 0$ measuring the ratio, it states that the distance between solutions of the original and average systems is of order $\epsilon$ on a time-interval of order $1/\epsilon$. See for instance [1], [2].

Here, motivated by low-thrust orbital transfer in space mechanics (see also [3], [4], [5], [6], [7]), we consider minimum-time optimal control for a control system (1), and use averaging on the Hamiltonian equation resulting from Pontryagin Maximum Principle. Averaging for optimal control was also studied in [8].

As already noted in [9], one has to be careful in sorting fast and slow variables among the adjoint variables in this equation, and averaging only holds with properly normalized adjoint vectors. Then, convergence is obviously in a weak sense because, contrary to Cauchy problems, two point boundary value problems do not have unique solutions (and in fact, in the problem with fixed slow variables and free fast variables, more and more solutions appear as $\epsilon$ goes to zero).

We study this convergence and give a procedure to recover the solutions of the optimal control problem from solving the average equations only. We implement transformations between initial and final conditions of the slow variables of the average and original system; these are the same as in second order averaging mentioned above (in the context of optimal control they were already discussed in [5]). We include a reconstruction of the fast variable in the average system.

Although there is no proof per se, we do witness very clearly, in our numerical experiments, a reconstruction of the fast variable, with an error that is possibly not of order 1

with respect to $\epsilon$, but that tends to zero. Although a second order term in the average dynamics would be needed to achieve this in the classical Cauchy problem, simulations are extremely clear. We then describe an original procedure to recover the initial conditions of the adjoint vectors as well as of the fast variables from solving the average system. This procedure relies on the assumption that reconstruction of the fast variables is effective. One of its interest is that it explains very well the fact that there are more and more solutions as $\epsilon \to 0$. Another and more practical interest is that it provides very accurate initial guesses for the original two point boundary value problem.

## II. AVERAGING FAST-OSCILLATING OPTIMAL CONTROL SYSTEMS

The time-optimal problem of slow-fast dynamical systems, which is referred to as “original problem” hereafter, is introduced first. Then, its average counterpart is defined.

### A. Original problem

We consider the following slow-fast controlled dynamics where $\epsilon$ is to be understood as a small parameter, scaling the slow part:

$$ 
\begin{aligned}
\dot{I} &= \epsilon \left[ F_0(I, \varphi, \epsilon) + \sum_{i=1}^{m} u_i F_i(I, \varphi, \epsilon) \right], \\
\dot{\varphi} &= \omega(I) + \epsilon \left[ G_0(I, \varphi, \epsilon) + \sum_{i=1}^{m} u_i G_i(I, \varphi, \epsilon) \right], \\
|u| &= \sqrt{u_1^2 + \cdots + u_m^2} \le 1,
\end{aligned}
\quad (1) $$

The slow state variables, $I$, live on some manifold $M$ of dimension $n$, and there are fixed endpoint conditions $I_0, I_f$. Conversely, $\varphi$ is the fast variable, and all the data (that is the vector fields $F_i$ and $G_i$) are $2\pi$-periodic with respect to it. There are no endpoint conditions neither on $\varphi \in \mathbf{R}$, nor on its angle class in $\mathbf{S}^1 \simeq \mathbf{R}/2\pi\mathbf{Z}$. We assume that the pulsation $\omega(I)$ is uniformly bounded on $M$ by some positive constant. The control $u$ is prescribed to the unit Euclidean ball of $\mathbf{R}^m$. We are interested in the minimum time control of this system. Pontryagin maximum principle implies that time minimizing trajectories are projections on $M$ of integral curves of the Hamiltonian system defined on $T^*(M \to \mathbf{R})$ by

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

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

and

$$H_i(I, \varphi, p_I, p_\varphi, \varepsilon) := p_I F_i(I, \varphi, \varepsilon) + p_\varphi G_i(I, \varphi, \varepsilon)$$

for $i = 0, \dots, m$. Namely,

$$\begin{aligned}
\dot{I} &= \varepsilon \frac{\partial K}{\partial p_I}, & \dot{p}_I &= -p_\varphi \frac{\partial \omega}{\partial I} - \varepsilon \frac{\partial K}{\partial I}, \\
\dot{\varphi} &= \omega(I) + \varepsilon \frac{\partial K}{\partial p_\varphi}, & \dot{p}_\varphi &= -\varepsilon \frac{\partial K}{\partial \varphi}.
\end{aligned} \quad (2)$$

The initial and final phases $\varphi(0)$ and $\varphi(t_f)$ are free, so that we look for integral curves of $H$ such that $p_\varphi(0) = p_\varphi(t_f) = 0$ (transversality condition of Pontryagin maximum principle). Hence, the solution of the two-point boundary problem is achieved by finding a zero $(t_f, p_{I0}, \varphi_0)$ of the shooting function

$$S_{\varepsilon}(t_f, p_{I0}, \varphi_0, I_0, I_f) = \begin{bmatrix} I(t_f, I_0, \varphi_0, p_{I0}, 0, \varepsilon) - I_f \\ p_\varphi(t_f, I_0, \varphi_0, p_{I0}, 0, \varepsilon) \\ |p_{I0}| - 1 \end{bmatrix} \quad (3)$$

parameterized by $I_0$, $I_f$ and $\varepsilon$. Here, $I(t, I_0, \varphi_0, p_{I0}, 0, \varepsilon)$ and $p_\varphi(t, I_0, \varphi_0, p_{I0}, 0, \varepsilon)$ denote the unique solution at time $t$ of the Hamiltonian system, Eq. (2), with initial conditions $I_0, \varphi_0, p_{I0}, p_\varphi(0) = 0$. The constraint $|p_{I0}| - 1$ is associated to the unknown $t_f$, i.e., the time needed to reach $I_f$ and $p_\varphi(t_f) = 0$ starting from the aforementioned initial conditions. Such normalization is always possible by homogeneity. The compact notation $S_\varepsilon(t_f, p_{I0}, \varphi_0)$ is used in the reminder whenever boundary conditions $I_0$ and $I_f$ are understood.

## B. Average problem

The term $-p_\varphi \frac{\partial \omega}{\partial I}$ in the dynamics of the adjoint state $p_I$ is not a priori small, but rather of order zero with respect to $\varepsilon$. Its presence therefore prevents us from using standard averaging techniques to approximate this Hamiltonian system when $\varepsilon$ goes to zero. Instead, let us introduce the re-scaling

$$s := \varepsilon t, \quad \psi := \varepsilon \varphi,$$

and define

$$\bar{H}(J, \psi, p_J, p_\psi) := p_\psi \omega(J) + \bar{K}(J, p_J)$$

where $\bar{K} = \bar{H}_0 + \bar{K}_0$ with

$$\bar{K}_0(J, p_J) := \frac{1}{2\pi} \int_0^{2\pi} \sqrt{\sum_{i=1}^{m} H_i^2(J, \varphi, p_J, p_\varphi = 0, \varepsilon = 0)} d\varphi,$$

and

$$\bar{H}_0(J, p_J) := \left\langle p_J, \frac{1}{2\pi} \int_0^{2\pi} F_0(J, \varphi, \varepsilon = 0) d\varphi \right\rangle.$$

The "average" Hamiltonian $\bar{H}$ is smooth on $\Omega$, the open complement of

$$\bar{\Sigma} := \varpi(\Sigma), \quad \varpi : T^*M \times \mathbf{R} \to T^*M,$$

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

So is the original Hamiltonian $H$, for small enough $\varepsilon$, on the complementary of a thickening of $\Sigma$.

The variable $\psi$ does not appear in $\bar{H}$, so for any value of the first integral $p_\psi$, one can consider the two point boundary value problem in the rescaled time $s$ (we denote $s' = d/ds$),

$$J' = \frac{\partial \bar{K}}{\partial p_J}, \quad p'_J = -p_\psi \frac{\partial \omega}{\partial J} - \frac{\partial \bar{K}}{\partial J},$$

for fixed endpoints $J_0, J_f$ in $M$. Solving this problem is equivalent of finding a zero $(s_f, p_{J0})$ of the average shooting function

$$S(s_f, p_{J0}, J_0, J_f, p_\psi) = \begin{bmatrix} J(s_f, J_0, p_{J0}, p_\psi) - J_f \\ |p_{J0}| - 1 \end{bmatrix}$$

parameterized by $J_0, J_f$ and $p_\psi$. In analogy with the notation used in Eq. (3), $J(s, J_0, p_{J0}, p_\psi)$ denotes the unique solution (assumed to be defined up to this point) at time $s$ of the average Hamiltonian system with initial conditions $J_0, p_{J0}$ and $p_\psi$. Given a solution $(s_f, p_{J0})$, a simple integration allows to compute

$$\Delta\psi(J_0, J_f, p_\psi) = \psi(s_f) - \psi(0)$$

because $\psi$ is cyclic.

Note that, in general, the shooting problem has several solutions that need to be compared (e.g., by using additional optimality conditions), so that we should explicit the dependence of $\Delta\psi$ with respect to the particular solution retained. For the sake of simplicity, we assume that for a given $p_\psi$, there is only one solution, hence the dependence on $p_\psi$ of $(s_f, p_{J0})$, and so of $\Delta\psi$.

## III. CONVERGENCE OF THE TRAJECTORIES

After providing a practical motivation of this study emphasizing difficulties in the analysis of the convergence of a given solution of Problem (3) for $\varepsilon \to 0$, we discuss how the small parameter, $\varepsilon$, the average adjoint to the fast variable, $p_\psi$, and the boundary phases, $\varphi(0)$ and $\varphi(t_f)$, relate each others. Then, we outline an algorithm that is capable of approximating all the different solutions of the original problem for a given $\varepsilon$ by leveraging on solutions of its average counterpart only. This approximation can serve as a reliable initial guess to solve the original problem.

### A. Motivation of the study

The numerical example discussed in Section IV is now used to disclose the fundamental questions that motivated our study (we are only interested to qualitative behaviors in this paragraph; all details of the case study are provided in Section IV):

Let $(t_f^*, p_{I0}^*, \varphi_0^*)$ be a solution of the original problem, namely $S_\varepsilon(t_f^*, p_{I0}^*, \varphi_0^*) = 0$. Continue $(t_f^*, p_{I0}^*, \varphi_0^*)$ for decreasing $\varepsilon$. Does a solution of the average problem, $(s_f^*, p_{J0}^*, p_\psi^*)$, exist such that $t_f^*/\varepsilon$ and $p_{I0}^*$ converge to $s^*$ and $p_{J0}^*$, respectively?
<figure><img src="image_1.png" /><figcaption>Fig. 1. Differential continuation of a solution of the original system for decreasing ε (colored curve). The black-dashed lines denote solutions of the average problem.</figcaption></figure>

First, we note that the realization of this experiment entails numerical difficulties because $t_f^* = O(\epsilon^{-1})$ and small $\epsilon$ requires a stringent tolerance of the integrator. Hence, it is not possible to approach $\epsilon = 0$ too closely. However, Fig. 1 suggests that the convergence to a single solution of the average problem does not occur, and that the continued branch seems to approach to a subset of solutions of the average problem (not a single one) parameterized by $p_\psi$. This is because $\phi_0^*$, circulates in $\mathbb{S}^1$. This observation lends itself to the following question.

What $p_\psi^*$ should be used, and how does it relate to $\epsilon$ and $\phi_0^*$?

Inspection of Fig. 1 reveals that average solutions can be parameterized by $p_\psi^*$. In addition, the closest average solution to $(t_f^*/\epsilon, p_{I0}^*)$, i.e., the one that should serve as initial guess to the original problem, depends on the value of $\epsilon$. Modeling the relationship between $\epsilon, p_\psi^*$, and $\phi_0^*$ is the key to predict the evolution of the continued trajectory of Fig. 1.

Should we consider the same boundary conditions for both the original and the average problems?

The relation between $\epsilon, p_\psi^*$, and $\phi_0$ (discussed in the reminder of this section) requires the refinement of the boundary conditions of the average problem. In fact, using $(J_0, J_f) = (I_0, I_f)$ does not yield an accurate enough approximation of the original trajectory, as depicted in Fig. 2. Here, the change of the initial phase is sufficient to obtain remarkably different trajectories of the original problem. The near-identity transformation discussed in Section III-B is mandatory to generate average trajectories that remain sufficiently close to the original ones.

## B. Considerations

Using the same endpoint conditions for the original and average problems, i.e., $(J_0, J_f) = (I_0, I_f)$, does not provide with satisfactory matching between original and average

<figure><img src="image_2.png" /><figcaption>Fig. 2. Trajectory of the orbital eccentricity for $\epsilon = 10\epsilon_0$. Solid lines are solutions of the original problem. The black-dashed line is the solution of the average TPBVP for $p_\psi = 0$. Initial conditions are the same for all trajectories (only $\phi(0)$ is different).</figcaption></figure>

trajectories as emphasized above. A crucial step in our approach is to refine the approximation by making an $O(\epsilon)$ correction on these endpoints. To this end, let us define near identity transformations

$$ 
\begin{aligned}
\hat{I} &= J + \epsilon\nu_I(J, p_J, \varphi), \\
\hat{p}_I &= p_J + \epsilon\nu_{p_I}(J, p_J, \varphi), \\
\hat{p}_\varphi &= \epsilon p_\psi + \epsilon\nu_{p_\varphi}(J, p_J, \varphi),
\end{aligned}
 $$

using the usual ansatz: for each variable $X$ ($X$ equal to $I, p_I$ or $p_\varphi$) of the original system, $\nu_X$ is searched so that the reconstruction $\tilde{X}$ from the average system is an approximation of order two in $\epsilon$ of its dynamics:

$$ \dot{\tilde{X}}(t) = \dot{X}(t) + O(\epsilon^2). $$

This leads to the following choices, first for $\nu_{p_\varphi}$:

$$ \frac{\partial \nu_{p_\varphi}}{\partial \varphi} = - \frac{1}{\omega(J)} \frac{\partial K}{\partial \varphi} (J, \varphi, p_J, p_\varphi = 0, \epsilon = 0). $$

Note that $K$ is evaluated at $p_\varphi = 0$, which is compatible with the $O(\epsilon^2)$ discrepancy. Then $\nu_{p_\varphi}$ is chosen of zero average. Similarly,

$$ \frac{\partial \nu_{p_I}}{\partial \varphi} = \frac{1}{\omega(J)} \left[ - \nu_{p_\varphi}(J, p_J, \varphi) \frac{\partial \omega}{\partial J} - \frac{\partial K}{\partial I}(J, \varphi, p_J, p_\varphi = 0, \epsilon = 0) + \frac{\partial \bar{K}}{\partial J}(J, p_J) \right]. $$

The crucial point in the previous expression is the addition of the term

$$ -\nu_{p_\varphi}(J, p_J, \varphi) \frac{\partial \omega}{\partial J} $$

that cannot be neglected. Finally,

$$ \frac{\partial \nu_I}{\partial \varphi} = \frac{1}{\omega(J)} \left( \frac{\partial K}{\partial p_I}(J, \varphi, p_J, p_\varphi = 0, \epsilon = 0) - \frac{\partial \bar{K}}{\partial J}(J, p_J) \right). $$

Again, $\nu_{p_I}$ and $\nu_I$ can be chosen of zero average. We now use these transformations to construct an approximate solution of the original problem for a given value of $\epsilon$.
Transversality conditions ($p_\varphi(0) = p_\varphi(t_f) = 0$) are used to determine possible values for the endpoint phases by solving $\hat{p}_\varphi(0) = \hat{p}_\varphi(t_f) = 0$ that is, up to some order in $\varepsilon$,

$$p_{\psi} + \nu_{p_{\varphi}}(J_0, p_{J0}, \varphi_0) = p_{\psi} + \nu_{p_{\varphi}}(J_f, p_J(s_f), \varphi_f) = 0. \quad (4)$$

These two equations usually have several solutions, all parameterized by $p_\psi$. Let $(\varphi_0, \varphi_f)$ be one of them. We eventually want to improve our current estimation which is based on a solution of the average system for endpoints $J_0 = I_0$, and $J_f = I_f$. Actually, we should have solved for

$$J_0 = I_0 - \varepsilon \nu_I(I_0, p_{J0}, \varphi_0), \quad J_f = I_f - \varepsilon \nu_I(I_f, p_J(s_f), \varphi_f)$$

where $(s_f, p_{J0})$ is our solution of the average system. Computing the sensitivities of the shooting function $S(s_f, p_{J0}, J_0, J_f, p_\psi)$ with respect to $J_0$ and $J_f$, an update of order one in $\varepsilon$ of its solution can be computed. In particular, $\Delta\psi$ is updated according to

$$\Delta\psi(I_0, I_f, p_\psi) + \varepsilon\delta\Delta\psi(I_0, I_f, p_\psi, \varphi_0, \varphi_f).$$

Note that while we kept implicit the dependence of the update $\delta\Delta\psi$ on the solution $(s_f, p_{J0})$ of the average boundary value problem, we made explicit its dependence on the particular choice of solution for the phases $\varphi_0$ and $\varphi_f$. All in all, assuming that $\Delta\psi/\varepsilon$ provides an estimate of $\varphi(t_f) - \varphi(0)$, we obtain the following equation, parameterized by $p_\psi$:

$$\frac{\Delta\psi(I_0, I_f, p_\psi)}{\varepsilon} = \varphi_f(p_\psi) - \varphi_0(p_\psi) - \delta\Delta\psi(I_0, I_f, p_\psi, \varphi_0(p_\psi), \varphi_f(p_\psi)).$$

### C. Methodology

We now propose an algorithm for the generation of reliable initial guesses for the original problem that leverage on solutions of the average system:

#### 1) Solve the average problem

$$\forall p_{\psi} \in [p_{\psi}^{min}, p_{\psi}^{max}] \text{ find } s_{f}^{p_{\psi}}, p_{J0}^{p_{\psi}} \text{ s.t.:} \\ S(s_{f}^{p_{\psi}}, p_{J0}^{p_{\psi}}, p_{\psi}) = 0$$

#### 2) Find compatible endpoint phases

$$\forall p_{\psi} \in [p_{\psi}^{min}, p_{\psi}^{max}] \text{ find } \varphi_{0}^{(j)}, \varphi_{f}^{(k)} \text{ s.t.:} \\ \nu_{p_{\psi}}(J_{0}, p_{J0}^{p_{\psi}}, \varphi_{0}^{(j)}) = -p_{\psi}, \quad j = 1, \dots, n_{0}^{p_{\psi}} \\ \nu_{p_{\psi}}(J_{0}, p_{J}(s_{f}^{p_{\psi}}), \varphi_{f}^{(k)}) = -p_{\psi}, \quad k = 1, \dots, n_{f}^{p_{\psi}}$$

We omit the dependence on $p_\psi$ of $\varphi_0^{(j)}, \varphi_f^{(k)}$, and all other variables defined hereafter to simplify the notation.

#### 3) Compute sensitivities

$$\forall p_{\psi} \in [p_{\psi}^{min}, p_{\psi}^{max}] \text{ find } \delta s_{f}^{(jk)}, \delta p_{J0}^{(jk)}, \delta \Delta \psi^{(jk)} \text{ s.t.:} \\ \frac{\partial S}{\partial s_f} \delta s_{f}^{(jk)} + \frac{\partial S}{\partial p_{J0}} \delta p_{J0}^{(jk)} = -\frac{\partial S}{\partial J_0} \delta J_0^{(j)} - \frac{\partial S}{\partial J_f} \delta J_f^{(k)} \\ \delta \Delta \psi^{(jk)} = \frac{\partial \Delta \psi}{\partial s_f} \delta s_{f}^{(jk)} + \frac{\partial \Delta \psi}{\partial p_{J0}} \delta p_{J0}^{(jk)} + \frac{\partial \Delta \psi}{\partial J_0} \delta J_0^{(j)} \\ j = 1, \dots, n_0^{p_{\psi}}, \quad k = 1, \dots, n_f^{p_{\psi}}$$

where $\delta J_0^{(j)}$ and $\delta J_f^{(k)}$ are evaluated by means of the near-identity transformation as

$$\delta J_0^{(jk)} = -\nu_I \left( I_0, p_{J0}^{p_{\psi}}, \varphi_0^{(j)} \right), \\ \delta J_f^{(k)} = -\nu_I \left( I_f, p_J(s_f^{p_{\psi}}), \varphi_f^{(k)} \right).$$

#### 4) Identify solutions at $\varepsilon$

Given $\varepsilon$, find $p_\psi$ s.t.:

$$\mathrm{mod} \left( \frac{\Delta\psi(I_0, I_f, p_\psi)}{\varepsilon}, 2\pi \right) = \varphi_f^{(k)} - \varphi_0^{(j)} - \delta\Delta\psi^{(jk)} \quad (5)$$

Values of $p_\psi$ solution of Eq. (5) and the associated $\varphi_0^j, p_{J0}^{p_\psi} + \varepsilon\delta p_{J0}^{(jk)}$, and $s_f^{p_\psi} + \varepsilon\delta s_f^{(jk)}$ can be used to initialize the shooting of the original problem.

We stress that all steps only use the average system and the near-identity transformation. The small parameter, $\varepsilon$, only appears at the left-hand side of Eq. (5), whereas the right-hand side is independent of $\varepsilon$.

Concerning the complexity of the methodology: Step 1 implies first the solution of one average problem for $p_\psi = 0$ (this is the only shooting problem that has to be solved without a priori knowledge). Then, differential continuation can be used to obtain the solution in the range $[p_\psi^{min}, p_\psi^{max}]$. Step 2 involves the zero search of algebraic equations. Step 3 requires the solution of $(n+1)$-dimensional linear systems. Finally, step 4 is the identification of intersections of a planar closed curve with a discontinuous function of $p_\psi$. Graphical insight in this last step is provided in Section IV.

### IV. CASE STUDY: TIME-OPTIMAL ORBITAL TRANSFER

Orbital mechanics provides an ideal case study for this approach. We consider a Keplerian and co-planar time-optimal low-thrust transfer. Slow variables are some set of orbital elements defining the instantaneous shape and orientation of the orbit. Slow elements used in this work are

$$I_1 = a, \quad I_2 = e \cos(\nu), \quad I_3 = e \sin(\nu),$$

where $a, e$, and $\nu$ denote the semi-major axis, eccentricity, and argument of periapsis of the orbit, respectively. The fast variable is an angle positioning the satellite on a precise point of the orbit defined by $I$. The mean longitude is used for this purpose. Variational equations of these variables are provided by [10][Chap. 10.4]

$$\dot{I}_1 = 2\varepsilon \sqrt{\frac{I_1^3}{\mu(1 - I_2^2 - I_3^2)}} \left[ (I_3 \sin \sigma - I_2 \cos \sigma) u_1 + (1 + I_2 \cos \sigma + I_3 \sin \sigma) u_2 \right], \\ \dot{I}_2 = \varepsilon \sqrt{\frac{I_1(1 - I_2^2 - I_3^2)}{\mu}} \left[ \sin \sigma u_1 + \frac{I_2 + (2 + I_2 \cos \sigma + I_3 \sin \sigma) \cos \sigma}{1 + I_2 \cos \sigma + I_3 \sin \sigma} u_2 \right],$$
<table><caption>TABLE I
PARAMETERS OF THE CASE STUDY.</caption><thead><tr><th colspan="2">Constants</th></tr></thead><tbody><tr><td>Gravitational parameter, μ</td><td>3.986 · 10<sup>14</sup> m<sup>3</sup></td></tr><tr><td>Nominal thrust-to-mass ratio, ε<sub>0</sub></td><td>10<sup>-4</sup> N/m</td></tr><tr><th colspan="2">Initial conditions</th></tr><tr><td>Semi-major axis, I<sub>1</sub></td><td>7000 · 10<sup>3</sup> km</td></tr><tr><td>Eccentricity vector x, I<sub>2</sub></td><td>0.01</td></tr><tr><td>Eccentricity vector y, I<sub>3</sub></td><td>0</td></tr><tr><th colspan="2">Final conditions</th></tr><tr><td>Semi-major axis, I<sub>1</sub></td><td>7500 · 10<sup>3</sup> km</td></tr><tr><td>Eccentricity vector x, I<sub>2</sub></td><td>0</td></tr><tr><td>Eccentricity vector y, I<sub>3</sub></td><td>0</td></tr></tbody></table>

$$ \dot{I}_3 = \varepsilon \sqrt{\frac{I_1 (1 - I_2^2 - I_3^2)}{\mu}} \left[ -\cos \sigma u_1 + \frac{I_3 + (2 + I_2 \cos \sigma + I_3 \sin \sigma) \sin \sigma}{1 + I_2 \cos \sigma + I_3 \sin \sigma} u_2 \right], \\ \dot{\varphi} = \sqrt{\frac{I_1^3}{\mu}} + \varepsilon \sqrt{\frac{I_1 (1 - I_2^2 - I_3^2)}{\mu}} \left[ \left( \frac{I_2 \cos \sigma + I_3 \sin \sigma}{1 + \sqrt{1 - I_2^2 - I_3^2}} + 2 \frac{\sqrt{1 - I_2^2 - I_3^2}}{1 + I_2 \cos \sigma + I_3 \sin \sigma} \right) u_1 + \frac{I_2 \cos \sigma - I_3 \sin \sigma}{1 + \sqrt{1 - I_2^2 - I_3^2}} \frac{2 + I_2 \cos \sigma + I_3 \sin \sigma}{1 + I_2 \cos \sigma + I_3 \sin \sigma} u_2 \right]. $$

where the small parameter, $\varepsilon$, is the thrust-to-mass ratio magnitude, $\mu$ is the planetary gravitational constant, and $\sigma$ is the true longitude. Conversion from true to mean longitude, $\varphi = \varphi(\sigma, I_2, I_3)$, is achieved by evaluating the Kepler's law for equinoctial elements. Table I lists the parameters of the TPBVP. Controls $u_1$ and $u_2$ are radial and tangential thrust, respectively.

Figure 3 provides a graphical interpretation of the fourth step of the methodology. The black-dashed curve is the right-hand side of the equation, which is independent of $\varepsilon$. Blue curves depict the left-hand term. Colored dots indicate solutions of the original problem. This shows that intersections are excellent initial guesses.

Because $\varepsilon$ appears at the denominator of the left-hand term, the blue curve seems to circulate from the bottom to the top when $\varepsilon$ is arguably large. Then, the number of solutions (i.e., intersections) increases when $\varepsilon$ is very small (the slope of the curve becomes more and more vertical), as shown in Fig. 4, where a very small $\varepsilon$ is used. We stress that each red dot of this figure denotes a converged solution of the original problem.

The quality of the estimation of the solution can be appreciated in Fig. 5. Here, solid lines are solutions of the original TPBVP obtained by numerical continuation on the parameter $\varepsilon$ (starting from $\varepsilon = \varepsilon_0$ we want to reduce it). The Hampath software described in [11] is used for this purpose. Dots denote the approximation achieved by means of our approach. The continuation parameter is not monotonically

<figure><img src="image_3.png" /><figcaption>Fig. 3. Graphical interpretation of the 4th point of the methodology (solution of Eq. (5)). Black-dashed lines are the right-hand side of the equation, which is independent of $\varepsilon$. Solid-blue curve is the left-hand side, $\Delta\psi/\varepsilon$. Large dots are the solutions of the original problem.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 4. Solutions found for $\varepsilon = 0.08\varepsilon_0$. Refer to the caption of Fig. 3 for further details.</figcaption></figure>

decreasing as shown in the close up of Fig. 5(b). We note that this zoom is centered on the value $\varepsilon = 0.995\varepsilon_0$, which is used in Fig. 3(b). Here, four intersections are present, which correspond to the four solutions obtained with the numerical continuation for that specific value of $\varepsilon$. Two new solutions are generated when $\Delta\psi/\varepsilon$ is tangent to the right-hand side curve of Eq. (5). We note that colors of the two branches are consistent with the colors of dots in Fig. 3.

Finally, Fig 6 depicts the trajectory of the orbital eccentricity (i.e., $\sqrt{I_2^2 + I_3^2}$) for $\varepsilon = 10\varepsilon_0$. Both solutions of the original TPBVP and of the average counterpart are
<figure><img src="image_5.png" /><figcaption>Fig. 5. Convergence of the solution with respect to $\epsilon$. Solid lines are solutions of the original problem. Black dots denote the estimation of the solutions obtained by means of the proposed algorithm.</figcaption></figure>

depicted. Endpoints of the three trajectories coincide. The two solutions of the original problem differ in their boundary phases $\varphi_0$ and $\varphi_f$. We note that the black dashed trajectory is not the one used to generate the initial guess used to solve the shooting of the original problems. Conversely, initial guesses were generated by using average trajectories with endpoints slightly modified by means of Eq. 4. This correction is mandatory to avoid the deterioration of the guess which was pointed out in Fig. 2.

## V. CONCLUSION

This paper offered a methodology for the accurate approximation of time-optimal maneuvers of fast-oscillating systems. The practical interest of this contribution is the relationship that was established between the magnitude of the small parameter $\epsilon$ (read control) and the boundary phases, $\varphi_0$ and $\varphi_f$. Specifically, the number of possible solutions and their quantitative assessment for a specific $\epsilon$ can be estimated using only the information given by the average solution. In addition, the proposed methodology provides insight into how new combinations of $\varphi_0$ and $\varphi_f$ (i.e., multiple potential local optima) are generated when $\epsilon$ is decreased. This qualitative behavior cannot be fully understood by inspecting solutions of the original problem owing to numerical difficulties related to the integration of

<figure><img src="image_6.png" /><figcaption>Fig. 6. Trajectory of the orbital eccentricity for $\epsilon = 10\epsilon_0$. Solid lines are solutions of the original problem. The black-dashed line is the solution of the average TPBVP for $p_{\psi} = 0$.</figcaption></figure>

trajectories for very low $\epsilon$.

**Acknowledgements.** The authors thank Richard Epenoy (CNES) and Thierry Dargent (Thales Alenia Space) for many fruitful discussions.

## REFERENCES

[1] J. A. Sanders, F. Verhulst, and J. Murdock, *Averaging methods in nonlinear dynamical systems*, vol. 59 of *Applied Mathematical Sciences*. New York: Springer, second ed., 2007.

[2] V. I. Arnold, *Mathematical methods of classical mechanics*, vol. 60 of *Graduate Texts in Mathematics*. New York: Springer-Verlag, 2nd ed., 1989. Translated from Russian by K. Vogtmann and A. Weinstein.

[3] T. N. Edelbaum, “Optimum power-limited orbit transfer in strong gravity fields,” *AIAA J.*, vol. 3, pp. 921–925, 1965.

[4] S. Geffroy, *Généralisation des techniques de moyennation en contrôle optimal – Application aux problèmes de rendez-vous orbitaux en poussée faible*. Thèse de doctorat, Institut National Polytechnique de Toulouse, Toulouse, France, Oct. 1997.

[5] T. Dargent, “Averaging technique in T_3D an integrated tool for continuous thrust optimal control in orbit transfers,” in *4th AAS/AIAA Space Flight Mechanics Meeting*, (Santa Fe, New Mexico), Jan. 2014.

[6] B. Bonnard and J.-B. Caillau, “Riemannian metric of the averaged energy minimization problem in orbital transfer with low thrust,” *Ann. Inst. H. Poincaré Anal. Non Linéaire*, vol. 24, no. 3, pp. 395–411, 2007.

[7] A. Bombrun and J.-B. Pomet, “Asymptotic behavior of the time optimal orbital transfer for low thrust 2-body control system,” *DCDS supplements*, pp. 122–129, Sept. 2007.

[8] F. Chaplais, “Averaging and deterministic optimal control,” *SIAM J. Control Optim.*, vol. 25, no. 3, pp. 767–780, 1987.

[9] J. Rouot, *Méthodes géométriques et numériques en contrôle optimal et applications au transfert orbital à poussée faible et à la nage à faible nombre de Reynolds*. Thèse de doctorat, Univ. de Nice - Sophia Antipolis, Nov. 2016.

[10] R. H. Battin, *An Introduction to the Mathematics and Methods of Astrodynamics, Revised Edition*. American Institute of Aeronautics and Astronautics (AIAA), Jan. 1999.

[11] J.-B. Caillau, O. Cots, and J. Gergaud, “Differential continuation for regular optimal control problems,” *Optimization Methods and Software*, vol. 27, no. 2, pp. 177–196, 2012.