# Efficient Numerical Optimal Control for Highly Oscillatory Systems

Jakob Harzer, Jochem De Schutter, Moritz Diehl This research was supported by
DFG via Research Unit FOR 2401 and project 424107692 and by the EU via ELO-X
953348. Jakob Harzer and Jochem De Schutter are with the Department of
Microsystems Engineering (IMTEK) and Moritz Diehl is with the Department of
Microsystems Engineering (IMTEK) and Department of Mathematics, University of
Freiburg, Georges-Koehler-Allee 102, 79110 Freiburg, Germany. Emails:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We present an efficient transcription method for highly oscillatory optimal
control problems. For these problems, the optimal state trajectory consists of
fast oscillations that change slowly over the time horizon. Out of a large
number of oscillations, we only simulate a subset to approximate the slow
change by constructing a semi-explicit differential-algebraic equation that
can be integrated with integration steps much larger than one period. For the
solution of optimal control problems with direct methods, we provide a way to
parametrize and regularize the controls. Finally, we utilize the method to
find a fuel-optimal orbit transfer of a low-thrust satellite. Using the novel
method, we reduce the size of the resulting nonlinear program by more than one
order of magnitude.

##  I Introduction

Oscillating motion appears in a wide range of physical systems, from fast
oscillations in crystal oscillators, over the forced vibrations of an audio
speaker, up to the very slow tides. Although we might not expect it to, the
oscillating motion might not be periodic but also include a slow change, that,
compared to the fast oscillations, is barely noticeable. To notice this
change, one would have to observe the motion over a long horizon and possibly
a large number of oscillations.

We will call this behavior ‘highly oscillatory’ and we define ‘highly
oscillatory systems’ following [1] as “systems whose solutions may be
oscillatory in the sense that there is a fast solution which varies regularly
about a slow solution” under the condition that the timescale of the fast
solution is much shorter than the interval of interest. In the same fashion,
we define the property ‘highly oscillatory’ for a trajectory of a controlled
system.

When simulating highly oscillatory systems, a sufficiently small integration
step size is needed to simulate the system accurately over a single cycle. But
since we are interested in the long-term behavior over a large number of
cycles, the number of integration steps required becomes excessively large. We
can transfer this definition also to optimal control problems (OCP). A
continuous-time OCP is now called ‘highly oscillatory’ if the optimal state
and control trajectories are expected to be highly oscillatory.

When a highly oscillatory OCP is solved using direct multiple shooting, a
large number of steps and controls corresponds to a large number of variables
and constraints in the resulting nonlinear program (NLP), which in consequence
requires huge computational effort to solve.

The question we want to answer in this paper is:

> “How can we exploit the fact that the expected optimal trajectory of a
> system in an optimal control problem is highly oscillatory in order to speed
> up the optimization process?”

Consider the simulation of a system
$f:{\mathbb{R}}^{n_{x}}\times{\mathbb{R}}^{n_{u}}\rightarrow{\mathbb{R}}^{n_{x}}$.
Suppose $x(t)\in{\mathbb{R}}^{n_{x}},\ t\in[0,t_{\mathrm{f}}]$ is the highly
oscillatory solution of the initial value problem

$\displaystyle\dot{x}(t)$ $\displaystyle=f(x(t),u(t))$ (1a) $\displaystyle
x(0)$ $\displaystyle=x_{0}$ (1b)

where $x_{0}\in{\mathbb{R}}^{n_{x}}$, as visualized in Fig. 1. Let us assume
for a moment that the controls are fixed, that we know the period of the
oscillations $T$, and that we have a one-period map
$\Psi:{\mathbb{R}}^{n_{x}}\rightarrow{\mathbb{R}}^{n_{x}}$ that computes the
state $x_{k+1}=\Psi(x_{k})$ after one period from a given $x_{k}$.

Figure 1: Highly oscillatory trajectory and its envelope

Since the cycles only change slowly, the map $\Psi$ is a near-to-identity map
[2]. Starting from $x_{0}$ we are able to generate a series of points
$x_{1},x_{2},\dots$ that samples the trajectory periodically. In the manner of
[3], a trajectory that smoothly interpolates this slowly changing series is
called the ‘envelope’ $x_{\mathrm{e}}(t)$. We assume that $x_{\mathrm{e}}(t)$
is ‘sufficiently smooth’ in the sense that its derivatives in $t$ exist and
are bounded, and that it stems from some envelope dynamics

$\displaystyle\dot{x}_{\mathrm{e}}(t)=f_{e}(x_{\mathrm{e}}(t))$ (2)

whose solution $x_{\mathrm{e}}(t)$ is equal to the highly oscillatory solution
$x(t)$ at multiples of the period $kT,k\in\mathbb{N}$. More details can be
found in [2]. Since these envelope dynamics are assumed to generate a slowly
changing trajectory, we can integrate them with stepsizes much bigger than a
single cycle duration, while still retaining a decent accuracy.

###  I-A Contribution

In this work, we construct a differential-algebraic equation (DAE) that
approximates the envelope dynamics of the highly oscillatory system. Instead
of simulating all cycles of an oscillating trajectory, we only simulate a few
to gather information about the slow change that occurs over the time horizon.
To account for unknown and possibly changing cycle durations $T$ of the
oscillations, we will scale the independent variable $t$ and use phase
conditions to define a cycle. We provide a way to parametrize the controls and
use the method to efficiently solve highly oscillatory OCPs using direct
methods. This method can significantly reduce the number of variables in the
resulting NLP, but requires careful regularization of the state trajectories.

This simulation approach is closely related to the stroboscopic method for
highly oscillatory problems [2] where the envelope dynamics are also modeled
with a differential equation and then numerically approximated. Quite similar
as well are the multirevolution methods[4] and the envelope following methods
[3, 1, 4, 5] that instead of approximating the continuous system (2), model
the envelope dynamics using a difference equation.

We propose an envelope-following formalism specifically tailored for direct
optimal control. Whereas averaging methods have been used for highly
oscillatory optimal control before, such as in [6], to the authors’ knowledge,
ideas from envelope following methods, on the other hand, have not been used
for optimal control yet.

###  I-B Outline

After introducing a guiding example, we start with the definition of a single
cycle of the oscillations in Section II. In Section III we construct an
$N$-cycle OCP, that we then simplify by approximating the envelope dynamics
using a DAE in Section IV. We discuss the integration of the DAE and problems
of this approach in Section V. In Section VI we discuss a numerical example on
the model of a low-thrust satellite, and conclude the paper with Section VII.

All the coding of this work is done in Python 3.9, we use CasADi [7] to set up
the optimization problems and solve the nonlinear programs using IPOPT [8]
equipped with the linear solver MA27[9].

###  I-C Guiding Example

As a guiding example, we consider the nonlinear predator-prey system, also
known as the Lotka-Volterra equations. The system models the interaction
between a prey population of size $r$ and a predator population of size $s$.
We modify the original system by adding a control $u\in{\mathbb{R}}$ to
influence the growth of the predator species. The model is then given by

$\displaystyle x:=\begin{bmatrix}r\\\ s\end{bmatrix}\in{\mathbb{R}}^{2},$
$\displaystyle f(x,u):=\begin{bmatrix}\frac{2}{3}r-\frac{4}{3}rs\\\
rs-s\end{bmatrix}+\begin{bmatrix}0\\\ u\end{bmatrix}.$ (3)

For visualization, Fig. 2 shows a trajectory from the initial state
$\bar{x}_{0}=[1,1.5]^{\top}$ with constant control $u(t)=0.03$.

##  II Definition of a single cycle

We discretize the oscillatory trajectory into individual cycles by introducing
two affine phase conditions

$\displaystyle
q^{\top}x^{\scalebox{0.6}[0.7]{$-$}}=b^{\scalebox{0.6}[0.7]{$-$}}$ (4a)
$\displaystyle q^{\top}x^{\scalebox{0.7}{$+$}}=b^{\scalebox{0.7}{$+$}}$ (4b)

for the starting point $x^{\scalebox{0.6}[0.7]{$-$}}\in{\mathbb{R}}^{n_{x}}$
and end point $x^{\scalebox{0.7}{$+$}}\in{\mathbb{R}}^{n_{x}}$ of a cycle,
respectively, that are defined by a normalized vector
$q\in{\mathbb{R}}^{n_{x}}$ that satisfies $\|q\|=1$, and scalars
$b^{\scalebox{0.6}[0.7]{$-$}},b^{\scalebox{0.7}{$+$}}\in{\mathbb{R}}$. These
phase conditions have to be chosen depending on the system, and usually
involve prior knowledge of the expected oscillatory trajectories. We
distinguish two cases:

  1. (A)

The trajectory $x(t)$ is discretized by giving a condition on an oscillatory
state. For these types of systems we can choose
$b^{\scalebox{0.6}[0.7]{$-$}}=b^{\scalebox{0.7}{$+$}}$, which is equivalent to
defining a linear Poincaré section. For example, for the predator-prey system
(3) we will specify that “for the starting and end point of a cycle the first
state is 1” as the phase conditions

$\displaystyle q^{\top}$ $\displaystyle x^{\scalebox{0.6}[0.7]{$-$}}-1=0$ (5a)
$\displaystyle q^{\top}$ $\displaystyle x^{\scalebox{0.7}{$+$}}-1=0,$ (5b)

for which we choose $q=[1,\ 0]^{\top}$ and
$b^{\scalebox{0.7}{$+$}}=b^{\scalebox{0.6}[0.7]{$-$}}=1$. This condition is
visualized in Fig. 2.

Figure 2: Phase conditions (5) of type (A) discretize the state trajectory of
the Predator-Prey system into individual cycles.

  2. (B)

The trajectory $x(t)$ is discretized by giving a condition on a non-
oscillatory state, typically a ‘phase state’ such as a clockstate or a state
for a rotation angle. This is required if this state influences the dynamics
strongly and periodically over a single cycle, cf. [5, 2]. This approach is
used later in the example in Section VI.

Similar phase conditions have also been used by others to define a cycle of
the oscillations [10, 11, 12, 3, 5, 13].

To later integrate over $N$ cycles and also to handle changing cycle
durations, we introduce the numerical time $\tau\in[0,N]$, similar to [3, 1].
The time $t\in[0,t_{\mathrm{f}}]$ will now be called the ‘physical time’ and
can be included as an additional state in the model. We define, by a slight
abuse of notation, $x(\tau):=x(t(\tau))$. The relationship between $t$ and
$\tau$ is given by

$\displaystyle\frac{\ \mathrm{d}t(\tau)}{\ \mathrm{d}\tau}=T(\tau),$
$\displaystyle T(\tau)\coloneqq T_{k},\ \tau\in[k,k+1),$ (6)

where $T_{k}$ is the duration of the $k$-th cycle. This way, expressed in
numerical time, the same control grid can be applied to each cycle.

Let us parametrize the controls over one cycle by a matrix

$\displaystyle
U_{k}=[u_{k,0},u_{k,1},\dots,u_{k,N_{\mathrm{ctr}}-1}]\in{\mathbb{R}}^{n_{u}\times
N_{\mathrm{ctr}}}$ (7)

and apply $u(\tau)=\tilde{u}(\tau;U_{k})$ for $\tau\in[k,k+1)$. For example,
as depicted in Fig. 3, for the cycle $k$, this could be a piecewise constant
control parameterization

$\displaystyle\tilde{u}(\tau;U_{k})=u_{k,m},\ \text{ for }\ \tau\ \text{mod}\
k\in\left[\tfrac{m}{N_{\mathrm{ctr}}},\tfrac{m+1}{N_{\mathrm{ctr}}}\right),$
(8)

and for $m=0,..,N_{\mathrm{ctr}}-1$. A single cycle is divided into
$N_{\mathrm{ctr}}$ equidistant intervals and the $m$-th control is applied
over the $m$-th interval. The scaled dynamics for a cycle with control
parameterization $U_{k}$ now read

$\displaystyle\frac{\ \mathrm{d}x}{\
\mathrm{d}\tau}(\tau)=T(\tau)f(x(\tau),\tilde{u}(\tau;U_{k})),\qquad\tau\in[k,k+1).$
(9)

Let $F$ be a numerical integration routine that simulates these dynamics over
the numerical cycle duration $\Delta\tau=1$. In the fashion of [2], we call
this integrator the ‘micro-integrator’.

To numerically simulate a single cycle with a given control $U_{k}$ that
starts at a point $x_{k}$, we solve

$\displaystyle 0$
$\displaystyle=Q^{\top}(x^{\scalebox{0.6}[0.7]{$-$}}_{k}-x_{k})$ (10a)
$\displaystyle 0$ $\displaystyle=C(z_{k},U_{k}),$ (10b)

where the matrix $Q$ is defined such that $\bar{Q}=[q|Q]$ is an orthonormal
basis of the state space.

Here, the algebraic equations

$\displaystyle
C(z,U)\coloneqq\begin{cases}0&=q^{\top}x^{\scalebox{0.6}[0.7]{$-$}}-b^{\scalebox{0.6}[0.7]{$-$}}\\\
0&=x^{\scalebox{0.7}{$+$}}-F(x^{\scalebox{0.6}[0.7]{$-$}},U,T)\\\
0&=q^{\top}x^{\scalebox{0.7}{$+$}}-b^{\scalebox{0.7}{$+$}}\end{cases}\qquad\in{\mathbb{R}}^{n_{x}+2}$
(11)

define one cycle of duration $T$ with starting point
$x^{\scalebox{0.6}[0.7]{$-$}}$ and end point $x^{\scalebox{0.7}{$+$}}$ where
the variables are collected in
$z\coloneqq(x^{\scalebox{0.6}[0.7]{$-$}},x^{\scalebox{0.7}{$+$}},T)$. These
equations can be solved numerically using for example a Newton-type root-
finding method. By construction system (11) has multiple periodic solution,
hence a sufficiently good initial guess is required to simulate exactly one
cycle.

##  III $N$-Cycle Optimal Control Problem

By discretizing the state and control trajectories as explained above, we can
formulate an optimal control problem over a number of $N$ cycles as
w∑_i=0^N-1 L(z_k,U_k) + E(x^ $+$ N-1) 0= Q^⊤(x^ $-$ _0 - x_0 ) 0= C(z_k,U_k) &
k = 0,..,N-1 0≤H(z_k,U_k) & k = 0,..,N-1 0= Q^⊤(x^ $-$ _k - x^ $+$ _k-1) & k =
1,..,N-1 where $w=(z_{0},\dots,z_{N-1},U_{0}\dots,U_{N-1})$, $L$ is a cycle
cost and $E$ a terminal cost. Essentially, we simulate $N$ cycles with (III)
that are constrained to start at the initial point by (III) and are connected
by (III). Equation (III) enforces state and control constraints. Fig. 3
visualizes an example of such an $N$-cycle OCP.

Figure 3: Exemplary optimal solution of a 5-cycles OCP of the predator-prey
system (3).

As mentioned in the introduction, solving such a problem becomes increasingly
costly for larger $N$.

An alternative control parameterization that accommodates the highly
oscillatory nature of the system, derives the cycle controls from a low-order
polynomial

$\displaystyle U(\tau;U_{c})=\sum_{i=1}^{d_{u}}U_{\mathrm{c},i}\ell_{i}(\tau)$
(12)

that we evaluate at the cycle $k$ as $U_{k}:=U(k+\frac{1}{2};U_{c})$. We
construct this polynomial from $d_{u}$ construction points where $\ell_{i}$
are the orthogonal basis polynomials built from the construction times
$\tau_{\mathrm{c},1},\dots,\tau_{\mathrm{c},d_{u}}$. The construction controls
$U_{c}:=(U_{\mathrm{c},1},\dots,U_{\mathrm{c},d_{u}})\in{\mathbb{R}}^{d_{u}\times
n_{u}\times N_{\mathrm{ctr}}}$ are included as optimization variables. For a
constant periodic control parameterization $U(\tau)=U_{\mathrm{const}}$,
$d_{u}=1$, i.e. the same periodic control is then applied in each cycle.

##  IV Envelope Slope Approximation with a DAE

The core idea of this work is that instead of integrating the oscillating ODE
(1a) with a small stepsize, we integrate the DAE

$\displaystyle\frac{\ \mathrm{d}x}{\ \mathrm{d}\tau}$
$\displaystyle=f_{\mathrm{e}}(x(\tau),z(\tau),U(\tau))$ (13a) $\displaystyle
0$ $\displaystyle=g_{\mathrm{e}}(x(\tau),z(\tau),U(\tau))$ (13b)

that approximates the envelope dynamics with a large stepsize. In this paper,
we present the ‘central differences’ DAE that approximates the dynamics of the
envelope at the point $x$ as

$\displaystyle
f_{\mathrm{e}}(x,z):=x^{\scalebox{0.7}{$+$}}-x^{\scalebox{0.6}[0.7]{$-$}}$
(envelope approx.) (14a) with the algebraic equations $g_{\mathrm{e}}(x,z)$
given by $\displaystyle 0$
$\displaystyle=Q^{\top}\left(\frac{x^{\scalebox{0.7}{$+$}}+x^{\scalebox{0.6}[0.7]{$-$}}}{2}-x\right)$
(connecting cond.) (14b) $\displaystyle 0$ $\displaystyle=C(z,U)$ (cycle
cond.). (14c)

Essentially, at a point $x$, we simulate one cycle around it to estimate the
slow change of the highly oscillatory trajectory. Algebraic equation (14c)
corresponds to the micro-integration of a single cycle. The approximation of
the envelope dynamics is inspired by a central difference scheme and is
sketched in Fig. 4. The differential state $x$ is constrained to ‘lie in the
middle’ of the algebraic states $x^{\scalebox{0.6}[0.7]{$-$}}$ and
$x^{\scalebox{0.7}{$+$}}$ by the connecting condition (14b).

Figure 4: The ’Central Difference’ DAE approximates the dynamics of the
envelope of the state $s$ of the predator-prey model.

Multiplying by $Q^{\top}$ in (14b) projects the condition into the null space
of $q$. This prevents the DAE from being overdetermined since the direction
$q$ is already constrained by the phase conditions inside (14c).

In a slightly different formulation, a similar method was already used in [10]
to simulate electrical oscillators where the authors employ forward-difference
and backward-difference schemes and also suggest using a central-difference
scheme. Also noteworthy is the approach by [2] where the authors approximate
the envelope dynamics using the previous and the next cyclic point, that in
consequence requires the simulation of two cycles.

##  V Macro-Integration of the Envelope DAE

For the simulation of the oscillatory system over a timescale much larger than
the duration of a single cycle, instead of simulating a large number of $N$
cycles in detail, we now just integrate the DAE whose solution approximates
the envelope of the oscillations. Again, in the manner of [2], this
integration procedure will be called the ‘macro-integration’. Since we now use
the numerical time $\tau$ instead of the physical time, integrating the DAE
with a stepsize of $N$ simulates the highly oscillatory dynamics over $N$
cycles. In practice, for the macro-integration of the DAE we use implicit RK
methods.

Similar to [14], we can list the following three sources of error for the
overall simulation procedure:

  1. 1.

Errors of the micro-integration over a single cycle.

  2. 2.

Errors of the central difference approximation of the envelope dynamics.

  3. 3.

Errors of the macro-integration of the constructed envelope DAE.

To analyze these errors analytically is outside the scope of this paper.

###  V-A Simulation Error and Regularization

The simulation method presented here is tailored to highly oscillatory
problems where we expect that the trajectories are composed of fast
oscillations that change very little over the time horizon. Thus, we expect
the cycle duration to only drift slowly, the shape of the cycles to remain
similar and the envelope to be sufficiently smooth. If this is not the case,
we introduce a large simulation error contribution from error source 2.

Dependent on the cost functions $L$ and $E$, the optimizer might prefer
trajectories with large values in the higher-order derivatives. This is
especially a problem for the presented method that becomes inaccurate if the
envelope is not sufficiently smooth. One way to ensure this is to add
additional constraints or regularization terms to the objective function that
penalize the derivatives of the envelope or control trajectories in the macro-
integration.

###  V-B Regularization Experiments

Suppose we want to drive the predator prey system (3) from an initial state
$\bar{x}_{0}=[1,2]$ into a target state $\bar{x}_{\mathrm{target}}=[1,0.7]$
within $N=20$ cycles that are characterized by the phase conditions (4a). The
controls are very small and bounded by

$\displaystyle-0.005\leq u(\tau)\leq 0.005,\forall\tau,$ (15)

such that even with an optimal control, the system cannot be driven to the
target point. Instead of solving the $N$-cycle OCP (III), we instead solve an
OCP  x_e(⋅),z_e(⋅),U_c∥x_e(N) - ¯x_N∥^2_2 + α∫_0^N ∥ d3xdτ3(τ) ∥^2_2 dτ 0=
x_e(0) - ¯x_0 ˙x_e= f_e(x_e(τ),z_e(τ),U(τ; U_c)), & τ∈[0,N] 0=
g_e(x_e(τ),z_e(τ),U(τ; U_c)), & τ∈[0,N] 0≤h(x_e(τ),z_e(τ),U(τ; U_c)), &
τ∈[0,N] that uses the DAE dynamics. We search for an optimal polynomial
control $U(\tau;U_{c})$ of degree 2, i.e. $d_{u}=3$, with
$N_{\mathrm{ctrl}}=3$ controls per cycle. In each micro-integration, we chain
10 RK4 integration steps to simulate the model dynamics over each of the 3
control intervals, i.e. we perform 30 RK4 steps per cycle. To keep the
envelope sufficiently smooth, we regularize the third derivative of the
envelope with a small positive weight $\alpha>0$. We discretize the OCP using
a single interval of a 3-stage Gauss-Legendre collocation scheme of order 6
and solve the NLP using the NLP solver IPOPT[8].

Figure 5: Found Solution to the OCP with Regularization $\alpha=2\cdot
10^{-5}$ Figure 6: Found Solution to the OCP without Regularization $\alpha=0$

Fig. 5 shows the optimal trajectory of the envelope in the states $r$ and $s$,
the micro-integrations, as well as an a-posteriori performed simulation of the
systems with the found optimal controls. With $\alpha=2\cdot 10^{-5}$ we
observe that the trajectory and endpoint of the envelope matches the
simulation very well.

In contrast, Fig. 6 shows the found trajectories for $\alpha=0$. The computed
collocation trajectory reaches the target point $\bar{x}_{\mathrm{target}}$,
the tracking objective is thus zero, but the endpoint of the exact simulation
$x_{\mathrm{end},\mathrm{sim}}$ does not at all resemble the integration
endpoint $x_{\mathrm{end}}$. This is because the optimizer chooses the
optimization variables such that the envelope varies strongly. This then
results in a large simulation error of the method, which then drives the
integration endpoint to the target. The optimizer ‘used’ the integration error
to minimize the objective by driving the trajectories out of the allowed
regions.

The choice of the regularization factor is important since it is a trade-off
between smoother solutions that are not prone to the integration error and the
original optimization objective. To investigate this, we solve the OCP for a
range of $\alpha\in[10^{-7},10^{-2}]$. As a baseline, we use the solution of
the full OCP (III) with the same polynomial control parametrization but
without the regularization term whose optimal objective can be interpreted as
the lower bound.

Figure 7: Tracking objective of the DAE OCP compared to an exact simulation
result for the optimal controls for different $\alpha$.

Fig. 7 shows the relationship between the choice of the regularization factor
$\alpha$ and the discrepancy between the ‘optimal’ tracking objective and the
simulated a-posteriori tracking objective.

We observe that for large $\alpha$ the regularization impacts the optimal
trajectory considerably, such that the tracking objective is larger than
necessary. For decreasing $\alpha$, the objective of the DAE OCP approaches
the baseline objective; there is a ‘sweet-spot’ region where the regularized
and the full exact OCP coincide. If $\alpha$ is too small, the DAE OCP,
through the erroneous integration explained above, converges to even lower
objectives. In this region, the simulation does not match the found envelope
trajectories.

##  VI Fuel-Optimal Orbit Transfers of a Low-Thrust Satellite

As a second example, we consider a problem similar to the problem formulated
by [15]. We want to find the optimal thrust direction and magnitude to
transfer a satellite from a circular park orbit at
$500\text{\,}\mathrm{k}\mathrm{m}$ with an inclination of
$87.4$\mathrm{\SIUnitSymbolDegree}$$ to a target mission orbit at
$1400\text{\,}\mathrm{k}\mathrm{m}$ with the same inclination but an
eccentricity of $0.0011$. Since the satellite has a very-low thrust engine, it
requires a large number of orbits to reach the target orbit. The satellite has
an initial mass of $m_{0}=$160\text{\,}\mathrm{k}\mathrm{g}$$, a maximum
thrust of $T_{0}=$0.2\text{\,}\mathrm{N}$$ and a specific impulse of
$I_{\mathrm{sp}}=$1600\text{\,}\mathrm{s}$$, as in [15]. We model the
satellite’s state using modified equinoctial elements $p,f,g,h,k,L$, where
only the semi-late rectum $p$ and the true longitude $L$ have intuitive
physical meaning[6]. Their relationship to classical orbit elements is given
in [15]. Also including a state for the satellite mass $m$ as well as a
clockstate $t$, we obtain the state vector and control vector

$\displaystyle x=[p,f,g,h,k,L,m,t]\in{\mathbb{R}}^{8},$ $\displaystyle
u\in{\mathbb{R}}^{3}.$ (16)

The state dynamics are given by

$\displaystyle f(x,u):=\begin{bmatrix}A(x)\Delta(x,u)+b(x)\\\
-\frac{T_{0}}{g_{\mathrm{e}}I_{\mathrm{sp}}}\sqrt{u^{\top}u+\epsilon^{2}}\\\
1\end{bmatrix}$ (17)

where $g_{e}=$9.81\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-2}$$ is the
gravitational acceleration at sea level[6]. Information on the matrices $A$
and $b$ can be found in [15]. The perturbation acceleration

$\displaystyle\Delta(x,u)=\Delta_{\mathrm{G}}(x)+\frac{T_{0}}{m}u$ (18)

contains a term for the gravitational perturbations $\Delta_{\mathrm{G}}$ as
described by [15]. These perturbations are due to earth’s oblateness, they
dominate the satellite’s thrust[15] and give rise to the highly oscillatory
state trajectory. The acceleration by the thruster is characterized by the
thrust vector $u=[u_{1},u_{2},u_{3}]^{\top}\in{\mathbb{R}}^{3}$ that satisfies

$\displaystyle-1\leq$ $\displaystyle u_{i}$ $\displaystyle\leq 1,\qquad
i=1,2,3,$ (19a) $\displaystyle u^{\top}u$ $\displaystyle\leq 1.$ (19b)

For numerical reasons, we include a small parameter $\epsilon=10^{-5}$ to
smoothen the nonlinear square root in (17). We constrain the classical
elements at the endpoint of the trajectory to be within certain bounds of the
target values, as explained in [15].

We discretize the trajectory into orbits with the phase conditions

$\displaystyle L^{\scalebox{0.6}[0.7]{$-$}}$ $\displaystyle=L_{0}$ (20a)
$\displaystyle L^{\scalebox{0.7}{$+$}}$ $\displaystyle=L_{0}+2\pi$ (20b)

on the phase state $L$, where $L_{0}$ is the initial value of $L$ and
$q,b^{\scalebox{0.6}[0.7]{$-$}},b^{\scalebox{0.7}{$+$}}$ are chosen
accordingly. This type of conditions (B) is necessary since the dynamics (17)
are $2\pi$-periodic in the state $L$.

For the micro-integration of an orbit, we divide it into
$N_{\mathrm{ctrl}}=30$ piecewise constant control intervals that each
integrate the dynamics using one step of a 3-stage Gauss-Legendre (GL)
collocation scheme of order 6, which results in
$U(\tau)\in{\mathbb{R}}^{3\times 30}$. The algebraic variables now also
contain the variables of the collocation scheme.

We search for a fuel optimal transfer over $578$ orbits, that we split into
two intervals of $N=289$ orbits each. In each interval, we integrate the DAE
dynamics using a 5-stage GL collocation, and use a constant parametrization of
the controls $U(\tau)$ with $d_{u}=1$. The macro-integration variables are
initialized by a linear interpolation between the initial state and target
final state, the micro-integration variables by simulating a single orbit at
the interpolation point. The control variables for all intervals are
initialized as $[0,\ 0.1,\ 0]^{\top}$.

We solve the NLP with a total of $17\ 334$ variables and $17\ 222$ constraints
in about 20 seconds. The results are presented in Fig. 8 and 9. For reference,
we simulate the found optimal controls with high accuracy. We observe that the
envelope agrees with the simulation over the orbits and at the endpoint, with
a maximum relative error of $8.85\cdot 10^{-5}$ in the state $g$. The optimal
final mass is $155.48\text{\,}\mathrm{k}\mathrm{g}$, which almost coincides
with the value found by [15]. Due to memory limitations, we were not able to
solve the 32-times bigger full optimization problem that for the chosen
control and state discretization has about $555\ 000$ variables.

Figure 8: Semi-major axis and eccentricity of a simulation of the found
optimal controls with the target value indicated by the dashed line Figure 9:
The highly oscillatory trajectory of the state $f$ and $L$ as computed by the
DAE method as well as a simulation of the found optimal controls. The state
$L$ is subject to the phase conditions (20), in this state the macro-
integration values are not tied to values of the start and endpoint of the
micro-integration.

##  VII Conclusion

We presented a novel method to treat highly oscillatory optimal control
problems. The state and control trajectories of the problem are discretized
into whole cycles using phase conditions. Instead of simulating all cycles, we
instead integrate a differential-algebraic equation that approximates the
envelope dynamics of the state trajectory. As a numerical example, we used the
presented method to efficiently optimize the orbit transfer of a low-thrust
satellite. Future work will adress a more rigourous analyis and theoretical
justification of the method as well as the exploration of more use-cases.

## References

  * [1] L. R. Petzold, L. O. Jay, and J. Yen, “Numerical solution of highly oscillatory ordinary differential equations,” Acta Numerica, vol. 6, p. 437–483, 1997.
  * [2] M. P. Calvo, P. Chartier, A. Murua, and J. M. Sanz-Serna, “A stroboscopic numerical method for highly oscillatory problems,” in Numerical Analysis of Multiscale Computations (B. Engquist, O. Runborg, and Y.-H. R. Tsai, eds.), (Berlin, Heidelberg), pp. 71–85, Springer Berlin Heidelberg, 2012\.
  * [3] L. R. Petzold, “An efficient numerical method for highly oscillatory ordinary differential equations,” SIAM Journal on Numerical Analysis, vol. 18, no. 3, pp. 455–479, 1981.
  * [4] O. Graf, “Multirevolution methods for orbit integration,” in Proceedings of the Conference on the Numerical Solution of Ordinary Differential Equations (D. G. Bettis, ed.), (Berlin, Heidelberg), pp. 471–490, Springer Berlin Heidelberg, 1974.
  * [5] K. A. Gallivan, “Detection and integration of oscillatory differential equations with initial stepsize, order and method selection,” 12 1980.
  * [6] G. Yang, “Direct optimization of low-thrust many-revolution earth-orbit transfers,” Chinese Journal of Aeronautics - CHIN J AERONAUT, vol. 22, pp. 426–433, 08 2009.
  * [7] J. A. E. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl, “CasADi – A software framework for nonlinear optimization and optimal control,” Mathematical Programming Computation, vol. 11, no. 1, pp. 1–36, 2019.
  * [8] A. Wächter and L. Biegler, “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” Mathematical programming, vol. 106, pp. 25–57, 03 2006.
  * [9] “Hsl. a collection of fortran codes for large scale scientific computation..”  http://www.hsl.rl.ac.uk/.
  * [10] T. Mei and J. Roychowdhury, “A time-domain oscillator envelope tracking algorithm employing dual phase conditions,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 27, no. 1, pp. 59–69, 2008.
  * [11] D. Linaro, D. d. Giudice, A. Brambilla, and F. Bizzarri, “Application of envelope-following techniques to the shooting method,” IEEE Open Journal of Circuits and Systems, vol. 1, pp. 22–33, 2020.
  * [12] P. Maffezzoni, “A versatile time-domain approach to simulate oscillators in rf circuits,” IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 56, no. 3, pp. 594–603, 2009.
  * [13] C. W. Gear, “Automatic detection and treatment of oscillatory and/or stiff ordinary differential equations,” in Numerical Integration of Differential Equations and Large Linear Systems (J. Hinze, ed.), (Berlin, Heidelberg), pp. 190–206, Springer Berlin Heidelberg, 1982.
  * [14] M. P. Calvo, P. Chartier, A. Murua, and J. Sanz-Serna, “Numerical stroboscopic averaging for odes and daes,” Applied Numerical Mathematics, vol. 61, 11 2010.
  * [15] J. T. Betts, “Very low-thrust trajectory optimization using a direct sqp method,” Journal of Computational and Applied Mathematics, vol. 120, no. 1, pp. 27–40, 2000.