A derivative-free approach to optimal control problems with a piecewise
constant Mayer cost function NovembreNovember 2021 XX C. Audet, P.-Y. Bouchet,
L. Bourdin Bouchet, Audet, Bourdin Hamel, Benoit, Karine Hébert (2021). Un
exemple de citation , _Journal of Journals_ , vol. X issue Y, p.
n-mhttps://www.gerad.ca/fr Mai 2021

Pierre-Yves Bouchet LABEL:affil:geradLABEL:affil:polymtl Charles Audet
LABEL:affil:geradLABEL:affil:polymtl Loïc Bourdin LABEL:affil:unilim

affil:geradGERAD, Montréal (Qc), Canada, H3T 1J4 affil:polymtlPolytechnique
Montréal, Montréal (Qc), Canada, H3T 1J4 affil:unilimXLIM Research Institute,
UMR CNRS 7252, University of Limoges, France.

<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
Abstract A piecewise constant Mayer cost function is used to model optimal
control problems in which the state space is partitioned into several regions,
each having its own Mayer cost value. In such a context, the standard
numerical methods used in optimal control theory naturally fail, due to the
discontinuities and the null gradients associated with the Mayer cost
function. In this paper an hybrid numerical method, based on both derivative-
free optimization and smooth optimization techniques, is proposed to solve
this class of problems. Numerical simulations are performed on some standard
control systems to show the efficiency of the hybrid method, where Nomad and
Ipopt are used as, respectively, derivative-free optimization and smooth
optimization solvers.

#### Keywords:

optimal control, piecewise constant Mayer cost function, derivative-free and
blackbox optimization.

Résumé Une fonction de coût de Mayer constante par morceaux est requise pour
correctement modéliser des problèmes de contrôle optimal dans lesquels
l’espace des états est partitonné en différentes régions ayant chacune une
valeur de coût de Mayer. Dans un tel contexte, les méthodes numériques
communes dans la théorie du contrôle optimal échouent naturellement, à cause
des discontinuités et des gradients nuls associés à la fonction de coût de
Mayer. Dans cet article, nous proposons une méthode numérique hybride,
exploitant des techniques d’optimisation lisse et d’optimisation sans
dérivées, pour traiter cette classe de problèmes. Nous présentons des
simulations numériques sur plusieurs systèmes de contrôle classiques pour
montrer l’efficacité de notre méthode hybride, en exploitant les solveur Ipopt
pour l’optimisation lisse et Nomad pour l’optimisation sans dérivées.

#### Mots clés :

contrôle optimal, fonction de coût de Mayer constante par morceaux,
optimisation sans dérivées et de boîte noire.

Work of the first author is supported by NSERC Canada Discovery Grant
2020-04448.

##  1 Introduction

###  1.1 Motivation

The present work was initially motivated by numerically solving Bolza optimal
control problems of the form

$\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}x\in\mathcal{AC}([0,T],\mathbb{R}^{n})\\\
u\in\mathcal{L}^{\infty}([0,T],\mathbb{R}^{m})\end{array}}{\operatorname*{minimize}}&\begin{array}[]{l}\left\lfloor\lVert
x(T)-x_{\mathrm{end}}\rVert\right\rfloor+\displaystyle\int_{0}^{T}\ell(t,x(t),u(t))\,dt\end{array}\\\\[23.68048pt]
\mbox{subject to}&\begin{array}[t]{ll}\dot{x}(t)=f(t,x(t),u(t)),&\mbox{a.e.\
}t\in[0,T],\\\\[3.0pt] x(0)=x_{0},\\\\[3.0pt] c(t,x(t),u(t))\leq
0,&\mbox{a.e.\ }t\in[0,T],\\\\[3.0pt] c_{T}(x(T))\leq 0\end{array}\end{array}$
($\mathrm{P}_{\mathrm{Motiv}}$)

where standard notations are used (recalled in Sections 3 and 4) and in which
the dynamics, the constraint functions and the Lagrange cost function are
smooth, but where the Mayer cost is discontinuous, being defined as the
truncation (involving the floor function
$\left\lfloor\cdot\right\rfloor:\mathbb{R}\to\mathbb{Z}$) of the distance
between the final state $x(T)$ and a reference target $x_{\mathrm{end}}$. Note
that the corresponding Mayer cost function is piecewise constant with integer
values. In this paper we call plateaus the subsets of $\mathbb{R}^{n}$ on
which a piecewise constant function from $\mathbb{R}^{n}$ to $\mathbb{R}$ has
a given value. As illustrated in Figure 1(a), in the two-dimensional case
$n=2$ the plateaus of the Mayer cost function of Problem
($\mathrm{P}_{\mathrm{Motiv}}$) are concentric rings centered at
$x_{\mathrm{end}}$. Such a Mayer cost function is used to model practical
problems in which a kind of discrete distance to a target is evaluated, for
example in precision sports (shooting, darts, archery, precision landing of a
parachute, etc.). Figure 1(b) illustrates several properties of these
problems. For example, two trajectories ending on the same plateau are only
discriminated by the Lagrange cost. Also, a trajectory with a high Lagrange
cost ending on a given plateau may have a smaller Bolza cost than another with
a low Lagrange cost but ending on a higher-value plateau. Finally, the two
trajectories ending on the lowest-cost plateau show that reaching
$x_{\mathrm{end}}$ as close as possible may not be the optimal solution. These
examples show that one cannot remove the truncation in the Mayer cost function
without altering (possibly significantly) the optimal solutions.

Figure 1: (a) Representation of the Mayer cost function of Problem
($\mathrm{P}_{\mathrm{Motiv}}$) in the two-dimensional case $n=2$ with
$x_{\mathrm{end}}=(8,8)$. (b) Artificial trajectories and their associated
costs.

In optimal control theory, the numerical methods are usually classified as
either direct or indirect methods [30, 36]. An indirect method relies on the
first-order necessary optimality conditions provided by the Pontryagin maximum
principle (PMP). Most of the literature on PMP concerns smooth problems (see,
e.g., [12, 13, 29]), but extensions exist for problems satisfying (only)
Lipschitz-continuity assumptions, by using generalized notions of gradient
(see, e.g., [14, 21, 37]). Some papers derive a PMP for problems with
discontinuities, such as [22] in the context of hybrid dynamics, or [9] in the
context of time crisis problems where discontinuities occur as an indicator
function in the Lagrange cost function. However, deriving a relevant version
of the PMP in our context is unlikely because optimal trajectories are
inclined to end on the boundaries of the plateaus of the Mayer cost function,
where gradients and generalized gradients are not defined. Our numerical
simulations provided in Section 4 confirm this phenomenon. Therefore, in this
paper, we rather focus on direct methods which rely on a full discretization
of the optimal control problem, via a Runge-Kutta method [31] for example, in
order to recover a finite-dimensional optimization problem that can be
numerically solved using a relevant optimization algorithm.

However, a full discretization of Problem ($\mathrm{P}_{\mathrm{Motiv}}$)
produces a finite-dimensional optimization problem which may preserve the
discontinuity of the Mayer cost function, while most optimization algorithms
in the literature are dedicated to smooth problems (see, e.g., [27]).
Nevertheless, one can easily design numerical approaches attempting to
overcome this obstacle. A first approach consists in considering a smooth
approximation of the discretized problem that can be numerically solved using
a smooth optimization (SO) solver such as Ipopt [38]. However smoothing the
piecewise constant Mayer cost function leads to a smooth function with
gradients which are null (or almost null) on most of the domain. Hence they do
not provide useful information to the SO algorithm which thus focuses only on
the optimization of the Lagrange cost, while the latter may be less rewarding
than optimizing the Mayer cost, as illustrated in Figure 1(b). Our numerical
results provided in Section 4 confirm that this first approach is inefficient.
A second approach relies in solving directly the full discretization of
Problem ($\mathrm{P}_{\mathrm{Motiv}}$), but with a derivative-free and
blackbox optimization (DFBBO) algorithm. Indeed these algorithms are designed
to solve singular finite-dimensional optimization problems. DFBBO relies only
on a proximal analysis of candidate solutions close to the current incumbent,
and requires neither existence nor values of gradients [6, 15]. In this paper
we focus on the Mads algorithm [5] implemented in the Nomad solver [2], but
other techniques exist (see, e.g., [3, 8, 16, 24]). Mads can be applied to a
full discretization of Problem ($\mathrm{P}_{\mathrm{Motiv}}$), even if the
discontinuities of the original problem are preserved. However Mads is usually
efficient for problems with a few dozens of variables at most [6, Section
1.4], while a discretized optimal control problem usually has a much larger
scale. Our numerical results provided in Section 4 confirm that this second
approach is also inefficient.

###  1.2 Contributions

The difficulties mentioned above are actually present in a whole class of
optimal control problems. For example, any piecewise constant Mayer cost
function leads to a lack of information provided to SO algorithms, because its
gradients are zero (and remain almost-zero after smoothing). As discussed
above, these difficulties cannot be addressed easily with common techniques
from the literature. Therefore, the objective of the present paper is to
propose a numerical method to solve efficiently optimal control problems with
a piecewise constant Mayer cost function, in view of handling various
situations of possible plateaus (as illustrated in Figure 2 in the two-
dimensional case $n=2$). To illustrate polygonal plateaus as in Figure 2(c),
consider parking a car downtown where we look for a cheap parking lot, while
the position of the exact parking spot within this parking lot is less
relevant. In such a problem, the parking lots are the plateaus, each having
its own Mayer cost (for example the price to park in the lot), while all
parking spots in the same parking lot have the same quality. However note that
reaching a parking spot far from the entrance requires more displacement time.
Hence some local considerations (taken into account by the integral Lagrange
cost) may discriminate two parking spots, but this optimization has a
noticeably smaller impact than the choice of the parking.

Figure 2: Examples of possible plateaus in the two-dimensional case $n=2$. In
all figures of the present paper, the light-colored plateaus represent the
lowest-cost plateaus to reach, while the dark ones represent the highest-cost
plateaus to avoid.

To propose a numerical method to solve efficiently optimal control problems
with a piecewise constant Mayer cost function, we exploit the fact that the
inner difficulty depends entirely on the final state variable only. The main
idea is thus to optimize the final state alone via a DFBBO algorithm, while a
smooth subproblem, which consists in solving the initial problem with the
final state being fixed and with a minimal Lagrange cost, is solved via a SO
algorithm. The resulting method is called hybrid as it combines tools from
DFBBO and SO.

The idea to separate variables can be extended to handle cases where the final
state is not the only variable affecting a singularity. Consider for example a
variant of Problem ($\mathrm{P}_{\mathrm{Motiv}}$) where the Mayer cost is
defined as $\left\lfloor\lVert
x(T)-x_{\mathrm{end}}\rVert\right\rfloor+\left\lfloor\lVert
x(T/2)-x_{\mathrm{mid}}\rVert\right\rfloor$, requiring for the trajectory to
approach a target $x_{\mathrm{mid}}$ at the midtime and a target
$x_{\mathrm{end}}$ at the final time. Here the idea is to optimize the states
$x(T/2)$ and $x(T)$ via a DFBBO algorithm and solve via a SO algorithm two
smooth subproblems recovering a trajectory joining respectively $x(0)$ to
$x(T/2)$ and $x(T/2)$ to $x(T)$, with a minimal Lagrange cost.

Hence the main generic problem that we address in this work is given by

$\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}x\in\mathcal{AC}([0,T],\mathbb{R}^{n})\\\
u\in\mathcal{L}^{\infty}([0,T],\mathbb{R}^{m})\end{array}}{\operatorname*{minimize}}&\begin{array}[]{l}g(\\{x(t)\\}_{t\in\Omega})+\displaystyle\int_{0}^{T}\ell(t,x(t),u(t))\,dt\end{array}\\\\[23.68048pt]
\mbox{subject to}&\begin{array}[t]{ll}\dot{x}(t)=f(t,x(t),u(t)),&\mbox{a.e.\
}t\in[0,T],\\\\[3.0pt] x(0)=x_{0},\\\\[3.0pt] c(t,x(t),u(t))\leq
0,&\mbox{a.e.\ }t\in[0,T],\\\\[3.0pt] c_{T}(x(T))\leq 0\end{array}\end{array}$
($\mathrm{P}_{\mathrm{Generic}}$)

where the Mayer cost depends on the state of the control system at one or
several times (indexed by the set $\Omega\subset[0,T]$ with finite
cardinality) and the corresponding Mayer cost function is piecewise constant
in a generic sense, that is, without being restricted to a discrete distance
to some reference targets. Actually, even if we focus on piecewise constant
Mayer cost functions in this paper, our work can be considered for any Mayer
cost function presenting singularities such as discontinuities (when involving
the floor function $\left\lfloor\cdot\right\rfloor$ for example).

The hybrid method proposed in this paper relies on the knowledge of the state
variables which may affect the singularities. Precisely, once Problem
($\mathrm{P}_{\mathrm{Generic}}$) is fully discretized and seen as a finite-
dimensional optimization problem, we reformulate it as a main problem handling
the state variables which may affect the singularities, and involving
subproblems dealing with all remaining state-control variables. The main
problem is solved using a DFBBO algorithm. Each subproblem’s goal is to
recover a partial feasible trajectory-control pair joining with a minimal
partial Lagrange cost two consecutive states proposed by the main DFBBO
algorithm. Since the variables affecting the singularities are fixed, the
subproblems are smooth discrete optimal control problems with fixed terminal
states, solved efficiently via a SO algorithm. Hence the hybrid method applies
both DFBBO and SO algorithms,

The goal of this paper is to formally introduce the hybrid method. Section 2
develops a mathematical framework to express precisely the reformulation
discussed above, but for a generic singular finite-dimensional optimization
problem, allowing to discuss properties of the reformulation in a simple and
general context. The specification to discretized optimal control problems is
done in Section 3, with details on how the hybrid method works in that
context. Section 4 proposes numerical experiments showing that the hybrid
method is efficient for optimal control problems with a piecewise constant
Mayer cost function. Finally Section 5 provides an opening discussion on
future perspectives and possible improvements of the hybrid method.

Let us observe that the hybrid method has connections with trajectory
optimization [19, 28] in a context of driving an aircraft. However, in
contrary to the present work, these papers require existence and explicit
knowledge of first- and second-order information about the system and the
costs. The hybrid method also has connections with bilevel optimization [17,
33] in which the common framework defines an upper level trying to minimize
the optimal value of a lower level parameterized by the upper level. However
the bilevel paradigm involves a conflictual leader-follower structure, thus it
differs from this work which simply partitions the variables.

##  2 Reformulation of a generic singular optimization problem

This paper proposes a numerical approach dedicated to solve a class of
discretized optimal control problems. These problems can be written as finite-
dimensional constrained optimization problems of the form

$\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}x\in\mathbb{R}^{n}\end{array}}{\operatorname*{minimize}}&\begin{array}[]{l}\varphi(x)\end{array}\\\\[4.30554pt]
\mbox{subject to}&\begin{array}[t]{ll}c(x)\leq 0\end{array}\end{array}$
($\mathcal{P}$)

with $n\geq 1$ real variables $x\in\mathbb{R}^{n}$, where
$\varphi:\mathbb{R}^{n}\to\mathbb{R}$ is the objective function, and with
$d\geq 1$ inequality constraints expressed as a component-wise vectorial
inequality $c(x)\leq 0$ via the constraint function
$c:\mathbb{R}^{n}\to\mathbb{R}^{d}$. In some cases, the functions $\varphi$ or
$c$ may be singular (e.g., discontinuous) with only a known subset of
variables affecting the singularities. Exploiting this knowledge, this section
focuses on the reformulation of ($\mathcal{P}$) given by

$\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}y\in\mathbb{R}^{p}\end{array}}{\operatorname*{minimize}}&\left(\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}z\in\mathbb{R}^{q}\end{array}}{\inf}&\begin{array}[]{l}\varphi(y,z)\end{array}\\\\[4.30554pt]
\mbox{subject to}&\begin{array}[t]{ll}c(y,z)\leq
0\end{array}\end{array}\right)\end{array}$ ($\mathcal{R}$)

where the $n$ real variables $x\in\mathbb{R}^{n}$ are partitioned into two
disjoint groups, being expressed as $x=(y,z)$ and ordered so that
$y\in\mathbb{R}^{p}$ stands for the $p\in\llbracket 0,n\rrbracket$ variables
which may affect the singularities of ($\mathcal{P}$) (called singular
variables) and $z\in\mathbb{R}^{q}$ stands for the remaining $q=n-p$ variables
(called smooth variables). The content in the large parenthesis of
($\mathcal{R}$), depending on $y\in\mathbb{R}^{p}$, represents the (possibly
extended-) real value $\inf\\{\varphi(y,z):z\in\mathbb{R}^{q},\;c(y,z)\leq
0\\}$. The format includes situations where $p=0$ or $p=n$, but the case of
interest is $p\in\llbracket 1,n-1\rrbracket$ with $n\geq 2$. In view of
analysing Reformulation ($\mathcal{R}$), we introduce for any
$y\in\mathbb{R}^{p}$ the subproblem ($\mathcal{S}(y)$) given by

$\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}z\in\mathbb{R}^{q}\end{array}}{\operatorname*{minimize}}&\begin{array}[]{l}\varphi(y,z)\end{array}\\\\[4.30554pt]
\mbox{subject to}&\begin{array}[t]{ll}c(y,z)\leq 0.\end{array}\end{array}$
($\mathcal{S}(y)$)

Reformulation ($\mathcal{R}$) is introduced to numerically solve Problem
($\mathcal{P}$) when the latter has singularities (e.g., discontinuities)
affected only by a known subset of variables. In that context, solving
($\mathcal{P}$) numerically may be hard since it is a singular optimization
problem in which gradients are not available. This falls into the DFBBO
framework in which neither explicit knowledge nor existence of derivatives is
required. However a DFBBO algorithm is usually less efficient than a SO
algorithm using gradients when those are available [6, Section 1.4]. Therefore
Reformulation ($\mathcal{R}$) is relevant since it separates the variables
affecting the singularities of ($\mathcal{P}$) from the others. Hence, the
main idea of the hybrid method introduced in this paper is to numerically
solve ($\mathcal{R}$) instead of ($\mathcal{P}$), by optimizing the variables
$y$ at the singular unconstrained upper level via a DFBBO algorithm, while the
remaining variables $z$ are optimized by solving ($\mathcal{S}(y)$) via a SO
algorithm. This approach gives as few variables as possible to a DFBBO
algorithm and leaves all the others to a smooth subproblem solved with a more
efficient SO algorithm.

The section is organized as follows. First Section 2.1 recalls standard
notions, such as the definitions of a feasible and a well-defined optimization
problem. Section 2.2 analyzes relationships existing between ($\mathcal{P}$)
and ($\mathcal{R}$). Section 2.3 discusses issues arising in the numerical
resolution of ($\mathcal{R}$) via the hybrid method.

###  2.1 Feasible and well-defined optimization problems

Consider the optimization problem ($\mathcal{P}$). The points
$x\in\mathbb{R}^{n}$ such that $c(x)\leq 0$ form the feasible set. A feasible
point $x^{*}$ is said to be an optimal solution of ($\mathcal{P}$) if it
satisfies $\varphi(x^{*})\leq\varphi(x)$ for all feasible points $x$. In that
case $\varphi(x^{*})$ is called the optimal value of ($\mathcal{P}$).

A constraint violation function for ($\mathcal{P}$) quantifies the feasibility
or infeasibility of any point $x\in\mathbb{R}^{n}$. Writing
$c(x)=(c_{1}(x),\dots,c_{d}(x))$ for all $x\in\mathbb{R}^{n}$, the constraint
violation function $h_{c}$ considered in this work, as in [20], is given by

$h_{c}:\left\\{\begin{array}[]{ccl}\mathbb{R}^{n}&\to&\mathbb{R}\\\
x&\mapsto&\sum\limits_{k=1}^{d}\max\\{c_{k}(x),0\\}^{2}\end{array}\right.$

in which the square allows to avoid introducing additional singularities.

######  Definition 2.1.

For any $0\leq\varepsilon\leq+\infty$, the _$\varepsilon$ -feasible set_ of
($\mathcal{P}$) is defined by

$\Theta_{\varepsilon}=\\{x\in\mathbb{R}^{n}:h_{c}(x)\leq\varepsilon\\}$

and Problem ($\mathcal{P}$) is said to be $\varepsilon$_-feasible_ (denoted by
$\varepsilon$-F) if $\Theta_{\varepsilon}\neq\emptyset$, while it is said to
be $\varepsilon$_-infeasible_ otherwise. When ($\mathcal{P}$) is
$\varepsilon$-F for some $0\leq\varepsilon\leq+\infty$, it is said to be
$\varepsilon$_-feasible-undefined_ if $\varphi(\Theta_{\varepsilon})$ is not
bounded below, and $\varepsilon$_-feasible-defined_ ($\varepsilon$-FD) if
$\varphi(\Theta_{\varepsilon})$ is bounded below, and $\varepsilon$_-feasible-
well-defined_ ($\varepsilon$-FWD) if $\varphi(\Theta_{\varepsilon})$ is
bounded below and attains its lower bound.

In the sequel the $\varepsilon$ will be omitted in the acronyms of Definition
2.1 when it equals $0$ (e.g., a FWD problem is a $0$-FWD problem). Note that
($\mathcal{P}$) has an optimal solution if and only if ($\mathcal{P}$) is FWD.

######  Definition 2.2.

When it exists, the _minimal infeasibility_ of ($\mathcal{P}$) is defined as
the smallest value $0\leq\varepsilon\leq+\infty$ for which the set
$\Theta_{\varepsilon}$ is nonempty.

Note that ($\mathcal{P}$) has a minimal infeasibility if and only if the
optimization problem

$\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}x\in\mathbb{R}^{n}\end{array}}{\operatorname*{minimize}}&h_{c}(x)\end{array}$
($\mathcal{H}$)

is FWD. If so, the optimal value of ($\mathcal{H}$) equals the minimal
infeasibility of ($\mathcal{P}$).

######  Remark 2.3.

Consider Problem ($\mathcal{P}$). The set $\Theta_{0}$ is the feasible set,
$\Theta_{\infty}=\mathbb{R}^{n}$, and there necessarily exists
$\varepsilon\geq 0$ such that $\Theta_{\varepsilon}\neq\emptyset$. However
($\mathcal{P}$) may be infeasible with no minimal infeasibility. For example,
with $c:x\in\mathbb{R}\mapsto e^{x}\in\mathbb{R}$ and $\varphi\equiv 0$, there
is no smallest value $0\leq\varepsilon\leq+\infty$ for which
$\Theta_{\varepsilon}\neq\emptyset$.

The concepts introduced above for Problem ($\mathcal{P}$) are extended to
Reformulation ($\mathcal{R}$) and to Subproblems ($\mathcal{S}(y)$), replacing
$h_{c}$ by $h_{c(y,\cdot)}$, for all $y\in\mathbb{R}^{p}$.

###  2.2 Relationships between Problem ($\mathcal{P}$) and Reformulation
($\mathcal{R}$)

The equality

$\left(\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}x\in\mathbb{R}^{n}\end{array}}{\inf}&\begin{array}[]{l}\varphi(x)\end{array}\\\\[4.30554pt]
\mbox{subject to}&\begin{array}[t]{ll}c(x)\leq
0\end{array}\end{array}\right)\quad=\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}y\in\mathbb{R}^{p}\end{array}}{\inf}&\left(\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}z\in\mathbb{R}^{q}\end{array}}{\inf}&\begin{array}[]{l}\varphi(y,z)\end{array}\\\\[4.30554pt]
\mbox{subject to}&\begin{array}[t]{ll}c(y,z)\leq
0\end{array}\end{array}\right)\end{array}$

holds true, even when the infima are not attained, or even when they are equal
to $\pm\infty$. Hereafter we investigate when the infima are attained (and
thus are minima) and what relations can be obtained between the minimizers.
The next proposition (whose proof follows directly from Definition 2.1)
establishes an equivalence result on the FWD nature of ($\mathcal{P}$) and the
FWD nature of ($\mathcal{R}$), and, if so, derives some relations between
their respective optimal solutions.

######  Proposition 2.4.

The following statements hold:

  1. 1.

If ($\mathcal{P}$) is FWD, then, for any optimal solution
$x^{*}=(y^{*},z^{*})$ of ($\mathcal{P}$), $y^{*}$ is an optimal solution of
($\mathcal{R}$) (which is FWD) and $z^{*}$ is an optimal solution of
$\mathcal{S}(y^{*})$ (which is FWD).

  2. 2.

If ($\mathcal{R}$) is FWD and has an optimal solution $y^{*}$ for which
$\mathcal{S}(y^{*})$ is FWD, then, for any optimal solution $z^{*}$ of
$\mathcal{S}(y^{*})$, $x^{*}=(y^{*},z^{*})$ is an optimal solution of
($\mathcal{P}$) (which is FWD).

The second item of Proposition 2.4 guarantees that we can recover an optimal
solution of ($\mathcal{P}$) from an optimal solution $y^{*}$ of
($\mathcal{R}$), provided that $\mathcal{S}(y^{*})$ is FWD. The next
counterexamples illustrate this requirement, with two different obstructions.
In these examples (and the others of this section), note that the singularity
affecting the variable $y$ is not mandatory but allows to remain consistent
with our setting.

######  Example 2.5.

Consider $\varphi:(y,z)\in\mathbb{R}^{2}\mapsto
e^{\left\lvert{y}\right\rvert+z}\in\mathbb{R}$ and $c\equiv 0$. In that
context ($\mathcal{R}$) is FWD and any $y\in\mathbb{R}$ is an optimal solution
of ($\mathcal{R}$). However, for any $y\in\mathbb{R}$, ($\mathcal{S}(y)$) is
FD but not FWD. Note that ($\mathcal{P}$) is FD, but is not FWD.

######  Example 2.6.

Consider $\varphi:(y,z)\in\mathbb{R}^{2}\mapsto\left\lvert{y}\right\rvert
e^{z}+\left\lvert{y-1}\right\rvert z^{2}\in\mathbb{R}$ and $c\equiv 0$. In
that context ($\mathcal{P}$) is FWD and its unique optimal solution is
$x^{*}=(0,0)$. Furthermore ($\mathcal{R}$) is FWD and has two optimal
solutions $y^{*}_{1}=0$ and $y^{*}_{2}=1$. Note that $\mathcal{S}(y^{*}_{1})$
is FWD (with $z^{*}_{1}=0$ an optimal solution), while
$\mathcal{S}(y^{*}_{2})$ is FD but not FWD.

Proposition 2.4 provides only a theoretical result establishing relations
between the optimal solution $x^{*}=(y^{*},z^{*})$ of ($\mathcal{P}$), the
optimal solution $y^{*}$ of ($\mathcal{R}$) and the optimal solution $z^{*}$
of $\mathcal{S}(y^{*})$. These relations rely in particular on the FWD nature
of $\mathcal{S}(y^{*})$. Note that, in some cases, ($\mathcal{S}(y)$) may be
FWD if and only if $y=y^{*}$ (see Example 2.7). It is also possible that
($\mathcal{S}(y)$) is infeasible for some $y\in\mathbb{R}^{p}$ arbitrarily
close to $y^{*}$ (see Example 2.8).

######  Example 2.7.

Consider $\varphi:(y,z)\in\mathbb{R}^{2}\mapsto\left\lvert{y}\right\rvert
e^{z}\in\mathbb{R}$ and $c\equiv 0$. In that context ($\mathcal{P}$) has an
optimal solution $x^{*}=(y^{*},z^{*})=(0,0)$. However, when $y\neq y^{*}$,
$\mathcal{S}(y)$ is FD but is not FWD.

######  Example 2.8.

Consider
$\varphi:(y,z)\in\mathbb{R}^{2}\mapsto\left\lvert{y}\right\rvert+z^{2}\in\mathbb{R}$
and $c:(y,z)\in\mathbb{R}^{2}\mapsto(1-y-z,1-y+z)\in\mathbb{R}^{2}$. In that
context ($\mathcal{P}$) has an optimal solution $x^{*}=(y^{*},z^{*})=(1,0)$,
while ($\mathcal{S}(y)$) is infeasible for any $y<1$.

Examples 2.7 and 2.8 show that, when looking for a (numerical) optimal
solution of ($\mathcal{R}$), we may be confronted to subproblems
($\mathcal{S}(y)$) that are FD but not FWD, or even infeasible. These
situations have to be addressed carefully because any attempt to evaluate the
objective value of ($\mathcal{R}$) at some $y\in\mathbb{R}^{p}$ requires,
roughly speaking, a (numerical) optimal value of ($\mathcal{S}(y)$). A FD
subproblem ($\mathcal{S}(y)$) is not hard to handle, as we may obtain a
feasible point $z\in\mathbb{R}^{q}$ with objective value arbitrarily close to
the infimum, but an infeasible subproblem ($\mathcal{S}(y)$) is more
problematic. The next section discusses a dedicated numerical process to face
such an infeasibility. It relies on the next proposition, whose proof is
simple and left to the reader.

######  Proposition 2.9.

If $h_{c(y,\cdot)}$ admits a minimum for all $y\in\mathbb{R}^{p}$, then
$\mathcal{S}(y)$ has a minimal infeasibility, denoted by $\mu_{y}$, for all
$y\in\mathbb{R}^{p}$. If moreover ($\mathcal{P}$) is $\infty$-FD, then
$\mathcal{S}(y)$ is $\mu_{y}$-FD for all $y\in\mathbb{R}^{p}$.

###  2.3 Issues concerning the numerical resolution of Reformulation
($\mathcal{R}$)

In the hybrid method, the objective function of ($\mathcal{R}$) is treated as
a blackbox, that is, as a function with intractable expression which requires
here, roughly speaking, to solve a subproblem ($\mathcal{S}(y)$) depending on
the input $y\in\mathbb{R}^{p}$. This blackbox has to be handled carefully
since two main issues (discussed and illustrated in Section 2.2) need to be
addressed. First ($\mathcal{S}(y)$) may be infeasible and, second, denoting by
$\mu_{y}$ its minimal infeasibility (when it exists), ($\mathcal{S}(y)$) may
be $\mu_{y}$-FD but not $\mu_{y}$-FWD. Assuming that the conditions expressed
in Proposition 2.9 are satisfied, the following two-phases numerical process
is proposed to handle these drawbacks:

  * (i)

Solve with a SO algorithm the smooth unconstrained FWD problem

$\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}z\in\mathbb{R}^{q}\end{array}}{\operatorname*{minimize}}&h_{c(y,\cdot)}(z).\end{array}$
($\mathcal{H}_{y}$)

Denote by $\tilde{z}_{y}$ a (numerical) optimal solution of
($\mathcal{H}_{y}$) and get $\tilde{\mu}_{y}=h_{c(y,\cdot)}(\tilde{z}_{y})$ as
(numerical) approximation of the minimal infeasibility of ($\mathcal{S}(y)$).

  * (ii)

Solve with a SO algorithm the smooth constrained FD problem

$\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}z\in\mathbb{R}^{q}\end{array}}{\operatorname*{minimize}}&\begin{array}[]{l}\varphi(y,z)\end{array}\\\\[4.30554pt]
\mbox{subject
to}&\begin{array}[t]{ll}h_{c(y,\cdot)}(z)\leq\tilde{\mu}_{y}\end{array}\end{array}$
($\tilde{\mathcal{S}}_{y}$)

from the starting point $\tilde{z}_{y}$. Denote by $\tilde{z}_{y}^{*}$ the
(numerical) optimal solution of ($\tilde{\mathcal{S}}_{y}$) and get
$\varphi(y,\tilde{z}_{y}^{*})$ as (numerical) approximation of the optimal
value of ($\tilde{\mathcal{S}}_{y}$). In practice it is not restrictive that
($\tilde{\mathcal{S}}_{y}$) may be not FWD, since it is FD and solved up to
some numerical thresholds. Return $\tilde{\mu}_{y}$ and
$\varphi(y,\tilde{z}_{y}^{*})$.

The goal of the above two-phases process is the following. If
($\mathcal{S}(y)$) is feasible for a given $y\in\mathbb{R}^{p}$, then the
first step leads to a numerically feasible solution of ($\mathcal{S}(y)$), the
second step numerically solves ($\mathcal{S}(y)$) and the process returns a
numerically optimal value of ($\mathcal{S}(y)$). Otherwise, if
($\mathcal{S}(y)$) is infeasible, then an arbitrary choice has to be made to
assign a value to the blackbox at $y$. We could set the value $+\infty$, but
such a choice would not provide any information to pursue efficiently the
numerical resolution of ($\mathcal{R}$). Instead, via the above process, we
recover an infeasibility measure, as well as a corresponding finite objective
value. Recall that DFBBO algorithms can exploit these two (possibly
conflictual) metrics. For example the Mads algorithm [5] implemented in Nomad
[2] handles them via its progressive barrier formulation [4]. It constructs an
approximated Pareto front of undominated infeasible solutions on a diagram
showing objective value versus infeasibility. The front is truncated to
eliminate solutions with an infeasibility above a threshold decreasing with
the iterations. Each iteration of Mads starts from the solution with the
lowest objective function value among the truncated front, so that the
algorithm eventually finds a feasible solution with a low objective value.
Example 2.10 illustrates the relevancy of this approach.

######  Example 2.10.

Consider Example 2.8 and an attempt to numerically solve ($\mathcal{R}$) from
the initial guess $y=0$. Subproblem ($\mathcal{S}(y)$) is infeasible for any
$y<1$ and the above numerical two-phases process returns
$\tilde{\mu}_{y}=2(1-y)^{2}$ and $\varphi(y,\tilde{z}^{*}_{y})=y$. Hence any
value $0<y\leq 1$ has an higher objective function value than $y=0$ but a
lower infeasibility. An algorithm solving ($\mathcal{R}$) can be driven to
$y=1$, leading to a feasible subproblem ($\mathcal{S}(y)$), via an iterative
reduction of the maximal admissible infeasibility. Then it eventually reaches
$x=(y,\tilde{z}^{*}_{y})=(1,0)$ which is the optimal solution of
($\mathcal{P}$).

##  3 Numerical setup to solve Problem ($\mathrm{P}_{\mathrm{Generic}}$)

A full discretization of Problem ($\mathrm{P}_{\mathrm{Generic}}$), using a
Runge-Kutta method for example, can be written as

$\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}\\{x_{k}\\}_{k=1}^{N}\subset\mathbb{R}^{n},\\\
\\{u_{k}\\}_{k=0}^{N-1}\subset\mathbb{R}^{m}\end{array}}{\operatorname*{minimize}}&\begin{array}[]{l}g\left(\\{x_{k}\\}_{k\in\Omega}\right)+\sum\limits_{k=0}^{N-1}\ell_{k}\left(x_{k},u_{k}\right)\end{array}\\\\[23.68048pt]
\mbox{subject
to}&\begin{array}[t]{ll}x_{k+1}=f_{k}(x_{k},u_{k}),&k\in\llbracket
0,N-1\rrbracket,\\\ c_{k}(x_{k},u_{k})\leq 0,&k\in\llbracket
0,N-1\rrbracket,\\\ c_{N}(x_{N})\leq 0\end{array}\end{array}$ ($\mathrm{P}$)

where $n\geq 1$, $m\geq 1$ and $N\geq 1$ are fixed positive integers,
$\\{x_{k}\\}_{k=1}^{N}\subset\mathbb{R}^{n}$ are the state variables and
$\\{u_{k}\\}_{k=0}^{N-1}\subset\mathbb{R}^{m}$ are the control variables,
$x_{0}\in\mathbb{R}^{n}$ is the fixed initial state,
$\Omega\subseteq\llbracket 1,N\rrbracket$ is the set (with $\omega\geq 0$
elements) of the indices $k$ such that the piecewise constant Mayer cost
function $g:(\mathbb{R}^{n})^{\omega}\to\mathbb{R}$ actually depends on
$x_{k}$, each $\ell_{k}:\mathbb{R}^{n}\to\mathbb{R}$ is the Lagrange cost
function at time $k$, each
$f_{k}:\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}^{n}$ is the dynamics at
time $k$, each $c_{k}:\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}^{d_{k}}$
(with $d_{k}\geq 1$) is the mixed state-control constraints function at time
$k$ and $c_{N}:\mathbb{R}^{n}\to\mathbb{R}^{d_{N}}$ (with $d_{N}\geq 1$)
stands for the final state constraints function.

Note that Problem ($\mathrm{P}$) is a specific instance of Problem
($\mathcal{P}$). Therefore, following the approach developed in Section 2,
Section 3.1 proposes several possible reformulations of ($\mathrm{P}$) to
handle its singularities. Then Section 3.2 provides a numerical setup to
introduce the DFBBO-SO hybrid method based on these reformulations, and to
compare it with two non-hybrid approaches (one based on a SO algorithm
exclusively, and one based on a DFBBO algorithm exclusively). Some numerical
experiments, showing the efficiency of the hybrid method for optimal control
problems with a piecewise constant Mayer cost function, are provided in the
next Section 4.

###  3.1 Several reformulations of Problem ($\mathrm{P}$)

Take $(i,j)\in\llbracket 0,N\rrbracket^{2}$ with $j>i$ and
$(x_{i},x_{j})\in(\mathbb{R}^{n})^{2}$. The search of a trajectory linking
$x_{i}$ at time $i$ to $x_{j}$ at time $j$, with minimal partial Lagrange cost
and satisfying the state-control constraints of Problem ($\mathrm{P}$), leads
to the subproblem given by

$\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}\\{x_{k}\\}_{k=i+1}^{j-1}\subset\mathbb{R}^{n},\\\
\\{u_{k}\\}_{k=i\phantom{+1}}^{j-1}\subset\mathbb{R}^{m}\end{array}}{\operatorname*{minimize}}&\begin{array}[]{l}\sum\limits_{k=i}^{j-1}\ell_{k}\left(x_{k},u_{k}\right)\end{array}\\\\[23.68048pt]
\mbox{subject
to}&\begin{array}[t]{ll}x_{k+1}=f_{k}(x_{k},u_{k}),&k\in\llbracket
i,j-1\rrbracket,\\\ c_{k}(x_{k},u_{k})\leq 0,&k\in\llbracket
i,j-1\rrbracket.\end{array}\end{array}$ ($\mathrm{S}(i,x_{i},j,x_{j})$)

This subproblem is not always FWD in general. Hereafter the (possibly
extended-real) lower bound of its objective function over its feasible set is
denoted by $\underline{S}(i,x_{i},j,x_{j})$.

Now, following the approach developed in Section 2 and using the notation
introduced above, Problem ($\mathrm{P}$) has the reformulation

$\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}\\{x_{k}\\}_{k=1}^{N}\subset\mathbb{R}^{n}\end{array}}{\operatorname*{minimize}}&\begin{array}[]{l}g\left(\\{x_{k}\\}_{k\in\Omega}\right)+\sum\limits_{k=0}^{N-1}\underline{S}(k,x_{k},k+1,x_{k+1})\end{array}\\\\[4.30554pt]
\mbox{subject to}&\begin{array}[t]{ll}c_{N}(x_{N})\leq
0\end{array}\end{array}$ ($\mathrm{R}^{\llbracket 1,N\rrbracket}$)

where the variables are the entire trajectory, instead of the entire
trajectory-control pair in Problem ($\mathrm{P}$). For any given trajectory
$\\{x_{k}\\}_{k=1}^{N}$, the search of a control $\\{u_{k}\\}_{k=0}^{N-1}$
driving $x_{0}$ along this trajectory is implicit. The search of each control
$u_{k}$ linking $x_{k}$ to $x_{k+1}$, while minimizing $\ell_{k}(x_{k},u_{k})$
and satisfying the constraint $c_{k}(x_{k},u_{k})\leq 0$, is done in the
smooth Subproblem $(\mathrm{S}(k,x_{k},k+1,x_{k+1}))$. However there is still
an important number of variables (precisely $nN$) left to Reformulation
($\mathrm{R}^{\llbracket 1,N\rrbracket}$).

It is possible to reduce this number of variables. The variant

$\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}x_{N}\in\mathbb{R}^{n}\end{array}}{\operatorname*{minimize}}&\begin{array}[]{l}g\left(x_{N}\right)+\underline{S}\left(0,x_{0},N,x_{N}\right)\end{array}\\\\[4.30554pt]
\mbox{subject to}&\begin{array}[t]{ll}c_{N}(x_{N})\leq
0\end{array}\end{array}$ ($\mathrm{R}^{\\{N\\}}$)

has only the final state $x_{N}$ as a variable. For any given
$x_{N}\in\mathbb{R}^{n}$, the search of a control linking $x_{0}$ to $x_{N}$,
with minimal Lagrange cost and satisfying the state-control constraints of
Problem ($\mathrm{P}$), is left to the unique subproblem
$(\mathrm{S}(0,x_{0},N,x_{N}))$.

When $N$ is even and $\Omega\subseteq\\{N/2,N\\}$, another variant in which
the variables are the final state $x_{N}$ and the midterm state $x_{N/2}$ is
possible. In that case there are two subproblems. The first one links $x_{0}$
to $x_{N/2}$ and the second one links $x_{N/2}$ to $x_{N}$. Denoting by
$\lambda_{0}=0$, $\lambda_{1}=N/2$, $\lambda_{2}=N$, the corresponding
reformulation is given by

$\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}\\{x_{\lambda_{1}},x_{\lambda_{2}}\\}\subset\mathbb{R}^{n}\end{array}}{\operatorname*{minimize}}&\begin{array}[]{l}g\left(\\{x_{k}\\}_{k\in\Omega}\right)+\sum\limits_{i=0}^{1}\underline{S}\left(\lambda_{i},x_{\lambda_{i}},\lambda_{i+1},x_{\lambda_{i+1}}\right)\end{array}\\\\[4.30554pt]
\mbox{subject to}&\begin{array}[t]{ll}c_{N}(x_{\lambda_{2}})\leq
0.\end{array}\end{array}$ ($\mathrm{R}^{\left\\{N/2,N\right\\}}$)

Finally, from a general point of view, it is possible to choose a priori which
states remain variables of the reformulated problem. To this aim, take
$r\in\llbracket 1,N\rrbracket$ and
$\Lambda=\\{\lambda_{1},\dots,\lambda_{r}\\}\subseteq\llbracket 1,N\rrbracket$
the set of all indices of the states given as variables to the reformulated
problem, with increasing values and satisfying $\lambda_{r}=N$ and
$\Omega\subseteq\Lambda$. Taking $\lambda_{r}=N$ is mandatory because the
constraint $c_{N}(x_{N})\leq 0$ cannot appear in the last subproblem, and
taking $\Omega\subseteq\Lambda$ is required because the Mayer cost cannot
appear in any subproblem. Define also for convenience $\lambda_{0}=0$. Finally
the general reformulation of Problem ($\mathrm{P}$) is given by

$\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}\\{x_{\lambda_{i}}\\}_{i=1}^{r}\subset\mathbb{R}^{n}\end{array}}{\operatorname*{minimize}}&\begin{array}[]{l}g\left(\\{x_{k}\\}_{k\in\Omega}\right)+\sum\limits_{i=0}^{r-1}\underline{S}\left(\lambda_{i},x_{\lambda_{i}},\lambda_{i+1},x_{\lambda_{i+1}}\right)\end{array}\\\\[4.30554pt]
\mbox{subject to}&\begin{array}[t]{ll}c_{N}(x_{N})\leq
0\end{array}\end{array}$ ($\mathrm{R}^{\Lambda}$)

with $nr$ variables and admitting $r$ subproblems, each trying to find a
control linking $x_{\lambda_{i}}$ to $x_{\lambda_{i+1}}$ with a minimal
partial Lagrange cost while satisfying the state-control constraints of
Problem ($\mathrm{P}$). Note that Reformulations ($\mathrm{R}^{\llbracket
1,N\rrbracket}$), ($\mathrm{R}^{\\{N\\}}$) and
($\mathrm{R}^{\left\\{N/2,N\right\\}}$) are particular instances of
Reformulation ($\mathrm{R}^{\Lambda}$).

The choices of $r\in\llbracket 1,N\rrbracket$ and
$\Lambda=\\{\lambda_{1},\dots,\lambda_{r}\\}\subseteq\llbracket 1,N\rrbracket$
follow the knowledge of which states $x_{k}$ may affect the singularities of
Problem ($\mathrm{P}$). Furthermore note that the smaller $r$ is, then the
smaller the number of variables left to Reformulation ($\mathrm{R}^{\Lambda}$)
solved with a DFBBO algorithm is. In that context note that ($\mathrm{P}$) and
($\mathrm{R}^{\Lambda}$) are specific instances of ($\mathcal{P}$) and
($\mathcal{R}$) from Section 2. Indeed the $p$ singular variables are the
state variables indexed by $\Lambda$ and, for any input
$y=(x_{\lambda_{i}})_{i=1}^{p}$, Reformulation ($\mathcal{R}$) has a (unique)
Subproblem ($\mathcal{S}(y)$) obtained by aggregating all independent
subproblems
$(\mathrm{S}(\lambda_{i},x_{\lambda_{i}},\lambda_{i+1},x_{\lambda_{i+1}}))$
together.

###  3.2 Numerical setup

In the next Section 4, numerical simulations will be performed on examples of
Problem ($\mathrm{P}$) in which the Mayer cost function involves the floor
function $\left\lfloor\cdot\right\rfloor$, causing discontinuities. This
section describes three numerical methods that will be implemented and
compared. The first one is a SO method, the second one is a DFBBO method and
the third one is the DFBBO-SO hybrid method based on the reformulations
introduced in Section 3.1. Before coming to these methods, we introduce first
a simple initialization used for each method.

#### Uncontrolled trajectory $\mathcal{T}_{\mathrm{uncontrolled}}$

This trajectory results from the induction $x_{k+1}=f_{k}(x_{k},u_{k})$ for
all $k\in\llbracket 0,N-1\rrbracket$, with each control $u_{k}$ set to zero.
It is used as initialization of the three numerical methods introduced below.

#### The $\mathcal{M}_{\mathrm{SO}}$ method

The purpose of the $\mathcal{M}_{\mathrm{SO}}$ method is to solve Problem
($\mathrm{P}$) by using a SO algorithm, via the JuMP modeling language [18]
and the Ipopt solver, despite the discontinuities due to the floor function
$\left\lfloor\cdot\right\rfloor$. Unfortunately, Ipopt cannot deal with
discontinuous problems. To overcome this obstacle, we use a smooth
approximation of the floor function (see Remark 3.1 below) which leads to a
smoothed version of ($\mathrm{P}$) accepted by the code. However the
$\mathcal{M}_{\mathrm{SO}}$ method remains likely inefficient. Indeed, even
after a smoothing, the plateaus resulting from the floor function remain and
lead to partial derivatives close to zero. Hence, the
$\mathcal{M}_{\mathrm{SO}}$ method may fail to detect the significant gains on
the objective value it would obtain by making some variables cross the floor
function plateaus. The presumed ineffectiveness of the
$\mathcal{M}_{\mathrm{SO}}$ method is confirmed by the numerical results in
Section 4.

#### The $\mathcal{M}_{\mathrm{DFBBO}}$ method

The $\mathcal{M}_{\mathrm{DFBBO}}$ method aims to solve Problem ($\mathrm{P}$)
thanks to the DFBBO algorithm Mads via the JuMP interface of the Nomad solver
[26]. However Nomad is usually restricted to small-scale problems [6, Section
1.4] while ($\mathrm{P}$) usually has a large scale. Fortunately a dimension
reduction is possible, as we can replace all state variables $x_{k}$ by the
trajectory generated by the control variables $u_{k}$ and the equality
constraints $x_{k+1}=f_{k}(x_{k},u_{k})$. However, since each control $u_{k}$
remains a variable, the scale of the reduced problem remains too large for
Nomad, in which a critical value of $50$ is indicated as the maximal number of
variables that can be easily handled [2]. Hence, the
$\mathcal{M}_{\mathrm{DFBBO}}$ method is likely inefficient, which is
confirmed by the numerical results obtained in Section 4. Note that, for fair
comparison with the $\mathcal{M}_{\mathrm{SO}}$ method, the
$\mathcal{M}_{\mathrm{DFBBO}}$ method is executed in Section 4 on the smoothed
version of Problem ($\mathrm{P}$), in addition to the above dimension
reduction.

#### The $\mathcal{M}^{\Lambda}_{\mathrm{hybrid}}$ method

The $\mathcal{M}^{\Lambda}_{\mathrm{hybrid}}$ method is the main contribution
of the present paper. Instead of numerically solving Problem ($\mathrm{P}$),
the $\mathcal{M}^{\Lambda}_{\mathrm{hybrid}}$ method focuses on Reformulation
($\mathrm{R}^{\Lambda}$) introduced in Section 3.1 and numerically solves it
via Nomad, while its objective function is evaluated by solving the smooth
subproblems
($\mathrm{S}(\lambda_{i},x_{\lambda_{i}},\lambda_{i+1},x_{\lambda_{i+1}})$)
via Ipopt. This evaluation relies on the numerical two-phases process
presented in Section 2.3, where the infeasibility and the corresponding finite
objective value are returned to Nomad to be used in its implementation of the
progressive barrier variant of the Mads algorithm [4]. Unlike the
$\mathcal{M}_{\mathrm{SO}}$ and $\mathcal{M}_{\mathrm{DFBBO}}$ methods, the
$\mathcal{M}^{\Lambda}_{\mathrm{hybrid}}$ applies each solver to a problem it
is designed for. Indeed Nomad solves a small-scale blackbox optimization
problem with complex internal structure, while Ipopt solves smooth
optimization problems. Note that, for fair comparison with the two methods
above, the $\mathcal{M}^{\Lambda}_{\mathrm{hybrid}}$ method is executed in
Section 4 on the reformulation of a smoothed version of Problem
($\mathrm{P}$), and not on the reformulation of ($\mathrm{P}$). Also, as in
Section 4 the examples have singularities affected by the state variables
$x_{N}$ and $x_{N/2}$ only, the two variants
$\mathcal{M}^{\\{N\\}}_{\mathrm{hybrid}}$ and
$\mathcal{M}^{\\{N/2,N\\}}_{\mathrm{hybrid}}$ of the
$\mathcal{M}^{\Lambda}_{\mathrm{hybrid}}$ method will be tested.

######  Remark 3.1.

As Ipopt cannot handle discontinuities from the floor function
$\left\lfloor\cdot\right\rfloor$, a smooth approximation is introduced as
follows. First, note that the Heaviside function
$\mathds{1}_{\mathbb{R}_{+}}(x)=0$ if $x<0$ and
$\mathds{1}_{\mathbb{R}_{+}}(x)=1$ if $x\geq 0$, can be approximated by
$\mathds{1}_{\mathbb{R}^{+}}(x)\approx exp\left(-e^{-x}\right)$ for any
$x\in\mathbb{R}$. Strengthening the approximation by a factor
$\tau\in\mathbb{R}_{+}^{*}$ and shifting the discontinuity from $0$ to any
value $q\in\mathbb{R}$, we obtain

$\forall(\tau,q)\in\mathbb{R}_{+}^{*}\times\mathbb{R},\quad\forall
x\in\mathbb{R},\quad\mathds{1}_{[q,+\infty[}(x)\approx
exp\left(-e^{-\tau(x-q)}\right).$

Hence a smooth approximation of $\left\lfloor\cdot\right\rfloor$ can be
derived on any interval $[a,b]\subset\mathbb{R}$, with
$(a,b)\in\mathbb{Z}^{2}$, $a<b$, and with any strengthening factor
$\tau\in\mathbb{R}_{+}^{*}$, as follows:

$\left\lfloor\cdot\right\rfloor_{\tau;a,b}:\left\\{\begin{array}[]{ccl}[a,b]&\to&[a,b]\\\
x&\mapsto&a+\sum\limits_{q=a+1}^{b-1}exp\left(-e^{-\tau(x-q)}\right).\end{array}\right.$
($\left\lfloor\cdot\right\rfloor_{\tau;a.b}$)

Note that $\left\lfloor\cdot\right\rfloor_{\tau;a.b}$ is of class ${\cal
C}^{\infty}$ for any admissible couple $(\tau;a,b)$. Its derivative is close
to $0$ at any $x$ far from $\mathbb{Z}$. Also note that
$\left\lfloor\cdot\right\rfloor_{\tau;a.b}$ is not constant on the intervals
$[q,q+1)$ with $q\in\mathbb{Z}$, and it does not reach integer values when
$x\in\mathbb{Z}$. Figure 3 shows $\left\lfloor\cdot\right\rfloor$, as well as
some approximations $\left\lfloor\cdot\right\rfloor_{\tau;a.b}$ for some
values of $\tau$, $a$ and $b$. The numerical tests presented in Section 4 are
performed with $\tau=75$.

Figure 3: Floor function $\left\lfloor\cdot\right\rfloor$ and smooth
approximations $\left\lfloor\cdot\right\rfloor_{\tau;a.b}$.

##  4 Computational study

The examples studied in this section are of the form of Problem
($\mathrm{P}_{\mathrm{Generic}}$) where the Mayer cost function involves the
floor function $\left\lfloor\cdot\right\rfloor$, causing discontinuities. We
discretize the timelapse $[0,T]$ (with $T>0$) with a constant timestep $T/N$
(with $N\in\mathbb{N}^{*}$), the control system via the standard fourth-order
Runge-Kutta method [31] and the integral Lagrange cost via the left rectangle
rule. This discretization procedure leads to Problem ($\mathrm{P}$) which
preserves the singularities due to the floor function. We implement and
compare the numerical methods introduced in Section 3.2 in the Julia language
using the modeling tool JuMP [18]. Both Ipopt and Nomad [26] are interfaced
with JuMP. As explained in Section 3.2, for fair comparisons between the
numerical methods, Problem ($\mathrm{P}$) is preconditioned using the
smoothing procedure described in Remark 3.1.

This section is organized as follows. Section 4.1 deals with a controlled
harmonic oscillator, where the Mayer cost involves the final kinetic energy
minus a truncation of the final potential energy. Section 4.2 discusses a
Zermelo-type problem representing a boat on a river with a stream, trying to
reach a relevant parking in an harbour. Section 4.3 addresses a similar
Zermelo-type problem but with two targets to reach, one at the midtime $T/2$
and one at the final time $T$. Finally Section 4.4 deals with a truncated
version of the Lotka-Volterra stabilization problem and shows that the best
choice of $\Lambda$ in the $\mathcal{M}^{\Lambda}_{\mathrm{hybrid}}$ method
may be nontrivial.

#### Notation

Hereafter $\mathcal{AC}([0,T],\mathbb{R}^{n})$ denotes the space of absolutely
continuous functions defined on $[0,T]$ with values in $\mathbb{R}^{n}$ and
$\mathcal{L}^{\infty}([0,T],\mathbb{R}^{m})$ denotes the usual Lebesgue space
of essentially bounded control functions defined almost everywhere on $[0,T]$
with values in $\mathbb{R}^{m}$. The acronym a.e. stands for for almost every.
In view of using different weighted Euclidean norms in Sections 4.3 and 4.4,
the notation $\lVert X\rVert_{\alpha}$ stands for the Euclidean norm of the
vector $(\alpha_{1}X_{1},\dots,\alpha_{n}X_{n})$ for all $X$,
$\alpha\in\mathbb{R}^{n}$.

###  4.1 Truncated variant of the harmonic oscillator problem

Consider

$\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}(x,v)\in\mathcal{AC}([0,T],\mathbb{R}^{2})\\\
u\in\mathcal{L}^{\infty}([0,T],\mathbb{R})\end{array}}{\operatorname*{minimize}}&\begin{array}[]{l}\displaystyle\frac{M}{2}v(T)^{2}-\frac{K}{2}\lfloor
x(T)\rfloor^{2}+\frac{1}{T}\int_{0}^{T}u(t)\,dt\end{array}\\\\[23.68048pt]
\mbox{subject to}&\begin{array}[t]{ll}\dot{x}(t)=v(t),&\mbox{a.e.\
}t\in[0,T],\\\\[3.0pt] \dot{v}(t)=\frac{1}{M}(-Kx(t)+u(t)),&\mbox{a.e.\
}t\in[0,T],\\\\[3.0pt] x(0)=1,\\\\[3.0pt] v(0)=0,\\\\[3.0pt]
u(t)\in[0,1],&\mbox{a.e.\ }t\in[0,T]\end{array}\end{array}$ ($\mathcal{HO}$)

where the control system mimics a spring-mass system
$M\ddot{x}(t)+Kx(t)=u(t)$, $x(t)$ and $v(t)$ denote the position and velocity
of a point of mass $M>0$ attached to a spring of constant $K>0$, and where
$u(t)\in[0,1]$ denotes a constrained external force. Here the objective is to
minimize the (smooth) final kinetic energy $\frac{M}{2}v(T)^{2}$, plus the
average control value $\frac{1}{T}\int_{0}^{T}u(t)\,dt$, minus the (singular)
final truncated potential energy $\frac{K}{2}\lfloor x(T)\rfloor^{2}$. In this
example, note that the Mayer cost function is not piecewise constant, but
involves the floor function $\left\lfloor\cdot\right\rfloor$, causing
discontinuities.

Figure 4 reports the numerical results obtained with the methods introduced in
Section 3.2 on Problem ($\mathcal{HO}$), with the parameters $T=60$, $N=600$,
$M=2$ and $K=1/2$. They follow our expectations from Section 3.2. Among the
listed methods, only the $\mathcal{M}^{\Lambda}_{\mathrm{hybrid}}$ method
performs well. The $\mathcal{M}_{\mathrm{SO}}$ method fails to detect the gain
resulting from a large final spring elongation, and the
$\mathcal{M}_{\mathrm{DFBBO}}$ method fails because Nomad faces too many
variables. Both methods $\mathcal{M}^{\\{N\\}}_{\mathrm{hybrid}}$ and
$\mathcal{M}^{\\{N/2,N\\}}_{\mathrm{hybrid}}$ succeed because they use Nomad
with few variables to detect a relevant final state and Ipopt to recover a
trajectory joining it.

Figure 4: Numerical results obtained on Problem ($\mathcal{HO}$). The
$\mathcal{M}_{\mathrm{SO}}$ and $\mathcal{M}_{\mathrm{DFBBO}}$ methods are
both represented on the left graph and remain close to their initialization.
Both are outperformed by the $\mathcal{M}^{\\{N\\}}_{\mathrm{hybrid}}$ and
$\mathcal{M}^{\\{N/2,N\\}}_{\mathrm{hybrid}}$ methods.

###  4.2 A first Zermelo-type problem

In this section we consider a variant of the Zermelo problem, with a piecewise
constant Mayer cost function, given by

$\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}(x,y)\in\mathcal{AC}([0,T],\mathbb{R}^{2})\\\
(u,\theta)\in\mathcal{L}^{\infty}([0,T],\mathbb{R}^{2})\end{array}}{\operatorname*{minimize}}&\begin{array}[]{l}\displaystyle-\left\lfloor\frac{x(T)}{10}\right\rfloor+\frac{1}{T}\int_{0}^{T}u(t)\,dt\end{array}\\\\[23.68048pt]
\mbox{subject
to}&\begin{array}[t]{ll}\dot{x}(t)=s_{1}(x(t),y(t))+u(t)\cos\theta(t),&\mbox{a.e.\
}t\in[0,T],\\\\[3.0pt]
\dot{y}(t)=s_{2}(x(t),y(t))+u(t)\sin\theta(t),&\mbox{a.e.\
}t\in[0,T],\\\\[3.0pt] x(0)=0,\\\\[3.0pt] y(0)=0,\\\\[3.0pt] 0\leq y(t)\leq
y_{\rm{max}},&\forall t\in[0,T],\\\\[3.0pt] y(T)=y_{\rm{max}},\\\\[3.0pt]
u(t)\in[0,1],&\mbox{a.e.\ }t\in[0,T]\end{array}\end{array}$
($\mathcal{Z}_{1}$)

where $(x(t),y(t))$ is the position of a boat on a river
$\mathbb{R}\times[0,y_{\rm{max}}]$ of width $y_{\rm{max}}>0$, navigating from
$(0,0)$ and trying to reach the upper shore
$\mathbb{R}\times\\{y_{\rm{max}}\\}$. The boat is equipped with a motor
delivering a driving force $u(t)(\cos\theta(t),\sin\theta(t))$ with
controllable norm $u(t)\in[0,1]$ and direction $\theta(t)\in\mathbb{R}$. The
river has a stream $s=(s_{1},s_{2})$ (defined in Figure 5). Note that the
stream and the driving force actually do not represent physical forces. Indeed
the stream is a vector field of velocities and the control $(u(t),\theta(t))$
induces an additive term to the velocity of the boat. There are parking lots
of width $10$ in the upper shore. The objective is to reach the rightmost
parking lot at the final time, by taking into account the average value of
$u(t)$.

Figure 5: The stream
$s=(s_{1},s_{2}):\mathbb{R}\times[0,y_{\rm{max}}]\to\mathbb{R}^{2}$ with
$y_{\rm{max}}=6$. It repels the boat from the shores. The uncontrolled
trajectory (with the motor turned off, that is, with $(u,\theta)\equiv(0,0)$)
goes to $x\to-\infty$, $y\to\frac{y_{\mathrm{max}}}{2}$. The horizontal
component $s_{1}$ of the stream is powerful in the middle of the river and
cannot be overcame by the driving force of the boat delivered by the motor.

Figure 6 reports the numerical results obtained with the methods introduced in
Section 3.2 on Problem ($\mathcal{Z}_{1}$), with the parameters $T=40$,
$N=400$ and $y_{\rm{max}}=6$. These results follow our expectations from
Section 3.2. Recall that all methods start from the uncontrolled trajectory
which is infeasible since $y_{N}\neq y_{\rm{max}}$ (see Figure 5). The
$\mathcal{M}_{\mathrm{SO}}$ first generates a feasible trajectory (with
$y_{N}=6$) but with a high objective value ($x_{N}$ is far on the negative
values). Then, it struggles to optimize it because the Mayer cost is constant
on each interval $[10(q-1),10q)$ with $q\in\mathbb{Z}$ (up to slight
variations due to the smoothing, see Remark 3.1). This causes the
$\mathcal{M}_{\mathrm{SO}}$ method to converge to a local optimum. The
$\mathcal{M}_{\mathrm{DFBBO}}$ method fails to find a feasible solution. Both
the $\mathcal{M}^{\\{N\\}}_{\mathrm{hybrid}}$ and
$\mathcal{M}^{\\{N/2,N\\}}_{\mathrm{hybrid}}$ methods recover a feasible
trajectory with a low objective function value. The
$\mathcal{M}^{\\{N\\}}_{\mathrm{hybrid}}$ and
$\mathcal{M}^{\\{N/2,N\\}}_{\mathrm{hybrid}}$ methods both outperform the
other methods by a significant amount.

Figure 6: Numerical results obtained on Problem ($\mathcal{Z}_{1}$). The
$\mathcal{M}^{\\{N\\}}_{\mathrm{hybrid}}$ and
$\mathcal{M}^{\\{N/2,N\\}}_{\mathrm{hybrid}}$ methods perform well. The
$\mathcal{M}_{\mathrm{SO}}$ method successfully reaches the feasible set but
remains in a locally optimal solution, and the $\mathcal{M}_{\mathrm{DFBBO}}$
method fails to find a feasible trajectory.

###  4.3 A second Zermelo-type problem

In this section our objective is to highlight the expected superiority of the
$\mathcal{M}^{\\{N/2,N\\}}_{\mathrm{hybrid}}$ method over the
$\mathcal{M}^{\\{N\\}}_{\mathrm{hybrid}}$ method when considering a problem
with two targets $A$ and $B$ to approach (with respect to a discrete
distance), one at the final time $T$ and one at the midtime $T/2$. To this aim
we consider a second variant of the Zermelo problem given by

$\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}X\in\mathcal{AC}([0,T],\mathbb{R}^{2})\\\
(u,\theta)\in\mathcal{L}^{\infty}([0,T],\mathbb{R}^{2})\end{array}}{\operatorname*{minimize}}&\begin{array}[]{l}\left\lfloor\lVert
X(T)-A\rVert_{\alpha}\right\rfloor+\left\lfloor\lVert
X(\frac{T}{2})-B\rVert_{\beta}\right\rfloor+\displaystyle\dfrac{1}{T}\int_{0}^{T}u(t)\,dt\end{array}\\\\[23.68048pt]
\mbox{subject
to}&\begin{array}[t]{ll}\dot{x}(t)=s_{1}(X(t))+u(t)\cos\theta(t),&\mbox{a.e.\
}t\in[0,T],\\\\[3.0pt] \dot{y}(t)=s_{2}(X(t))+u(t)\sin\theta(t),&\mbox{a.e.\
}t\in[0,T],\\\\[3.0pt] X(0)=(0,0),\\\\[3.0pt]
X(t)\in\mathbb{R}\times[0,y_{\rm{max}}],&\forall t\in[0,T],\\\\[3.0pt]
u(t)\in[0,1],&\mbox{a.e.\ }t\in[0,T]\end{array}\end{array}$
($\mathcal{Z}_{2}$)

with $X=(x,y)$, and using the control system and parameters from Section 4.2
and where $A$, $\alpha$, $B$, $\beta\in\mathbb{R}^{2}$.

Note that the $\mathcal{M}^{\\{N\\}}_{\mathrm{hybrid}}$ method is
theoretically unsuitable for this problem since, with notations from Section
3, we have $\Omega=\\{N/2,N\\}$ and $\Lambda=\\{N\\}$ while the inclusion
$\Omega\subseteq\Lambda$ is required. However, to illustrate the lack of
performance of the $\mathcal{M}^{\\{N\\}}_{\mathrm{hybrid}}$ method on this
problem, we will implement it by adding the term $\left\lfloor\lVert
X(T/2)-B\rVert_{\beta}\right\rfloor$ to the objective of Subproblem
($S(0,X_{0},N,X_{N})$) which is solved by Ipopt.

Figure 7 reports the numerical results obtained with the methods introduced in
Section 3.2 on Problem ($\mathcal{Z}_{2}$), with the parameters $T$, $N$,
$y_{\mathrm{max}}$ and $s$ as in Section 4.2 and with $A=(-325,0)$,
$\alpha=(\frac{1}{20},\frac{5}{6})$, $B=(-175,6)$ and
$\beta=(\frac{1}{20},\frac{1}{2})$. The plateaus associated with $B$ are
represented on the lower graph on the locus $]-225,0]\times[0,6]$, while the
plateaus associated with $A$ appear on the locus $[-350,-225{[}\times[0,6]$.

The $\mathcal{M}_{\mathrm{SO}}$ and $\mathcal{M}_{\mathrm{DFBBO}}$ methods
leave the boat uncontrolled and miss the targets. The
$\mathcal{M}^{\\{N\\}}_{\mathrm{hybrid}}$ and
$\mathcal{M}^{\\{N/2,N\\}}_{\mathrm{hybrid}}$ methods perform better, although
only the $\mathcal{M}^{\\{N/2,N\\}}_{\mathrm{hybrid}}$ method approaches the
midterm target. As expected, only the method explicitly handling both
singularities is reliable. This example shows a case where the dynamics and
the Lagrange cost repel the system from a target, since the stream repels the
boat from the shores while confronting it is penalized by the Lagrange cost.
Then, confronting the stream appears irrelevant for each method except the
$\mathcal{M}^{\\{N/2,N\\}}_{\mathrm{hybrid}}$ method.

Figure 7: Numerical results obtained on Problem ($\mathcal{Z}_{2}$). The
$\mathcal{M}_{\mathrm{SO}}$ and $\mathcal{M}_{\mathrm{DFBBO}}$ methods fail to
differ from the uncontrolled trajectory. The
$\mathcal{M}^{\\{N\\}}_{\mathrm{hybrid}}$ method reaches the lowest-cost
plateau at the final time but fails to approach the midterm reference state.
The $\mathcal{M}^{\\{N/2,N\\}}_{\mathrm{hybrid}}$ method reaches both lowest-
cost plateaus.

###  4.4 A truncated variant of the Lotka-Volterra stabilization problem

In the three previous problems, the analysis of the numerical methods follows
our expectations from Section 3.2. In particular the
$\mathcal{M}^{\Lambda}_{\mathrm{hybrid}}$ method performs better than the
others, and $\Lambda$ must be chosen accordingly to the singularities
appearing in the Mayer cost function. Precisely $\Lambda$ must contain the set
$\Omega$ and, at a first glance, there is no need to make it larger. However,
and surprisingly, some singular problems with $\Omega=\\{N\\}$ may make the
$\mathcal{M}^{\\{N/2,N\\}}_{\mathrm{hybrid}}$ method perform better than the
$\mathcal{M}^{\\{N\\}}_{\mathrm{hybrid}}$ method. In this section we enlighten
and analyze this phenomenon on an example, based on a simplified two-
dimensional Lotka-Volterra biological system, given by

$\begin{array}[]{cl}\underset{\scriptsize\begin{array}[]{c}X=(x,y)\in\mathcal{AC}([0,T],\mathbb{R}^{2})\\\
u\in\mathcal{L}^{\infty}([0,T],\mathbb{R})\end{array}}{\operatorname*{minimize}}&\begin{array}[]{l}\displaystyle\left\lfloor\lVert
X(T)-\overline{X}\rVert_{\alpha}\right\rfloor+\frac{1}{T}\int_{0}^{T}u(t)\,dt\end{array}\\\\[23.68048pt]
\mbox{subject
to}&\begin{array}[t]{ll}\dot{x}(t)=x(t)(y(t)-\overline{y}-u(t)),&\mbox{a.e.\
}t\in[0,T],\\\\[3.0pt] \dot{y}(t)=y(t)(\overline{x}-x(t)),&\mbox{a.e.\
}t\in[0,T],\\\\[3.0pt] x(0)=x_{0},\\\\[3.0pt] y(0)=y_{0},\\\\[3.0pt]
u(t)\in[0,\frac{1}{4}],&\mbox{a.e.\ }t\in[0,T]\end{array}\end{array}$
($\mathcal{LV}$)

where $X(t)=(x(t),y(t))$ represents the number of thousands of individuals in
a specie of predators and preys respectively,
$X_{0}=(x_{0},y_{0})\in(\mathbb{R}_{+}^{*})^{2}$ stands for the initial
condition,
$\overline{X}=(\overline{x},\overline{y})\in(\mathbb{R}_{+}^{*})^{2}$ denotes
the equilibrium of the Lotka-Volterra system, and $u(t)\in[0,1/4]$ is the
ratio of predators decimated by an human intervention at time $t$. With no
intervention (that is, with $u\equiv 0$), the Lotka-Volterra system has a
periodical motion. The objective of Problem ($\mathcal{LV}$) is to stabilize
the system by approaching the equilibrium $\overline{X}$ (with respect to a
discrete distance), and minimizing the average human intervention. We use the
parameters $T=60$, $N=600$, $X_{0}=(x_{0},y_{0})=(2.25,4.25)$ and
$\overline{X}=(\overline{x},\overline{y})=(1,1)$. They lead to the final state
of the uncontrolled trajectory $X_{N}\approx(0.05,1.59)$.

The controls and trajectories found by each method developed in Section 3.2
are shown on Figure 8. The $\mathcal{M}_{\mathrm{SO}}$ and
$\mathcal{M}_{\mathrm{DFBBO}}$ methods fail to differ from the uncontrolled
trajectory. Contrary to the prior problems, the
$\mathcal{M}^{\\{N\\}}_{\mathrm{hybrid}}$ method also fails. The
$\mathcal{M}^{\\{N\\}}_{\mathrm{hybrid}}$ method only pushes one control value
at $1/4$ while all the others remain at $0$, to slightly alter the final state
$X_{N}$ and make it cross one single plateau. This behavior is intriguing.
Since this method leaves only the variable $X_{N}$ to Nomad, at a first
glance, there is no reason to explain why the method fails to propose a final
state $X_{N}$ closer to $\overline{X}$. This issue is of prime importance for
further research works and is discussed in more details in the next paragraph.
Finally note that the $\mathcal{M}^{\\{N/2,N\\}}_{\mathrm{hybrid}}$ method
performs well since it recovers a trajectory reaching the lowest-cost plateau.

Figure 8: Numerical results obtained on Problem ($\mathcal{LV}$). The
$\mathcal{M}_{\mathrm{SO}}$ and $\mathcal{M}_{\mathrm{DFBBO}}$ methods are
both represented on the left graph and remain at their initialization. The
$\mathcal{M}^{\\{N\\}}_{\mathrm{hybrid}}$ method also remains on a trajectory
close to the uncontrolled one. The
$\mathcal{M}^{\\{N/2,N\\}}_{\mathrm{hybrid}}$ method performs well and
successfully reaches the lowest-cost plateau.

As suggested by the numerical results obtained, Problem
($\mathrm{R}^{\\{N\\}}$) is difficult to solve numerically in the context of
this section. Recall that, in Problem ($\mathrm{R}^{\\{N\\}}$), the only
variable is the final state $X_{N}$, while the search of a trajectory linking
$X_{0}$ to $X_{N}$ is performed via Subproblem
$(\mathrm{S}(0,X_{0},N,X_{N}))$. Recall also that Ipopt looks for a (feasible
and) optimal solution of the subproblem, but it has no guarantee to localize
neither an optimal solution nor a feasible one. As it aims primarily a
feasible solution before optimizing it, Ipopt may converge to an infeasible
solution in its attempt to solve a feasible subproblem, because this solution
is a local minimum of the infeasibility measure. Hence the variable $X_{N}$
may trigger an hidden constraint in Problem ($\mathrm{R}^{\\{N\\}}$), as it
may generate a subproblem for which Ipopt fails to find a feasible solution.
This constraint does not appear on the mathematical formulation of Problem
($\mathcal{LV}$), as it comes from the computer code implementing the
$\mathcal{M}^{\Lambda}_{\mathrm{hybrid}}$ method, but it has an impact on the
returned solution. This phenomenon explains why Nomad struggles to solve
Problem ($\mathrm{R}^{\\{N\\}}$) in the context of this section. Indeed,
consider Figure 9 which shows the infeasibility of the solution returned by
Ipopt on Subproblem ($\mathrm{S}(0,X_{0},N,X_{N})$) for each point $X_{N}$.
The set of $X_{N}$ leading to a subproblem for which the solution returned by
Ipopt has a high infeasibility measure surrounds the equilibrium
$\overline{X}$. Therefore, starting from the final state of the uncontrolled
trajectory $X_{N}\approx(0.05,1.59)$ in the large black ribbon of Figure 9,
Nomad sees its initialization as a feasible solution but also sees the colored
zone as a zone of high infeasibility, and thus remains in the black ribbon.
Finally, as the $\mathcal{M}^{\\{N/2,N\\}}_{\mathrm{hybrid}}$ method performs
better than the $\mathcal{M}^{\\{N\\}}_{\mathrm{hybrid}}$ method, we suspect
that Subproblem $(\mathrm{S}(0,X_{0},N,X_{N}))$ recovering the entire
trajectory is harder to solve numerically by Ipopt than the two Subproblems
$(\mathrm{S}(0,X_{0},N/2,X_{N/2}))$ and $(\mathrm{S}(N/2,X_{N/2},N,X_{N}))$
recovering half of the trajectory each. Furthermore we even found that Ipopt
fails to find a feasible solution of the feasible Subproblem
$(\mathrm{S}(0,X_{0},N,X_{N}))$ taking as $X_{N}$ the final state of the
trajectory returned by the $\mathcal{M}^{\\{N/2,N\\}}_{\mathrm{hybrid}}$
method. This subproblem is feasible (since the trajectory returned by the
$\mathcal{M}^{\\{N/2,N\\}}_{\mathrm{hybrid}}$ method is a feasible solution),
but Ipopt fails to detect a feasible solution.

Figure 9: The infeasibility of the solution returned by Ipopt attempting to
solve Subproblem $(\mathrm{S}(0,X_{0},N,X_{N}))$ for each
$X_{N}=(x,y)\in[0,5]^{2}$. It may differ from the true minimal infeasibility
because Ipopt may converge to a locally infeasible solution in its attempt to
solve a feasible Subproblem $(\mathrm{S}(0,X_{0},N,X_{N}))$.

##  5 Discussion

This work introduces an hybrid numerical method, relying on two algorithms
from DFBBO and SO, to solve optimal control problems with a piecewise constant
Mayer cost function. This hybrid method is based on a reformulation of the
discretized problem which splits the search of a full trajectory-control pair
into several subproblems searching for partial trajectory-control pairs. The
state variables that may impact the singularities of the Mayer cost function
are optimized on the main problem via a DFBBO algorithm, while the search of a
trajectory joining those is left to smooth subproblems solved via a SO
algorithm.

As evidenced in Section 4, the $\mathcal{M}^{\Lambda}_{\mathrm{hybrid}}$
method solves efficiently optimal control problems with a piecewise constant
Mayer cost function, while non-hybrid approaches fail. Although we voluntarily
restricted our work to optimal control problems with a piecewise constant
Mayer cost function only, our method can be extended to a larger class of
optimal control problems, including for example nondifferentiable cost
functions, piecewise continuous cost functions, extended-real valued cost
functions, etc., or cases where the control variables may also impact the
singularities. From a general point of view, any cost or constraint function,
or any dynamics, involving a singularity handled by DFBBO algorithms may be
handled by an adapted hybrid method. However, this last class of problems is
large and some other methods may be more appropriate for some of its
subclasses. A future development could rely in a better identification of the
class of problems on which hybrid methods work well.

As explained above, the $\mathcal{M}^{\Lambda}_{\mathrm{hybrid}}$ method can
be extended in several directions. It can also be improved in several aspects.
First, we know from Section 4.4 that a relevant choice of the subset $\Lambda$
does not simply consist in making it as small as possible, even if a subset
$\Lambda$ with small cardinality is preferred. Hence, a possible improvement
is to fix the cardinal of $\Lambda$ and to leave its elements as variables of
the main problem. Note that such a process would not require any algorithm
from combinatorial optimization, since DFBBO algorithms handle problems with
both discrete and continuous variables [7]. Preliminary numerical experiments
are promising and this extension is expected to be the topic of future
research works.

Second, we noticed that the feasibility of the subproblems highly depends on
the values of the states optimized on the main problem, to the point that it
may be difficult to find states leading to feasible subproblems with a DFBBO
algorithm (see Figure 9 for an example). A possible improvement would be to
relax the constraint forcing the subproblems to match exactly the states
proposed by the main problem. The trajectory-finding subproblems may be
relaxed to simply approach the states instead of matching them. This could be
done by adding a penalization term to the objective functions of the
subproblems, penalizing the distance between the states of the system at the
times indexed by $\Lambda$ and the proposed states at these times.

Third, the known problem structure may be exploited to guide the main DFBBO
algorithm. For example, because of the dynamics of the control system
constraining the trajectory, forcing the system to reach a given state at a
given time restricts the space of admissible states at future times. Hence, in
the $\mathcal{M}^{\Lambda}_{\mathrm{hybrid}}$ method, there are links between
the state variables optimized at the main level. Furthermore, we may construct
surrogates [11, 15, 16, 35] of the subproblems to obtain a fast approximation
of their optimal solutions. For example, we may launch a SO algorithm and stop
its execution after a few iterations, and consider the low-precision solution.
The DFBBO algorithm may use this procedure to estimate how relevant a
potential solution is, with no need to systematically solve the subproblems
with full accuracy. Such a development could significantly reduce the
computational time of the $\mathcal{M}^{\Lambda}_{\mathrm{hybrid}}$ method.

## References

  * [1] Alarie, S. and Amaioua, N. and Audet, C. and Le Digabel, S. and Leclaire, L.-A.: Selection of variables in parallel space decomposition for the mesh adaptive direct search algorithm. Les Cahiers du GERAD. G-2018-38 (2018)
  * [2] Alarie, S., Audet, C. and Le Digabel, S. and Rochon-Montplaisir, V. and Tribes, C.: NOMAD version 4: Nonlinear optimization with the MADS algorithm. Les Cahiers du GERAD. G-2021-23 (2021)
  * [3] Audet, C. and Conn, A.R. and Le Digabel, S. and Peyrega, M.: A progressive barrier derivative-free trust-region algorithm for constrained optimization. Computational Optimization and Applications. Vol. 71, pp. 307–329 (2018)
  * [4] Audet, C. and Dennis Jr., J.E.: A Progressive Barrier for Derivative-Free Nonlinear Programming. SIAM Journal on Optimization. Vol. 20, pp. 445-472 (2009)
  * [5] Audet, C. and Dennis Jr., J.E.: Mesh Adaptive Direct Search Algorithms for Constrained Optimization. SIAM Journal on Optimization. Vol. 17, pp. 188-217 (2006)
  * [6] Audet, C. and Hare, W.: Derivative-Free and Blackbox Optimization, Springer Series in Operations Research and Financial Engineering (2017)
  * [7] Audet, C. and Le Digabel, S. and Tribes, C.: The Mesh Adaptive Direct Search Algorithm for Granular and Discrete Variables. SIAM Journal on Optimization. Vol. 29, pp. 1164–1189 (2019)
  * [8] Bagirov, A.M. and Gaudioso, M. and Karmitsa, N. and Mäkelä, M.M. and Taheri, S.: Numerical Nonsmooth Optimization: State of the Art Algorithms, Springer International Publishing (2020)
  * [9] Bayen, T. and Pfeiffer, L.: Second order conditions for a control problem with discontinuous cost. In: 2019 IEEE 58th Conference on Decision and Control (CDC). pp. 1480–1485 (2019)
  * [10] Bonnard, B. and Faubourg, L. and Trélat, E.: Mécanique céleste et contrôle des véhicules spatiaux, Springer-Verlag Berlin Heidelberg (2006)
  * [11] Booker, A.J. and Dennis Jr., J.E. and Frank, P.D. and Serfaini, D.B. and Torczon, V. and Trosset, M.W.: A Rigorous Framework for Optimization of Expensive Functions by Surrogates. Structural and Multidisciplinary Optimization. Vol. 17, pp. 1–13 (1999)
  * [12] Bressan, A. and Piccoli, B.: Introduction to the Mathematical Theory of Control, American Institute of Mathematical Sciences (2007)
  * [13] Bryson, A. E. and Ho, Y.-C.: Applied Optimal Control: optimization, estimation and control, Taylor & Francis (1975)
  * [14] Clarke, F.: Functional Analysis, Calculus of Variations and Optimal Control, Springer-Verlag London (2013)
  * [15] Conn, A.R. and scheinberg, K. and Vicente, L.: Introduction to Derivative-Free Optimization, Society for Industrial and Applied Mathematics (2009)
  * [16] Custódio, A. and Scheinberg, K. and Vicente, L.: Chapter 37: Methodologies and Software for Derivative-Free Optimization. In: Terlaky, T., Anjos, M.F. and Ahmed, S. (eds.): Advances and Trends in Optimization with Engineering Applications, pp. 495–506. Society for Industrial and Applied Mathematics (2017)
  * [17] Dempe, S.: Bilevel Optimization, Springer Optimization and Its Applications, Springer International Publishing (2020)
  * [18] Dunning, I. and Huchette, J. and Lubin, M. : JuMP: A Modeling Language for Mathematical Optimization. SIAM Review Vol. 59, pp. 295–320 (2017)
  * [19] Fisch, F.:Development of a Framework for the Solution of High-Fidelity Trajectory Optimization Problems and Bilevel Optimal Control Problems. PhD thesis, Technische Universität München (2011)
  * [20] Fletcher, R. and Leyffer, S.: Nonlinear Programming Without a Penalty Function. Mathematical Programming A. Vol. 91, pp. 239–269 (2002)
  * [21] Frankowska, H.: The First Order Necessary Conditions for Nonsmooth Variational and Control Problems. SIAM Journal on Control and Optimization. Vol. 22, pp. 1–12 (1984)
  * [22] Haberkorn, T. and Trélat, E.: Convergence Results for Smooth Regularizations of Hybrid Nonlinear Optimal Control Problems. SIAM Journal on Control and Optimization. Vol. 49, pp. 1498–1522 (2011)
  * [23] Levèvre, C.: Optimal Control of a Birth and Death Epidemic Process. Operations Research. Vol. 29, 971–982 (1981)
  * [24] Liuzzi, G. and Lucidi, S. and Rinaldi, F. and Vicente, L.: Trust-Region Methods for the Derivative-Free Optimization of Nonsmooth Black-Box Functions. SIAM Journal on Optimization. Vol. 29, pp. 3012–3035 (2019)
  * [25] Mombaur, K.: Optimal Control for Applications in Medical and Rehabilitation Technology: Challenges and Solutions. In: Hiriart-Urruty, J.-B., Korytowski, A., Maurer, H. and Szymkat, M. (eds.): Advances in Mathematical Modeling, Optimization and Optimal Control, pp. 103–145. Springer International Publishing (2016)
  * [26] Montoison, A. and Pascal, P. and Salomon, L. : NOMAD.jl: A Julia interface for the constrained blackbox solver NOMAD. Published on GitHub https://github.com/bbopt/NOMAD.jl (2020)
  * [27] Nocedal, J., Wright, S.J.: Numerical Optimization, Springer Series in Operations Research and Financial Engineering, Berlin, second edition (2006)
  * [28] Petersen, F.M. and Cornick, D.E. and Bauer, G.L. and Rehder, J.R.: A two-Level trajectory decomposition algorithm featuring optimal intermediate target selection. Journal of Spacecraft and Rockets. Vol. 14, pp. 676–682 (1977)
  * [29] Pontryagin, L. S. and Boltyanskii, V. G. and Gamkrelidze, R. V. and Mishchenko, E. F.: The mathematical theory of optimal processes, Interscience Publishers John Wiley & Sons, Inc. New York-London (1962)
  * [30] Rao, A.: A Survey of Numerical Methods for Optimal Control. Advances in the Astronautical Sciences. Vol. 135, pp. 497–528 (2010)
  * [31] Schwartz, A. and Polak, E.: Consistent Approximations for Optimal Control Problems Based on Runge–Kutta Integration. SIAM Journal on Control and Optimization. Vol. 34, pp. 1235–1269 (1996)
  * [32] Sethi, S., and Thompson, G.: Optimal Control Theory: Applications to Management Science and Economics, Springer New York (2000)
  * [33] Sinha, A. and Malo, P. and Deb, K.: A Review on Bilevel Optimization: From Classical to Evolutionary Approaches and Applications. IEEE Transactions on Evolutionary Computation. Vol. 22, pp. 276–295 (2018)
  * [34] Sharomi, O. and Tufail, M.: Optimal control in epidemiology. Annals of Operations Research. Vol. 251, pp. 55–71 (2017)
  * [35] Talgorn, B. : sgtelib: Surrogate model library for Derivative-Free Optimization. Published on GitHub https://github.com/bbopt/sgtelib (2019)
  * [36] Trélat, E.: Contrôle optimal : théorie et applications, Vuibert supérieur - Mathématiques concrètes (2005)
  * [37] Vinter, R.: Optimal Control, Birkhäuser Boston (2010)
  * [38] Wächter, A. and Biegler, L.T. : On the Implementation of a Primal-Dual Interior Point Filter Line Search Algorithm for Large-Scale Nonlinear Programming. Mathematical Programming. Vol. 106, pp. 25–57 (2006)