# SympOCnet: Solving optimal control problems with applications to high-
dimensional multi-agent path planning problems††thanks: Submitted to the
editors . This work was supported by OSD/AFOSR MURI grant FA9550-20-1-0358.

Tingwei Meng222Division of Applied Mathematics, Brown University, Providence,
RI 02912, USA<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>George_Karniadakis@Brown.edu). 333Tingwei Meng and
Zhen Zhang contributed equally to this work. Zhen Zhang222Division of Applied
Mathematics, Brown University, Providence, RI 02912, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>George_Karniadakis@Brown.edu). 333Tingwei Meng and Zhen Zhang contributed
equally to this work. Jérôme Darbon222Division of Applied Mathematics, Brown
University, Providence, RI 02912, USA<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>George_Karniadakis@Brown.edu).
444Corresponding author. George Em Karniadakis222Division of Applied
Mathematics, Brown University, Providence, RI 02912, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>George_Karniadakis@Brown.edu). 555School of Engineering, Brown University,
Providence, RI 02912, USA.

###### Abstract

Solving high-dimensional optimal control problems in real-time is an important
but challenging problem, with applications to multi-agent path planning
problems, which have drawn increased attention given the growing popularity of
drones in recent years. In this paper, we propose a novel neural network
method called SympOCnet that applies the Symplectic network to solve high-
dimensional optimal control problems with state constraints. We present
several numerical results on path planning problems in two-dimensional and
three-dimensional spaces. Specifically, we demonstrate that our SympOCnet can
solve a problem with more than $500$ dimensions in 1.5 hours on a single GPU,
which shows the effectiveness and efficiency of SympOCnet. The proposed method
is scalable and has the potential to solve truly high-dimensional path
planning problems in real-time.

###### keywords:

deep neural networks, optimal control, path planning, physics-informed
learning

49M99, 68T07

##  1 Introduction

Optimal control methods are widely applied in many practical problems, which
include path planning [17, 83, 43, 25, 75, 59, 65], humanoid robot control
[30, 34, 56, 36, 33, 26], and robot manipulator control [61, 48, 54, 63, 12].
With the increasing impact of drones in everyday tasks as well as specialized
military tasks, the path planning problem of multiple drones has become
increasingly important. However, it is challenging to model this problem using
the optimal control formulation because the dimensionality is usually very
high, causing difficulty in applying traditional numerical methods, such as
the very accurate pseudospectral (PS) method [31, 32, 77, 80, 38], to solve
it. For instance, for a multi-agent path planning problem with $M$ drones, the
motion of each drone has to be described by a state variable in
$\mathbb{R}^{m}$, therefore the dimension of the corresponding optimal control
problem will be $Mm$, which is huge when the number of drones $M$ is large.
Computational complexity of traditional numerical methods may scale
exponentially with respect to the dimension, which makes the high-dimensional
problems intractable. This is an open issue called “curse-of-dimensionality”
(CoD) [7]. Hence, new approaches are required to tackle high-dimensional
optimal control problems efficiently.

In literature, there are some algorithms for solving high-dimensional optimal
control problems (or the corresponding Hamilton-Jacobi PDEs), which include
optimization methods [18, 24, 21, 23, 15, 14, 87, 60, 55], max-plus methods
[1, 2, 29, 35, 39, 67, 66, 68, 69], tensor decomposition techniques [28, 44,
86], sparse grids [10, 37, 53], polynomial approximation [51, 52], model order
reduction [4, 57], optimistic planning [9], dynamic programming and
reinforcement learning [13, 11, 3, 8, 89], as well as methods based on neural
networks [5, 6, 27, 47, 42, 45, 46, 58, 73, 78, 81, 84, 62, 20, 22, 19, 71,
72, 49, 50, 74]. However, path planning problems are in general hard to solve
and cannot be solved directly by applying most of the aforementioned
algorithms. One difficulty comes from the state constraints given by the
complicated environment. Therefore, solving high-dimensional optimal control
problems with complicated state constraints is an important yet challenging
open problem. The difficulty of high dimensionality can be avoided by solving
the sequential path planning problem instead of the simultaneous planning
problem [13, 11, 79]. The main idea of sequential planning is to assign for
different drones different levels of priorities and then sequentially solve
the low-dimensional path planning problem involving each drone according to
its priority. However, sequential planning only provides a feasible solution
to the original simultaneous planning problem, which may not be optimal.
Hence, how to solve high-dimensional simultaneous planning problems with
complicated state constraints is still an unsolved problem.

Recently, neural networks have achieved great success in solving high-
dimensional PDEs and inverse problems and bear some promise in overcoming the
CoD. We can make further progress if we take advantage of the knowledge about
the underlying data generation process and encode such information in the
neural network architectures and the loss functions [49, 88, 40, 70]. As a
specific example, the physics-informed neural network (PINN), which was first
introduced in [76], encodes physical information directly into the loss
function of neural networks using the residuals of ODEs/PDEs and automatic
differentiation. In this work, we leverage PINNs and propose a novel method
using a Symplectic optimal control network called SympOCnet to solve high-
dimensional optimal control problems. SympOCnet is designed based on the
Symplectic network (SympNet), which is a neural network architecture proposed
in [50] and is able to approximate arbitrary symplectic maps. In the original
paper, SympNet is used to solve the inverse parameter identification problem
by acting as a surrogate model for the data generation procedure. In this
paper, we present a novel way to utilize SympNet for solving Hamiltonian ODEs
in a forward manner. By applying the SympOCnet, we encode our knowledge about
optimal control problems and Hamilton-Jacobi theory in the neural network. We
then apply SympOCnet on multi-agent path planning problems and show its
effectiveness and efficiency in solving high-dimensional problems.
Specifically, we can solve a path planning problem involving $256$ drones,
which corresponds to a $512$-dimensional optimal control problem.

The organization of this paper is given as follows. In Section 2, we briefly
summarize some preliminary background which will be used in the paper,
including Hamiltonian systems and symplectic maps in Section 2.1 and SympNet
in Section 2.2. In Section 3, we present our proposed SympOCnet method for the
optimal control problems without state constraints, which serves as a building
block for the full version of the SympOCnet method. The full version is then
introduced in Section 4, which is proposed for solving the optimal control
problems with state constraints. To be specific, we propose four loss
functions, which are compared in the numerical experiments. Our numerical
experiments on multi-agent path planning problems are presented in Section 5.
The first example in Section 5.1 is used to compare the four loss functions.
The second example in Section 5.2 demonstrates the generalizability of our
method when only partial data is available in the training process. Then, in
Sections 5.3 and 5.4, we demonstrate the effectiveness and efficiency of
SympOCnet on high-dimensional problems, with dimensions ranging from $64$ to
$512$. A summary is provided in Section 6.

##  2 Preliminary background

In this paper, we solve optimal control problems with SympOCnet, a neural
network with SympNet architecture and loss function given by the corresponding
Hamiltonian ODEs. This section provides some mathematical materials used in
the remainder of the paper. In Section 2.1, we provide a brief summary about
symplectic maps, Hamiltonian ODEs and their relations. Then, in Section 2.2,
we review the SympNet architectures. For more details regarding the theory and
applications of symplectic methods, we refer the readers to [41].

###  2.1 Hamiltonian systems and symplectic maps

######  Definition 2.1 (Symplectic maps).

Let $U$ be an open set in $\mathbb{R}^{2n}$. A differentiable map
$\phi:U\rightarrow\mathbb{R}^{2n}$ is called symplectic if the Jacobian matrix
$\nabla\phi$ satisfies

(1)
$\nabla\phi^{T}(\bm{z})J\nabla\phi(\bm{z})=J\quad\forall\bm{z}\in\mathbb{R}^{2n},$

where $J$ is a matrix with $2n$ rows and $2n$ columns defined by

(2) $J:=\begin{pmatrix}\bm{0}&I_{n}\\\ -I_{n}&\bm{0}\end{pmatrix},$

and $I_{n}$ denotes the identity matrix with $n$ rows and $n$ columns.

The Hamiltonian ODE system is a dynamical system taking the form

$\dot{\bm{z}}(s)=J\nabla H(\bm{z}(s)),$

where $J$ is the matrix defined in (2) and $H:U\rightarrow\mathbb{R}$ is a
function called Hamiltonian. The Hamiltonian systems and symplectic maps are
highly related to each other. To be specific, the Hamiltonian structure is
preserved under the change of variable using any symplectic map. This result
is stated in the following theorem.

######  Theorem 2.2.

[41, Theorem 2.8 on p. 187] Let $U$ and $V$ be two open sets in
$\mathbb{R}^{2n}$. Let $\phi:U\rightarrow V$ be a change of coordinates such
that $\phi$ and $\phi^{-1}$ are continuously differentiable functions. If
$\phi$ is symplectic, the Hamiltonian ODE system $\dot{\bm{z}}(s)=J\nabla
H(\bm{z}(s))$ can be written in the new variable $\bm{w}=\phi(\bm{z})$ as

(3) $\dot{\bm{w}}(s)=J\nabla\tilde{H}(\bm{w}(s)),$

where the new Hamiltonian $\tilde{H}$ is defined by

(4) $\tilde{H}(\bm{w})=H(\bm{z})=H(\phi^{-1}(\bm{w}))\quad\forall\bm{w}\in V.$

Conversely, if $\phi$ transforms every Hamiltonian system to another
Hamiltonian system by (3) and (4), then $\phi$ is symplectic.

Theorem 2.2 indicates that a symplectic map can transform a Hamiltonian ODE
system to another system which is again Hamiltonian and potentially in a
simpler form. This is the starting point of our proposed method. It is well
known from literature that the optimal trajectory of an optimal control
problem is related to a Hamiltonian ODE system under certain assumptions. Such
systems may be solved more easily if written in the right coordinates. We
propose to learn the corresponding coordinate transformation through a
parameterized family of symplectic maps $\phi_{\theta}$, where $\theta$
represents the unknown parameters to be learned, and solve the Hamiltonian
ODEs with new coordinates. The solution to the original problem can be
obtained by mapping the trajectory back to original phase space through
$\varphi_{\theta}=\phi_{\theta}^{-1}$.

###  2.2 Symplectic networks (SympNets)

SympNet is a neural network architecture proposed in [50] to approximate
symplectic transformations. There are different kinds of SympNet
architectures. In this paper, we use the G-SympNet as the class of our
symplectic transformation $\varphi_{\theta}$. For other architectures, we
refer the readers to [50].

Define a function
$\hat{\sigma}_{K,\bm{a},\bm{b}}\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ by
$\hat{\sigma}_{K,\bm{a},\bm{b}}(\bm{x}):=K^{T}(\bm{a}\odot\sigma(K\bm{x}+\bm{b}))$
for any $\bm{x}\in\mathbb{R}^{n}$, where $\sigma$ is the activation function,
and $\odot$ denotes the componentwise multiplication. Here $\bm{a},\bm{b}$ are
vectors in $\mathbb{R}^{l}$, $K$ is a matrix with $l$ rows and $n$ columns.
Any G-SympNet is an alternating composition of the following two parameterized
functions:

(5) $\begin{split}&\mathcal{G}_{up}\begin{pmatrix}\bm{x}\\\
\bm{p}\end{pmatrix}=\begin{pmatrix}\bm{x}\\\
\bm{p}+\hat{\sigma}_{K,\bm{a},\bm{b}}(\bm{x})\end{pmatrix}\quad\forall\bm{x},\bm{p}\in\mathbb{R}^{n},\\\
&\mathcal{G}_{low}\begin{pmatrix}\bm{x}\\\
\bm{p}\end{pmatrix}=\begin{pmatrix}\hat{\sigma}_{K,\bm{a},\bm{b}}(\bm{p})+\bm{x}\\\
\bm{p}\end{pmatrix}\quad\forall\bm{x},\bm{p}\in\mathbb{R}^{n},\end{split}$

where the learnable parameters are the matrix $K$ and the vectors
$\bm{a},\bm{b}\in\mathbb{R}^{l}$. The dimension $l$ (which is the dimension of
$\bm{a},\bm{b}$ as well as the number of rows in $K$) is a positive integer
that can be tuned, called the width of SympNet. In [50], it is proven that
G-SympNets are universal approximators within the family of symplectic maps.
Note that it is easy to obtain the inverse map of a G-SympNet, since we have
explicit formulas for $\mathcal{G}_{up}^{-1}$ and $\mathcal{G}_{low}^{-1}$
given as follows

(6) $\mathcal{G}_{up}^{-1}\begin{pmatrix}\bm{x}\\\
\bm{p}\end{pmatrix}=\begin{pmatrix}\bm{x}\\\
\bm{p}-\hat{\sigma}_{K,\bm{a},\bm{b}}(\bm{x})\end{pmatrix},\quad\mathcal{G}_{low}^{-1}\begin{pmatrix}\bm{x}\\\
\bm{p}\end{pmatrix}=\begin{pmatrix}\bm{x}-\hat{\sigma}_{K,\bm{a},\bm{b}}(\bm{p})\\\
\bm{p}\end{pmatrix}.$

##  3 SympOCnet for optimal control problems without state constraints

In this section, we consider the following optimal control problem without
state constraints

(7)
$\begin{split}&\min\left\\{\int_{0}^{T}\ell(\bm{x}(s),\bm{v}(s))ds\colon\bm{x}(0)=\bm{x}_{0},\bm{x}(T)=\bm{x}_{T},\dot{\bm{x}}(s)=\bm{v}(s)\in
U\,\forall s\in(0,T)\right\\}\\\ =\
&\min\left\\{\int_{0}^{T}\ell(\bm{x}(s),\dot{\bm{x}}(s))ds\colon\bm{x}(0)=\bm{x}_{0},\bm{x}(T)=\bm{x}_{T},\dot{\bm{x}}(s)\in
U\,\forall s\in(0,T)\right\\},\end{split}$

where $\bm{v}\colon[0,T]\to U$ is the control function taking values in the
control set $U\subseteq\mathbb{R}^{n}$, $\bm{x}\colon[0,T]\to\mathbb{R}^{n}$
is the corresponding trajectory, and $\ell\colon\mathbb{R}^{n}\times
U\to\mathbb{R}$ is called running cost.

It is well-known in the literature that the solution $\bm{x}$ to the optimal
control problem (7) satisfies the following Hamiltonian ODE system

(8)
$\begin{dcases}\dot{\bm{x}}(s)=\nabla_{\bm{p}}H(\bm{x}(s),\bm{p}(s)),&s\in[0,T],\\\
\dot{\bm{p}}(s)=-\nabla_{\bm{x}}H(\bm{x}(s),\bm{p}(s)),&s\in[0,T],\end{dcases}$

where the function $H\colon\mathbb{R}^{2n}\to\mathbb{R}$ is called Hamiltonian
defined by

(9) $H(\bm{x},\bm{p})=\sup_{\bm{v}\in
U}\\{\langle\bm{v},\bm{p}\rangle-\ell(\bm{x},\bm{v})\\},\quad\forall\bm{x},\bm{p}\in\mathbb{R}^{n}.$

A method to solve the optimal control problem (7) in high dimensions is to
solve the corresponding Hamiltonian ODE system (8) instead. However, when the
Hamiltonian is too complicated, the Hamiltonian ODEs may be hard to solve. To
tackle this difficulty, we propose a method called SympOCnet, which uses the
SympNet architecture to solve the optimal control problem (7) in high
dimensions. The key idea is to perform a change of variables in the phase
space using a symplectic map represented by a SympNet, and then solve the
Hamiltonian ODEs in the new coordinate system.

Given a symplectic map $\phi\colon\mathbb{R}^{2n}\to\mathbb{R}^{2n}$, we
define the change of variables $(\bm{y},\bm{q})=\phi(\bm{x},\bm{p})$ for any
$\bm{x},\bm{p}\in\mathbb{R}^{n}$. If the phase trajectory
$t\mapsto(\bm{x}(t),\bm{p}(t))$ is the solution to the Hamiltonian ODEs (8),
then by Theorem 2.2, the phase trajectory $t\mapsto(\bm{y}(t),\bm{q}(t))$
under the new coordinate system solves the Hamiltonian ODEs with a new
Hamiltonian $\tilde{H}\colon\mathbb{R}^{2n}\to\mathbb{R}$ defined by

$\tilde{H}(\bm{y},\bm{q})=H(\phi^{-1}(\bm{y},\bm{q}))\quad\forall\bm{y},\bm{q}\in\mathbb{R}^{n}.$

In other words, the function $t\mapsto(\bm{y}(t),\bm{q}(t))$ satisfies

(10)
$\begin{dcases}\dot{\bm{y}}(s)=\nabla_{\bm{q}}\tilde{H}(\bm{y}(s),\bm{q}(s)),&s\in[0,T],\\\
\dot{\bm{q}}(s)=-\nabla_{\bm{y}}\tilde{H}(\bm{y}(s),\bm{q}(s)),&s\in[0,T].\end{dcases}$

Here, we assume that the new Hamiltonian $\tilde{H}$ has a simpler form such
that the corresponding Hamiltonian ODE system (10) is easier to solve. To be
specific, we assume that $\tilde{H}$ does not depend on the variable $\bm{y}$,
and hence the Hamiltonian ODEs become

(11) $\begin{dcases}\dot{\bm{y}}(s)=\nabla\tilde{H}(\bm{q}(s)),&s\in[0,T],\\\
\dot{\bm{q}}(s)=0,&s\in[0,T].\end{dcases}$

The solution $s\mapsto\bm{q}(s)$ is a constant function, and the solution
$s\mapsto\bm{y}(s)$ is an affine function. Therefore, the solution
$s\mapsto(\bm{y}(s),\bm{q}(s))$ can be represented using three parameters in
$\mathbb{R}^{n}$, which are denoted by $\bm{q}_{0}=\bm{q}(0)$,
$\bm{y}_{0}=\bm{y}(0)$ and $\bm{u}=\dot{\bm{y}}(s)$. With these three
parameters, the solution to (11) is given by $\bm{y}(s)=\bm{y}_{0}+s\bm{u}$
and $\bm{q}(s)=\bm{q}_{0}$ for any $s\in[0,T]$.

Under this framework, we explain our proposed method in the remainder of this
section. In Section 3.1, we apply the SympOCnet method to approximate the
unknown function $\phi$ and the unknown parameters $\bm{u},\bm{y}_{0}$ and
$\bm{q}_{0}$. Then, after the training process, to further enhance the
accuracy we can also apply the pseudospectral method for post-processing (for
low-dimensional problems), as we explain in Section 3.2.

###  3.1 SympOCnet method

In this section, we describe how to apply SympOCnet method to solve the
optimal control problem (7). SympOCnet uses a SympNet architecture to
represent the inverse of the unknown function $\phi$, which is also a
symplectic map. In other words, we approximate $\phi^{-1}$ using a SympNet
denoted by $\varphi=(\varphi_{1},\varphi_{2})$. Here, $\varphi_{1}$ denotes
the first $n$ components in $\varphi$, which is a map from $(\bm{y},\bm{q})$
to $\bm{x}$, and $\varphi_{2}$ denotes the last $n$ components in $\varphi$,
which is a map from $(\bm{y},\bm{q})$ to $\bm{p}$. Other learnable parameters
include $\bm{u},\bm{y}_{0},\bm{q}_{0}\in\mathbb{R}^{n}$ mentioned above, which
respectively correspond to $\dot{\bm{y}}(s)$, $\bm{y}(0)$ and $\bm{q}(s)$ in
the new Hamiltonian ODE system (11). With the SympNet $\varphi$ and the
variables $\bm{u},\bm{y}_{0},\bm{q}_{0}$, the solution to the original
Hamiltonian ODEs (8) is given by
$s\mapsto\varphi(\bm{y}_{0}+s\bm{u},\bm{q}_{0})$.

To learn the SympNet $\varphi$ and the parameters $\bm{u},\bm{y}_{0}$ and
$\bm{q}_{0}$, we apply the PINN loss [76] on the Hamiltonian ODEs (8), and get

$\begin{split}\mathcal{L}_{res}&=\sum_{j=1}^{N}\left\lVert\frac{d\bm{x}(s_{j})}{ds}-\nabla_{\bm{p}}H(\bm{x}(s_{j}),\bm{p}(s_{j}))\right\rVert^{2}+\sum_{j=1}^{N}\left\lVert\frac{d\bm{p}(s_{j})}{ds}+\nabla_{\bm{x}}H(\bm{x}(s_{j}),\bm{p}(s_{j}))\right\rVert^{2}\\\
&=\sum_{j=1}^{N}\left\lVert\nabla_{\bm{y}}\varphi_{1}(\bm{y}_{0}+s_{j}\bm{u},\bm{q}_{0})\bm{u}-\nabla_{\bm{p}}H(\varphi(\bm{y}_{0}+s_{j}\bm{u},\bm{q}_{0}))\right\rVert^{2}\\\
&\quad\quad+\sum_{j=1}^{N}\left\lVert\nabla_{\bm{y}}\varphi_{2}(\bm{y}_{0}+s_{j}\bm{u},\bm{q}_{0})\bm{u}+\nabla_{\bm{x}}H(\varphi(\bm{y}_{0}+s_{j}\bm{u},\bm{q}_{0}))\right\rVert^{2},\end{split}$

where the second equality holds by the following chain rule

$\begin{split}\frac{d\bm{x}(s_{j})}{ds}&=\nabla_{\bm{y}}\bm{x}(\bm{y}_{0}+s_{j}\bm{u},\bm{q}_{0})\frac{d\bm{y}(s_{j})}{ds}+\nabla_{\bm{q}}\bm{x}(\bm{y}_{0}+s_{j}\bm{u},\bm{q}_{0})\frac{d\bm{q}(s_{j})}{ds}\\\
&=\nabla_{\bm{y}}\varphi_{1}(\bm{y}_{0}+s_{j}\bm{u},\bm{q}_{0})\bm{u},\\\
\frac{d\bm{p}(s_{j})}{ds}&=\nabla_{\bm{y}}\bm{p}(\bm{y}_{0}+s_{j}\bm{u},\bm{q}_{0})\frac{d\bm{y}(s_{j})}{ds}+\nabla_{\bm{q}}\bm{p}(\bm{y}_{0}+s_{j}\bm{u},\bm{q}_{0})\frac{d\bm{q}(s_{j})}{ds}\\\
&=\nabla_{\bm{y}}\varphi_{2}(\bm{y}_{0}+s_{j}\bm{u},\bm{q}_{0})\bm{u}.\end{split}$

Here, the data points are given by $s_{1},\dots,s_{N}\in[0,T]$. The gradients
$\nabla_{\bm{y}}\varphi_{1}(\bm{y}_{0}+s_{j}\bm{u},\bm{q}_{0})$ and
$\nabla_{\bm{y}}\varphi_{2}(\bm{y}_{0}+s_{j}\bm{u},\bm{q}_{0})$ are calculated
by the automatic differentiation of the SympNet $\varphi$. Moreover, we have
another loss term for the initial and terminal conditions of $\bm{x}$, defined
by

$\begin{split}\mathcal{L}_{bd}&=\left\lVert\bm{x}(0)-\bm{x}_{0}\right\rVert^{2}+\left\lVert\bm{x}(T)-\bm{x}_{T}\right\rVert^{2}\\\
&=\left\lVert\varphi_{1}(\bm{y}_{0},\bm{q}_{0})-\bm{x}_{0}\right\rVert^{2}+\left\lVert\varphi_{1}(\bm{y}_{0}+T\bm{u},\bm{q}_{0})-\bm{x}_{T}\right\rVert^{2}.\end{split}$

The total loss function $\mathcal{L}$ is a weighted sum of the two loss terms,
defined by

(12) $\mathcal{L}=\mathcal{L}_{res}+\lambda\mathcal{L}_{bd},$

where $\lambda$ is a tunable positive hyperparameter. After training, the
optimal trajectory is obtained from

(13) $\bm{x}(s)=\varphi_{1}(\bm{y}_{0}+s\bm{u},\bm{q}_{0})\quad\forall
s\in[0,T],$

where the function $\varphi_{1}$ contains the first $n$ components of the
SympNet $\varphi$. An illustration of our proposed SympOCnet method is shown
in Figure 1.

Figure 1: An illustration of the SympOCnet method. We propose to use SympNet
$\phi$ to map curvy trajectory in the phase space to simpler trajectory in new
coordinates and solve the corresponding Hamiltonian system of the optimal
control problem.

###  3.2 Post-processing using the pseudospectral method

In Section 3.1, we introduced the SympOCnet to solve the optimal control
problem (7). In this section, we explain how to apply the pseudospectral
method to postprocess SympOCnet results to improve feasibility and optimality.
The pseudospectral method is a popular method proposed to solve a general
optimal control problem [31, 32, 77, 80, 38]. It applies the quadrature rule
to discretize the integration and the ODE in a continuous optimal control
problem. In this way, the continuous optimal control problem is converted to a
finite-dimensional optimization problem with constraints. Then, some
optimization algorithms such as sequential quadratic programming are applied
to solve the finite-dimensional problem.

When the dimension is high, it is, in general, hard to solve the converted
finite-dimensional optimization problem. The difficulty comes from the number
of variables and constraints in the optimization problem, as well as the non-
convexity of the problem. However, there are some theoretical guarantees when
the initialization is good enough. Therefore, it can be applied as a
postprocessing procedure by taking the trajectory obtained from Section 3.1 as
an initialization.

We will apply this procedure in the numerical experiments in Sections 5.1 and
5.2. From the numerical results, we observe that the pseudospectral method
preserves the shape of the trajectories obtained from the training step of the
SympOCnet method, and it makes the trajectories smoother while improving the
optimality. Therefore, the pseudospectral method performs well as a post-
processing method. It is also observed in the experiments that pseudospectral
method itself with a linear initialization may not converge to the correct
solution, which justifies the choice of SympOCnet as a good initializer. Note
that we do not apply this post-process to the experiments in Sections 5.3 and
5.4 due to the high dimensionality.

##  4 SympOCnet for optimal control problems with state constraints

In this section, we focus on the following optimal control problem with state
constraints

(14)
$\begin{split}&\min\Bigg{\\{}\int_{0}^{T}\ell(\bm{x}(s),\bm{v}(s))ds\colon
h(\bm{x}(s))\geq 0,\dot{\bm{x}}(s)=\bm{v}(s)\in U\forall s\in[0,T],\\\
&\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\bm{x}(0)=\bm{x}_{0},\bm{x}(T)=\bm{x}_{T}\Bigg{\\}}\\\
=\,&\min\Bigg{\\{}\int_{0}^{T}\ell(\bm{x}(s),\dot{\bm{x}}(s))ds\colon
h(\bm{x}(s))\geq 0,\dot{\bm{x}}(s)\in U\forall s\in[0,T],\\\
&\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\bm{x}(0)=\bm{x}_{0},\bm{x}(T)=\bm{x}_{T}\Bigg{\\}},\end{split}$

where $\ell\colon\mathbb{R}^{n}\times U\to\mathbb{R}$ is the running cost, and
$h\colon\mathbb{R}^{n}\to\mathbb{R}^{m}$ is a function providing the
constraint on the state variable $\bm{x}$. For instance, to avoid obstacles,
$h$ can be defined using the signed distance function to the obstacles.

In order to apply the SympOCnet method in Section 3, we need to convert the
original problem (14) to an unconstrained optimal control problem, which is
achieved using the soft penalty method described in Section 4.1. The soft
penalty method will enforce the constraints in an asymptotic sense, i.e., when
the hyperparameters converge to zero under some conditions, which is
impractical in experiments. In the training process, when the hyperparameters
in the soft penalty are fixed, to improve the feasibility of the output
trajectory, we introduce extra terms in the loss function in Sections 4.2 and
4.3. The selection of these extra loss terms may differ from problem to
problem. Later in Section 5.1, we compare the performance of different loss
functions under our path planning problem setup to choose a suitable loss
function for the other experiments.

###  4.1 Soft penalty method for the optimal control problem

To convert a constrained problem to an unconstrained one, the soft penalty
method replaces the hard constraint by a soft penalty term in the objective
function. Here, we consider the log penalty function [82, 85], denoted by
$\beta_{a}\colon\mathbb{R}\to\mathbb{R}$ and defined by

(15) $\beta_{a}(x):=\begin{dcases}-\log(x),&\text{if}\ x>a,\\\
-\log(a)+\frac{1}{2}\left(\left(\frac{x-2a}{a}\right)^{2}-1\right),&\text{if}\
x\leq a,\end{dcases}$

where $a$ is a positive hyperparameter. The high-dimensional log penalty
function is the summation of the one-dimensional function acting on each
component, i.e., for any $m>1$, the log penalty function
$\beta_{a}\colon\mathbb{R}^{m}\to\mathbb{R}$ is defined by

(16)
$\beta_{a}(\bm{x}):=\sum_{i=1}^{m}\beta_{a}(x_{i})\quad\forall\bm{x}=(x_{1},\dots,x_{m})\in\mathbb{R}^{m}.$

With this log penalty term, the problem (14) is converted to

(17)
$\begin{split}\min\Bigg{\\{}\int_{0}^{T}\ell(\bm{x}(s),\dot{\bm{x}}(s))+\epsilon\beta_{a}(h(\bm{x}(s))ds\colon\dot{\bm{x}}(s)\in
U\forall s\in[0,T],\quad\quad\\\
\bm{x}(0)=\bm{x}_{0},\bm{x}(T)=\bm{x}_{T}\Bigg{\\}},\end{split}$

where $\epsilon$ and $a$ are positive hyperparameters. This problem (17) is in
the form of (7), and hence can be solved using the SympOCnet method in Section
3. By straightforward calculation, for a sequence of parameters $\epsilon_{k}$
and $a_{k}$ such that they both converge to zero and the following equations
hold

$\lim_{k\to\infty}\epsilon_{k}\log
a_{k}=0,\quad\quad\lim_{k\to\infty}\frac{a_{k}^{2}}{\epsilon_{k}}=0,$

the penalty term $\epsilon_{k}\beta_{a_{k}}$ converges to the indicator
function of $[0,+\infty)$, which takes the value $0$ in $[0,+\infty)$ and
$+\infty$ in $(-\infty,0)$. Therefore, the objective function in (17) with
parameters $\epsilon_{k}$, $a_{k}$ converges to the objective function in (14)
with the hard constraint given by $h$.

With this penalty term, the Hamiltonian corresponding to the optimal control
problem (17) equals

(18) $H_{\epsilon,a}(\bm{x},\bm{p})=\sup_{\bm{v}\in
U}\\{\langle\bm{v},\bm{p}\rangle-\ell(\bm{x},\bm{v})-\epsilon\beta_{a}(h(\bm{x}))\\}=H(\bm{x},\bm{p})-\epsilon\beta_{a}(h(\bm{x})),$

where $H$ is the Hamiltonian corresponding to the running cost $\ell$ defined
in (17).

######  Remark 4.1.

Note that another widely-used method to covert a constrained optimization
problem to an unconstrained one is the augmented Lagrangian method. In our
problem, since the constraint is enforced for any time $s\in[0,T]$, to apply
the augmented Lagrangian method, the Lagrange multiplier needs to be a
function of time. Then, with the added Lagrangian term, the new running cost
also depends on the time variable, which requires a more complicated
symplectic method to handle time-dependent Hamiltonians. Therefore, we do not
apply the augmented Lagrangian method for the constrained optimal control
problem (14) in this paper. Instead, we add extra terms in the loss function
(see Sections 4.2 and 4.3) to enforce the constraint.

###  4.2 Penalty method in the training process of SympOCnet

To enforce the state constraint $h(\bm{x}(s))\geq 0$, one way is to add the
penalty term in the loss function of the training process. Here, we introduce
two penalty terms, which include the log penalty $\beta_{a}$ defined in (16)
and the quadratic penalty. The loss function with the log penalty term
corresponding to the state constraint $h(\bm{x}(s))\geq 0$ is defined by

(19)
$\begin{split}\mathcal{L}_{log}(\varphi,\bm{y}_{0},\bm{q}_{0},\bm{u})&=\mathcal{L}(\varphi,\bm{y}_{0},\bm{q}_{0},\bm{u})+\frac{\tilde{\lambda}}{N}\sum_{j=1}^{N}\epsilon\beta_{a}(h(\varphi_{1}(\bm{y}_{0}+s_{j}\bm{u},\bm{q}_{0})))\\\
&\quad\quad+\epsilon\beta_{a}(\bm{x}_{0}-\varphi_{1}(\bm{y}_{0},\bm{q}_{0}))+\epsilon\beta_{a}(\varphi_{1}(\bm{y}_{0},\bm{q}_{0})-\bm{x}_{0})\\\
&\quad\quad+\epsilon\beta_{a}(\bm{x}_{T}-\varphi_{1}(\bm{y}_{0}+T\bm{u},\bm{q}_{0}))+\epsilon\beta_{a}(\varphi_{1}(\bm{y}_{0}+T\bm{u},\bm{q}_{0})-\bm{x}_{T}),\end{split}$

where $\beta_{a}$ is the penalty function defined in (16), $\mathcal{L}$ is
the loss function defined in Section 3.1, and $\tilde{\lambda}$ is a positive
hyperparameter. Similarly, the loss function with a quadratic penalty is
defined by

(20)
$\mathcal{L}_{quad}(\varphi,\bm{y}_{0},\bm{q}_{0},\bm{u})=\mathcal{L}(\varphi,\bm{y}_{0},\bm{q}_{0},\bm{u})+\frac{\tilde{\lambda}}{N}\sum_{j=1}^{N}\|\min\\{h(\varphi_{1}(\bm{y}_{0}+s_{j}\bm{u},\bm{q}_{0})),0\\}\|^{2},$

where we do not add the penalty term for the initial and terminal conditions,
since the quadratic penalty for them are the same with $\mathcal{L}_{bd}$ in
the loss function $\mathcal{L}$.

###  4.3 Augmented Lagrangian method in the training process of SympOCnet

Another method to enforce the state constraint $h(\bm{x}(s))\geq 0$ is to
apply augmented Lagrangian method in the training process. The loss function
with the Lagrangian term is defined as follows

(21)
$\begin{split}&\mathcal{L}_{aug}(\varphi,\bm{y}_{0},\bm{q}_{0},\bm{u},\bm{\mu},\bm{\lambda}_{1},\bm{\lambda}_{2})\\\
=\,&\mathcal{L}(\varphi,\bm{y}_{0},\bm{q}_{0},\bm{u})+\frac{1}{2\rho_{1}}\sum_{j=1}^{N}\|\max\\{\bm{0},\bm{\mu}(s_{j})-\rho_{1}h(\varphi_{1}(\bm{y}_{0}+s_{j}\bm{u},\bm{q}_{0}))\\}\|^{2}\\\
&\quad+\frac{1}{2\rho_{2}}\left(\|\bm{\lambda}_{1}-\rho_{2}(\bm{x}_{0}-\varphi_{1}(\bm{y}_{0},\bm{q}_{0}))\|^{2}+\|\bm{\lambda}_{2}-\rho_{2}(\bm{x}_{T}-\varphi_{1}(\bm{y}_{0}+T\bm{u},\bm{q}_{0}))\|^{2}\right),\end{split}$

where $\bm{\mu}(s_{j})\in\mathbb{R}^{m}$,
$\bm{\lambda}_{1},\bm{\lambda}_{2}\in\mathbb{R}^{n}$ are the augmented
Lagrange multipliers for the constraints $h(\bm{x}(s_{j}))\geq 0$,
$\bm{x}(0)=\bm{x}_{0}$ and $\bm{x}(T)=\bm{x}_{T}$, respectively, and
$\rho_{1},\rho_{2}>0$ are positive hyperparameters in the augmented Lagrangian
method. The training process becomes an iterative scheme where the $k$-th
iteration contains the following two steps

  * •

Train the SympNet $\varphi$ and parameters
$\bm{y}_{0},\bm{q}_{0},\bm{u}\in\mathbb{R}^{n}$ for several iterations using
the loss function
$\mathcal{L}_{aug}(\varphi,\bm{y}_{0},\bm{q}_{0},\bm{u},\bm{\mu}^{k},\bm{\lambda}_{1}^{k},\bm{\lambda}_{2}^{k})$;

  * •

Update the Lagrange multiplier by

$\begin{split}\bm{\mu}^{k+1}(s_{j})&=\max\\{\bm{0},\bm{\mu}^{k}(s_{j})-\rho_{1}h(\varphi_{1}(\bm{y}_{0}+s_{j}\bm{u},\bm{q}_{0}))\\}\quad\forall
j=1,\dots,N,\\\
\bm{\lambda}_{1}^{k+1}&=\bm{\lambda}_{1}^{k}-\rho_{2}(\bm{x}_{0}-\varphi_{1}(\bm{y}_{0},\bm{q}_{0})),\\\
\bm{\lambda}_{2}^{k+1}&=\bm{\lambda}_{2}^{k}-\rho_{2}(\bm{x}_{T}-\varphi_{1}(\bm{y}_{0}+T\bm{u},\bm{q}_{0})).\end{split}$

Note that in practice, we need to tune the hyperparameters $\rho_{1}$ and
$\rho_{2}$. In our implementation, these hyperparameters are tuned
automatically using a similar strategy as in [16].

##  5 Applications in path planning problems with obstacle and collision
avoidance

In this section, we apply our method to path planning problems with multiple
drones. We assume that each drone is represented by a ball in the physical
space $\mathbb{R}^{m}$ with radius $C_{d}$. Note that the meaning of the
notation $m$ in this section (i.e., the dimension of the physical space) is
different from $m$ in Section 4 (which is the number of state constraints). We
set the state variable to be
$\bm{x}=(\bm{x}_{1},\dots,\bm{x}_{M})\in\mathbb{R}^{Mm}$, where $M$ is the
number of drones, and each $\bm{x}_{j}\in\mathbb{R}^{m}$ denotes the position
of the center of each drone. In other words, the state variable is the
concatenation of the position variables of all drones. Similarly, the control
variable $\bm{v}$ is the concatenation of the velocity variables for all
drones, i.e., $\bm{v}$ equals $(\bm{v}_{1},\dots,\bm{v}_{M})\in U^{M}$, where
each $\bm{v}_{j}$ is the velocity of the $j$-th drone, and
$U\subseteq\mathbb{R}^{m}$ is the space for all admissible velocities for each
drone. In our numerical experiments, the control space $U$ is the ball in
$\mathbb{R}^{m}$ with zero center and radius $C_{\bm{v}}>0$, i.e., the norm of
the velocity of each drone has the upper bound $C_{\bm{v}}$.

We want to solve for the optimal trajectory with minimal energy, and hence the
running cost is set to be

$\ell(\bm{x},\bm{v})=\frac{1}{2}\|\bm{v}\|^{2}\quad\forall\bm{x}\in\mathbb{R}^{Mm},\bm{v}\in
U^{M}.$

The corresponding Hamiltonian equals

$\begin{split}H(\bm{x},\bm{p})&=\sup_{\bm{v}\in
U^{M}}\left\\{\langle\bm{v},\bm{p}\rangle-\frac{1}{2}\|\bm{v}\|^{2}\right\\}=\sum_{i=1}^{M}\sup_{\bm{v}_{i}\in
U}\left\\{\langle\bm{v}_{i},\bm{p}_{i}\rangle-\frac{1}{2}\|\bm{v}_{i}\|^{2}\right\\}\\\
&=\sum_{i=1}^{M}\begin{dcases}\frac{1}{2}\|\bm{p}_{i}\|^{2}&\text{if
}\|\bm{p}_{i}\|\leq C_{\bm{v}},\\\
C_{\bm{v}}\|\bm{p}_{i}\|-\frac{C_{\bm{v}}^{2}}{2}&\text{if
}\|\bm{p}_{i}\|>C_{\bm{v}},\end{dcases}\end{split}$

for any $\bm{p}=(\bm{p}_{1},\cdots,\bm{p}_{M})\in\mathbb{R}^{Mm}$, where each
$\bm{p}_{i}$ is the momentum variable in $\mathbb{R}^{m}$ for the $i$-th
drone.

To avoid obstacles and collisions among drones, we set the constraint function
$h$ to be $h=(h_{1},h_{2})$, where $h_{1}$ is for obstacle avoidance, and
$h_{2}$ is for avoiding collisions among drones. If there are $n_{o}$
obstacles, denoted by $E_{1},\dots,E_{n_{o}}$, then the function
$h_{1}\colon\mathbb{R}^{Mm}\to\mathbb{R}^{Mn_{o}}$ is defined by

(22)
$h_{1}(\bm{x}_{1},\dots,\bm{x}_{M})=(d(\bm{x}_{1}),d(\bm{x}_{2}),\dots,d(\bm{x}_{M}))\quad\forall\bm{x}_{1},\dots,\bm{x}_{M}\in\mathbb{R}^{m},$

where the function $d\colon\mathbb{R}^{m}\to\mathbb{R}^{n_{o}}$ is defined by

(23)
$d(\bm{x})=(d_{1}(\bm{x}),\dots,d_{n_{o}}(\bm{x}))\quad\forall\bm{x}\in\mathbb{R}^{m},$

and each function $d_{j}$ is defined such that $d_{j}(\bm{x})<0$ implies the
collision of the drone whose center is at the position $\bm{x}$ with the
$j$-th obstacle $E_{j}$. For instance, $d_{j}$ can be defined using the signed
distance function to $E_{j}$. The definition of the function $d_{j}$ depends
on the shape of the constraint set $E_{j}$, and hence we do not specify it now
but postpone its definition to each example in later sections.

The function $h_{2}\colon\mathbb{R}^{Mm}\to\mathbb{R}^{M(M-1)/2}$ is the
constraint function for collision avoidance; each component of the constraint
function gives a constraint for avoiding collision between a pair of drones.
Since each drone can be seen as a ball centered at a point in
$\mathbb{R}^{m}$, two drones collide if and only if the distance between their
centers is less than the sum of their radii. Hence, we set the $k$-th
component of $h_{2}$ to be

(24)
$(h_{2})_{k}(\bm{x}_{1},\dots,\bm{x}_{M})=\|\bm{x}_{i}-\bm{x}_{j}\|^{2}-(2C_{d})^{2}\quad\forall\bm{x}_{1},\dots,\bm{x}_{M}\in\mathbb{R}^{m},$

where $C_{d}$ is the radius of a drone, and $k=i+(j-1)(j-2)/2$ for any $1\leq
i<j\leq M$ is the corresponding constraint index for the pair $(i,j)$. Note
that the constraints may be duplicated to simplify the implementation. As a
result, $(h_{2})_{k}(\bm{x}_{1},\dots,\bm{x}_{M})<0$ is equivalent to the
collision between the $i$-th and $j$-th drones. In other words, such defined
constraint function $h_{2}$ can avoid collisions among drones.

In the following sections, we present several numerical results. In Section
5.1, we compare the performance of different loss functions $\mathcal{L}$,
$\mathcal{L}_{log}$, $\mathcal{L}_{quad}$ and $\mathcal{L}_{aug}$ to choose a
suitable loss function for the remaining experiments. In Section 5.2, we
present an offline-online training strategy for an optimal control problem
whose initial conditions are not known in the training process of the SympNet.
This experiment shows the generalizability of SympOCnet. Then, in Section 5.3,
we demonstrate the performance of SympOCNet on a high-dimensional problem
reported in [79]. From the results in Section 5.3, we observe that SympOCnet
can handle the path planning problem whose state space has dimension $512$,
and hence the method can potentially mitigate the CoD. In Section 5.4, we
apply SympOCnet to a swarm path planning problem in [74], and demonstrate good
performance and efficiency in path planning problems, where the agents move in
a three-dimensional space.

In Sections 5.1 and 5.2 we use a $6$-layer SympNet with $60$ hidden neurons in
each layer and ReLU activation function. In Sections 5.3 and 5.4 we use a
$6$-layer SympNet with $200$ hidden neurons in each layer and ReLU activation
function. In all these experiments, the parameter $\lambda$ in (12) is set to
be $600$, and the parameter $\tilde{\lambda}$ in (19) and (20) is set to be
$200$. By default, the time horizon $T$ is set to be $1$, and the speed limit
$C_{\bm{v}}$ is $25$. All the numerical experiments are run using a shared
NVIDIA GeForce RTX 3090 GPU and a shared NVIDIA RTX A6000 GPU. In each
experiment, we only show the figure for a part of the testing cases. More
numerical results and codes are provided in
https://github.com/zzhang222/SympOCNet. Video animations of these examples are
available online at
https://github.com/zzhang222/SympOCNet/tree/main/saved_results_animation.

###  5.1 Comparison among different loss functions

In this section, we consider the path planning problem with two obstacles and
four drones in $\mathbb{R}^{2}$. Each obstacle $E_{j}\subset\mathbb{R}^{2}$ is
a convex set containing points whose distance to a line segment (denoted by
$l_{j}\subset\mathbb{R}^{2}$) is less than a constant $C_{o}$. The initial
positions of the four drones are $(-2,-2)$, $(2,-2)$, $(2,2)$, $(-2,2)$, and
their terminal positions are $(2,2),(-2,2),(-2,-2),(2,-2)$. Therefore, the
optimal control problem (14) is an $8$-dimensional problem, whose initial
state $\bm{x}_{0}$ is $(-2,-2,2,-2,2,2,-2,2)$ and terminal state $\bm{x}_{T}$
is $(2,2,-2,2,-2,-2,2,-2)$. The obstacles and the initial positions of the
four drones are shown in Figure 2.

Figure 2: Initial positions of drones and obstacles. The four drones located
at their initial positions are represented by the four colored circles. The
two obstacles are represented by round cornered rectangles in the middle.

NN, $t=\frac{1}{3}$

NN, $t=\frac{2}{3}$

NN, $t=1$

PS, $t=\frac{1}{3}$

PS, $t=\frac{2}{3}$

PS, $t=1$

Figure 3: Output trajectories with loss $\mathcal{L}_{aug}$ and pseudospectral
method. The left three figures show the trajectories of our proposed SympOCnet
method at different time $t=\frac{1}{3},\frac{2}{3},1$, and the right three
figures show the outputs of the post-process (pseudospectral method) at
different time $t=\frac{1}{3},\frac{2}{3},1$. In each figure, the four colored
circles show the positions of the four drones at time $t$, and the curves
connected to the circles show the trajectories before time $t$.

To avoid the collisions of drones with obstacles, we define the first part of
the constraint function $h$ by (22), where each component $d_{j}$ of the
function $d$ in (23) is defined by

$d_{j}(\bm{x})=\min_{\bm{y}\in
l_{j}}\|\bm{x}-\bm{y}\|^{2}-(C_{o}+C_{d})^{2}\quad\forall\bm{x}\in\mathbb{R}^{2}.$

The term $\min_{\bm{y}\in l_{j}}\|\bm{x}-\bm{y}\|^{2}$ in $d_{j}$ gives the
squared distance between $\bm{x}$ and the line segment $l_{j}$. Here, we use
the squared distance instead of distance function to improve the smoothness of
the constraint function $h$. Since the obstacle $E_{j}$ contains all points
whose distance to the line segment $l_{j}$ is less than $C_{o}$, a drone
collide with $E_{j}$ if and only if the distance between its center and the
line segment $l_{j}$ is less than $C_{o}+C_{d}$. Therefore, such defined
function $d_{j}$ provides a constraint which avoids the collisions of drones
with the $j$-th obstacle $E_{j}$.

Under this problem setup, we compare four different loss functions, i.e., the
function $\mathcal{L}$ defined in (12), $\mathcal{L}_{log}$ in (19),
$\mathcal{L}_{quad}$ in (20), and $\mathcal{L}_{aug}$ in (21). With each loss
function, we train the neural network for $100,000$ iterations, and then run
the pseudospectral method to improve the results. By comparing the performance
of the four loss functions in this example, we choose a better loss function
for the other experiments in the remainder of this paper.

The outputs of the four loss functions are shown in Table 1, where the
hyperparameters in the penalty term $\epsilon\beta_{a}(h(\bm{x}(s))$ are set
to be $a=0.004$ and $\epsilon=0.0004$. Results for different loss functions
are shown in different columns. For each loss function, we repeat the
experiment $10$ times with different random seeds, and then take average to
obtain the statistics in the table. We observe the convergence of the
pseudospectral method in all repeated experiments, which shows the robustness
of our SympOCnet method, providing optimal or suboptimal solutions in this
example. The second to fourth lines in Table 1 show the minimal constraint
value, the cost value in the optimal control problem, and the running time of
the training process of our proposed SympOCnet method. The last four lines in
Table 1 show the minimal constraint value, the cost value in the optimal
control problem, the number of iterations, and the running time of the post-
process using pseudospectral method.

loss function | $\mathcal{L}$ | $\mathcal{L}_{log}$ | $\mathcal{L}_{quad}$ | $\mathcal{L}_{aug}$
---|---|---|---|---
SympOCnet | min constraint | -0.834 | -0.00704 | -0.0824 | -4.84E-04
cost | 79.62 | 94.15 | 76.50 | 93.98
running time (s) | 1229 | 1776 | 1361 | 1708
PS | min constraint | -4.70E-09 | -9.87E-12 | -2.04E-12 | -8.02E-13
cost | 73.47 | 91.00 | 76.95 | 80.91
# iterations | 49 | 19 | 17 | 23
running time (s) | 104 | 44 | 38 | 47
Table 1: The comparison among results of four loss functions with SympOCnet method and pseudospectral post-process. We show the minimal constraint value $\min_{s}h(\bm{x}(s))$, the cost $\int_{s=0}^{T}\frac{||\bm{v}(s)||^{2}}{2}ds$, and the running time of the SympOCnet as well as those statistics after the post-process. The loss function $\mathcal{L}_{aug}$ provides the solution with least amount of constraint violations and reasonable cost values. It can be seen that in all cases, the SympOCnet provides a good initialization for the pseudospectral method.  minimal constraint | cost | number of iterations | running time (s)
---|---|---|---
-1.00E+00 | 64.00 | 5 | 44
Table 2: The result of the pseudospectral method with the linear
initialization. We show the minimal constraint value $\min_{s}h(\bm{x}(s))$,
the cost $\int_{s=0}^{T}\frac{||\bm{v}(s)||^{2}}{2}ds$, the number of
iterations, and the running time of the pseudospectral method with the linear
initialization (i.e., the initial trajectory is an affine function of the time
variable, and the initial velocity is a constant function). The minimal
constraint value is $-1.00$, which shows that the constraint is not satisfied,
and the pseudospectral method with the linear initialization does not converge
to a feasible solution.

Comparing the results of the four loss functions, we observe that the penalty
term and the augmented Lagrangian term in the loss function indeed improve the
feasibility of the obtained trajectory. The constraints are better satisfied
using the loss $\mathcal{L}_{aug}$ in (21), while the corresponding cost value
may be a little bit higher than the the results of $\mathcal{L}_{quad}$. In
the other experiments, the problems are more complicated, and the constraints
are harder to satisfy. Therefore, we use the loss function $\mathcal{L}_{aug}$
for the experiments in the remainder of the paper to help enforce the state
constraints. For simpler problems, other loss function may be more preferable
for faster running time. How to choose the loss function adaptively may be a
possible future direction. Although the augmented Lagrangian method and log
penalty loss function are slightly slower than the other two loss functions,
it takes less than $30$ minutes to run $100,000$ training iterations, which
shows the efficiency of SympOCnet.

To show the improvement using our proposed SympOCnet method as an
initialization of the pseudospectral method, we solve the same problem using
the pseudospectral method with the linear initialization, i.e., we set the
initial guess of the pseudospectral method to be
$\bm{v}(t)=\frac{\bm{x}_{T}-\bm{x}_{0}}{T}$ and
$\bm{x}(t)=\frac{(\bm{x}_{T}-\bm{x}_{0})t}{T}$ for any $t\in[0,T]$. The
minimal constraint value, the cost value in the optimal control problem, the
number of iterations, and the running time are shown in Table 2. The minimal
constraint value in Table 2 is $-1.00$, which shows that the pseudospectral
method with the linear initialization does not converge to a feasible
solution. Therefore, compared with the linear initialization, our SympOCnet
method, no matter which loss function in $\mathcal{L}$, $\mathcal{L}_{log}$,
$\mathcal{L}_{quad}$, and $\mathcal{L}_{aug}$ we choose, provides a better
initial guess for the pseudospectral method.

Moreover, we show the obtained results in one trial using SympOCnet with loss
$\mathcal{L}_{aug}$ and the post-process using the pseudospectral method in
Figure 3. The left three figures show the trajectories of our proposed
SympOCnet method at different times, while the right three figures correspond
to the final results from the pseudospectral method at different times. In
each figure, the four colored circles show the positions of the four drones at
the specific time, and the curves connected to the circles show the
trajectories before the current time. We observe that the shape of the
trajectory given by our SympOCnet method does not change too much after the
post-process, which shows that SympOCnet provides a suboptimal solution to the
optimal control problem in this example.

###  5.2 Generalizability and offline-online training

case 1

case 2

case 3

case 4

case 5

case 6

Figure 4: Initial positions of drones, obstacles, and training data. In the
six figures, we show the initial positions of the four drones (represented by
colored circles), the positions of the two obstacles (represented by the round
cornered rectangles), and the randomly sampled training data (represented by
the yellow dots) in the six testing cases in Section 5.2.

NN, $t=\frac{1}{3}$

NN, $t=\frac{2}{3}$

NN, $t=1$

PS, $t=\frac{1}{3}$

PS, $t=\frac{2}{3}$

PS, $t=1$

NN, $t=\frac{1}{3}$

NN, $t=\frac{2}{3}$

NN, $t=1$

PS, $t=\frac{1}{3}$

PS, $t=\frac{2}{3}$

PS, $t=1$

Figure 5: Output trajectories of offline-online training and post-process. The
first and second rows correspond to the first and fourth case, respectively.
The left three columns are the output from L-BFGS, and the right three columns
are the output from the pseudospectral method. In each figure, the four
colored circles show the positions of the four drones at time $t$, and the
curves connected to the circles show the trajectories before time $t$.

In this section, we consider the same problem setup as in Section 5.1. To test
the generalizability of SympOCnet, we assume that the exact initial position
is unknown in the training process. We train one SympNet $\varphi$ using
multiple initial positions sampled from a uniform distribution centered at the
initial position in Section 5.1. With multiple initial positions, we train one
SympNet $\varphi$ and several variables
$\\{(\bm{y}_{0k},\bm{u}_{k},\bm{q}_{0k})\\}_{k}$, where the $k$-th trainable
tuple $(\bm{y}_{0k},\bm{u}_{k},\bm{q}_{0k})$ corresponds to the $k$-th sampled
initial position. We refer to this training process of the SympNet and all
trainable variables as “offline training”. After $100,000$ steps of offline
training, the exact initial positions are provided. To obtain an improved
solution in reasonable running time using the trained SympNet $\varphi$ as
well as the new information, we fix the SympNet $\varphi$ and apply the L-BFGS
method to train the variables $(\bm{y}_{0},\bm{u},\bm{q}_{0})$ with the exact
initial and terminal positions. We call this training process using L-BFGS
method the “online training”. The optimal trajectory is then computed using
(13), where $\varphi_{1}$ comes from the SympNet $\varphi$ trained in the
offline training, and the parameters $(\bm{y}_{0},\bm{u},\bm{q}_{0})$ are
outputted from the online training. After the online training, we also apply
the pseudospectral method to improve the quality of the results.

To generate training data for the initial positions, we sample $100$ points
from the uniform distribution in $[x-1,x+1]\times[y-1,y+1]$ where $(x,y)$ is
the initial position of a drone in Section 5.1. To test the generalizability
of our method, we have six cases. In the first three cases, the exact initial
positions are randomly sampled from the same distribution used to generate the
training data. In other words, we assume the exact data is covered in the area
of the training data. In the last three cases, we set the initial conditions
for the state variable to be $(-2,-4,2,-4,2,4,-2,4)$, $(-4,-4,0,-4,0,4,-4,4)$,
and $(0,-4,0,-2,0,2,0,4)$, respectively. These three initial positions are not
in the range of the training data. Therefore, the results with these three
initial positions show the generalizability of our SympOCnet method when the
test data is not covered by the area of the training data. The initial
positions of the drones, the positions of the obstacles, and the training data
are shown in Figure 4, where the four colored circles denote the drones, the
two round cornered rectangles denote the obstacles, and the yellow dots denote
the sampled training data.

In our experiments, we found that the performance is improved if we increase
the dimension of the problem by adding some latent variables. We increase the
dimension of the state space from $8$ to $16$, and define the Hamiltonian on
the larger space by

$H(\bm{x}_{1},\bm{x}_{2},\bm{p}_{1},\bm{p}_{2})=H_{\epsilon,a}(\bm{x}_{1},\bm{p}_{1})+\frac{1}{2}\|\bm{p}_{2}\|^{2}\quad\forall\bm{x}_{1},\bm{x}_{2},\bm{p}_{1},\bm{p}_{2}\in\mathbb{R}^{8},$

where $(\bm{x}_{1},\bm{p}_{1})$ are the variables in the original problem,
$(\bm{x}_{2},\bm{p}_{2})$ are the latent variables we add, and
$H_{\epsilon,a}$ is the function defined in (18). The minimal constraint, cost
in the optimal control problem, and running time of the L-BFGS method and
pseudospectral method are shown in Table 3. As in Section 5.1, we run $10$
repeated experiments and put their average values in the table. In the online
training, we run $1,000$ L-BFGS iterations to train the parameters of all the
six cases simultaneously. Note that we can train them all at once, since they
share the same SympNet, and the loss function is the summation of the loss
functions in all cases. Then, we run the pseudospectral method on each case
until it converges. The statistics for the L-BFGS training and the
pseudospectral method for the six cases are shown in each line of Table 3. The
minimal constraint value of L-BFGS method is not as good as the results in
Section 5.1, which is expected, since the initial positions are not known when
the SympNet is trained. However, from the minimal constraint values of the
pseudospectral post-process in the six cases, we conclude that our offline-
online training still provides a good initialization for the pseudospectral
method, which shows a good generalizability of our proposed method.

| minimal constraint | cost | time (s)
---|---|---|---
L-BFGS | -1.27E-04 | 231.44 | 158
PS (case 1) | -1.97E-11 | 102.00 | 59
PS (case 2) | -1.74E-11 | 94.59 | 77
PS (case 3) | -3.36E-11 | 90.95 | 79
PS (case 4) | -6.00E-09 | 141.22 | 65
PS (case 5) | -7.19E-08 | 140.47 | 111
PS (case 6) | 1.64E-13 | 112.41 | 67
Table 3: The results given by the offline-online training (L-BFGS) and the post-process (pseudospectral method) in six generalizability tests. The offline-online training using L-BFGS method provides a good initialization for the pseudospectral method. With this initialization, the pseudospectral method in the post-process improves the feasibility and optimality of the solution in all the six cases. | minimal constraint | cost | time (s)
---|---|---|---
PS (linear initialization, case 1) | -7.71E-09 | 85.19 | 134
PS (linear initialization, case 2) | 2.95E-07 | 83.98 | 173
PS (linear initialization, case 3) | -1.94E-10 | 79.67 | 178
PS (linear initialization, case 4) | -5.55E-16 | 80.44 | 157
PS (linear initialization, case 5) | -8.08E-14 | 136.16 | 147
PS (linear initialization, case 6) | -1.00 | 60.00 | 28
Table 4: The results given by the pseudospectral method with the linear
initialization in six generalizability tests. Although the pseudospectral
method with the linear initialization performs slightly better than our
proposed method in the first five cases, it violates the constraint in the
last case. Therefore, our proposed method provides a more robust
initialization for the pseudospectral method compared with the linear
initialization.

We also compare our results with the pseudospectral method whose
initialization is given by the linear interpolation (as descibed in Section
5.1) in the six generalizability tests. The minimal constraint value, cost
value in the optimal control problem, and running time of the pseudospectral
method with the linear initialization are shown in Table 4. Comparing Table 3
with Table 4, we see that although our proposed method has slightly larger
cost values in the first five cases, the pseudospectral method with the linear
initialization does not converge to a feasible solution in the last case.
Therefore, our proposed method provides a more robust initialization for the
pseudospectral method compared with the linear initialization.

Moreover, we plot the output trajectories of one trial of the first and fourth
cases computed using our proposed method in Figure 5. The left three columns
show the output trajectory given by the offline-online training, and the right
three columns are the output trajectory from the post-process. In each figure,
the yellow dots are the initial positions of the sampled training data.
Similar as in Section 5.1, the four colored circles show the current positions
of the four drones, and the curves connected to them illustrate the
trajectories before the current time. The first row in Figure 5 is for the
first case where the exact initial positions are in the area of the sampled
training data. The second row is for the fourth case where the exact initial
positions are outside the area of the training data. We observe that the
offline-online training process still provides feasible trajectories in both
cases, based on which the post-process improves the smoothness of the
trajectories to provide more optimal solutions. Therefore, our proposed
SympOCnet method has the generalizability to handle unseen data which are
close to the training data.

###  5.3 Multiple drones in a two-dimensional room

In this section, we consider the path planning problem for multiple drones in
a room, which is inspired by [79]. To be specific, we consider $M$ drones
moving in a room $[-C_{r},C_{r}]^{2}$ for some $C_{r}>0$. In our experiments,
we set the constant $C_{r}=5$ (i.e., the room has size $10\times 10$). We need
to avoid the collisions among the drones as well as the collisions between
drones and the walls. The constraint function $h=(h_{1},h_{2})$ contains two
parts: the function $h_{2}$ is defined by (24), and the function $h_{1}$ is
the obstacle constraint provided by the four walls. Hence, we set $h_{1}$ in
the form of (22), where the function $d$ is defined by
$d(\bm{x})=(x_{1}+C_{r},C_{r}-x_{1},x_{2}+C_{r},C_{r}-x_{2})$ for any
$\bm{x}=(x_{1},x_{2})\in\mathbb{R}^{2}$. The initial positions of the drones
are near the boundary of the room, and the terminal positions are the opposite
locations, i.e., we set $\bm{x}_{T}=-\bm{x}_{0}$. This is a high-dimensional
example (the dimension of the state space is $n=2M$, which ranges from $64$ to
$512$ in our experiments), and hence we apply our SympOCnet method without the
post-process.

First, we set the drone radius $C_{d}$ to be $0.3$. We run $10$ repeated
experiments (with $100,000$ training iterations in each experiment) with drone
numbers $M=32$ and $64$ to test the robustness of our SympOCnet method. The
results are shown in Table 5. In each trial, we compute the minimal normalized
distance of any pair of drones at any time grid, i.e., we compute
$\mathcal{D}$ defined by

$\mathcal{D}=\frac{1}{2C_{d}}\min_{1\leq i<j\leq M}\min_{1\leq k\leq
N}\|\bm{x}_{i}(t_{k})-\bm{x}_{j}(t_{k})\|,$

where $\bm{x}_{i}(t_{k})$ denotes the center position for the $i$-th drone at
time $t_{k}$. Then, we compute the average and standard deviation of these
minimal normalized distances $\mathcal{D}$ in the $10$ trials, which are shown
in the second and third columns in Table 5. We also compute the scaled cost in
the optimal control problem, which is defined by
$\frac{1}{M}\int_{0}^{T}\frac{1}{2}\|\dot{\bm{x}}(t)\|^{2}dt$, where
$\bm{x}(\cdot)$ is the output trajectory for the state variable $\bm{x}$ in
the optimal control problem. The average of the scaled cost is shown in the
fourth column in Table 5. In the last column of Table 5, we show the averaged
running time of the training process. From the results, we observe that our
SympOCnet method provides feasible results in reasonable running time for this
high-dimensional problem. Note that the collision between two drones happens
whenever the minimal normalized distance is less than $1$. Therefore, in
general the collision does not happen for $32$ drones, but it may happen for
$64$ or more drones. To provide a more feasible solution, in the following
experiments, we put a safe zone outside each drone by setting $C_{d}$ to be a
little bit larger than the actual drone radius.

# drones | $\mathbb{E}(\mathcal{D})$ | std($\mathcal{D}$) | scaled cost | time (s)
---|---|---|---|---
32 | 1.001 | 0.00390 | 85.84 | 1905
64 | 0.986 | 0.00347 | 100.43 | 2214
Table 5: The results for the path planning problem with $32$ and $64$ drones
in a room. We show the expected value and standard deviation of the minimal
normalized distances. The collision does not happen for 32 drones, but it may
happen for 64 or more drones. Therefore, we need to put a safe zone near each
drone to avoid collision when the number of drones is large. The computation
time increases only by 309 seconds when scaling from 32 to 64 drones.

We also tested our SympOCnet method on this path planning problem with $128$
drones (whose radius is $0.075$) and four obstacles. The initial and terminal
positions are similar to the cases of $32$ or $64$ drones (see Figure 6 (a)
and (d)). The four obstacles are balls inside the room, which are represented
by the black circles in Figure 6. The constraint function $d_{j}$ in (23)
corresponding to the obstacle $E_{j}$ is defined by
$d_{j}(\bm{x})=\|\bm{z}_{j}-\bm{x}\|^{2}-(C_{o}+C_{d})^{2}$ for any
$\bm{x}\in\mathbb{R}^{n}$, where $\bm{z}_{j}\in\mathbb{R}^{2}$ and $C_{o}>0$
are the center and radius of the obstacle $E_{j}$, respectively. The results
are shown in Figure 6. Each figure shows the current positions (represented by
colored circles) and the previous trajectories (represented by grey curves) of
the drones at different time $t=0,\frac{1}{3},\frac{2}{3},1$. To provide a
feasible solution, we set $C_{d}$ in the constraint function to be $0.1$
instead of $0.075$. It takes $4161$ seconds to finish $100,000$ training
iterations. From the numerical results, we found that there are no collisions
between obstacles and drones, and hence SympOCnet performs well and
efficiently in this high-dimensional problem.

$t=0$

$t=\frac{1}{3}$

$t=\frac{2}{3}$

$t=1$

Figure 6: Path planning of 128 drones with four obstacles on the plane. We
plot the predicted positions of 128 drones at time
$t=0,\frac{1}{3},\frac{2}{3},1$. The four black circles represent the four
obstacles, and the colored circles represent the drones. The paths from the
initial conditions to the current positions of all drones are plotted as gray
lines. There are no collisions in the three clipped snapshots, and the planned
trajectories are all inside the $10\times 10$ square.

Then, we tested our proposed method on a higher dimensional problem. We set
the drone number $M=256$, and hence the dimension of the state space is
$n=2M=512$. We solve this path planning problem with $256$ drones (whose
radius is $0.09$) in the two-dimensional room $[-5,5]^{2}$. The initial
positions are shown in Figure 7 (a), and the terminal condition is
$\bm{x}_{T}=-\bm{x}_{0}$. The results are shown in Figure 7. To provide a
feasible solution, we add a safe zone and set $C_{d}$ in the state constraint
to be $0.1$ instead of $0.09$. From this result, we observe no collisions
among the drones. It takes $5481$ seconds to obtain this feasible solution in
$100,000$ training iterations. Therefore, SympOCnet is able to efficiently
solve this $512$-dimensional path planning problem.

$t=0$

$t=\frac{1}{3}$

$t=\frac{2}{3}$

$t=1$

Figure 7: Path planning of 256 drones on the plane. We plot the predicted
positions of 256 drones at time $t=0,\frac{1}{3},\frac{2}{3},1$. The drones
are represented by colored circles. The paths from the initial positions to
the current positions of all drones are plotted as gray lines. There are no
collisions in the three clipped snapshots, and the planned trajectories are
inside the $10\times 10$ square.

###  5.4 Multiple drones with obstacle avoidance in a three-dimensional space

$t=\frac{1}{3}$

$t=\frac{2}{3}$

$t=1$

Figure 8: Path planning of 100 drones in a 3d space. We plot the predicted
positions of 100 drones at time $t=\frac{1}{3},\frac{2}{3},1$. The paths from
the initial positions to the current positions of all drones are plotted as
colored lines. The destinations are marked as red crosses. The drones reach
their destinations without any collision.

We consider the three-dimensional swarm path planning example in [74]. To be
specific, we consider $M=100$ drones with radius $0.18$. Their goal is to
avoid two obstacles which are located between their initial positions and
terminal positions, as illustrated in Figure 8. The total dimension of the
state space is $n=3M=300$. The $j$-th obstacle $E_{j}$ is defined to be the
rectangular set
$[C^{j}_{11},C^{j}_{12}]\times[C^{j}_{21},C^{j}_{22}]\times[C^{j}_{31},C^{j}_{32}]$,
where $C^{j}_{ik}$ is a constant scalar. We set the constraint function
$h=(h_{1},h_{2})$ in the same way as described before, where the function
$d_{j}$ in (23) is defined by

$\begin{split}&d_{j}(\bm{x})=\max\\{C^{j}_{i1}-C_{d}-x_{i},x_{i}-C^{j}_{i2}-C_{d}\colon
i=1,2,3\\},\end{split}$

for each $j=1,2$ and $\bm{x}=(x_{1},x_{2},x_{3})\in\mathbb{R}^{3}$. Note that
$d_{j}(\bm{x})\leq 0$ if and only if
$x_{i}\in[C^{j}_{i1}-C_{d},C^{j}_{i2}+C_{d}]$ for each $i=1,2,3$, which holds
if the ball centered at $\bm{x}$ with radius $C_{d}$ intersects with the
obstacle $E_{j}$. Therefore, the constraint $d_{j}(\bm{x})\geq 0$ enforces the
drones to avoid collisions with the obstacle $E_{j}$. In our experiment, we
set

$\begin{split}(C^{1}_{11},C^{1}_{12},C^{1}_{21},C^{1}_{22},C^{1}_{31},C^{1}_{32})&=(-1.8,1.8,-0.3,0.3,0.2,6.8),\\\
(C^{2}_{11},C^{2}_{12},C^{2}_{21},C^{2}_{22},C^{2}_{31},C^{2}_{32})&=(2.2,3.8,-0.8,0.8,0.2,3.8).\end{split}$

In other words, the constraint sets are
$[-1.8,1.8]\times[-0.3,0.3]\times[0.2,6.8]$ and
$[2.2,3.8]\times[-0.8,0.8]\times[0.2,3.8]$. We apply our proposed SympOCnet
method to solve this $300$-dimensional problem, and plot the result in Figure
8. Due to the high dimensionality, we do not apply the post-process. To
provide a feasible solution, we set $C_{d}$ in the state constraint to be
$0.2$ instead of $0.18$. From the numerical results, we observe no collisions
among the drones and the obstacles. The minimal distance between every pair of
drones is $0.3994$, which is bigger than twice of the radius. It takes $3633$
seconds to finish $150,000$ training iterations. Therefore, our proposed
SympOCnet method provides a feasible solution to this high-dimensional swarm
path planning problem in reasonable time.

##  6 Summary

We have proposed a novel SympOCnet method for efficiently solving high-
dimensional optimal control problems with state constraints. We applied
SympOCnet to several multi-agent simultaneous path planning problems with
obstacle avoidance. The numerical results show SympOCnet’s ability to solve
high-dimensional problems efficiently with dimension more than $500$. These
first results reveal the potential of SympOCnet for solving high-dimensional
optimal control problems in real-time. In future work, we are going to
consider possible combinations with the DeepONet architecture [64] to avoid
any training cost during inference and to endow the predictive system with
uncertainties, such as uncertain initial and terminal positions in path
planning problems.

## Acknowledgement

The simulations were run on GeForce RTX 3090 and RTX A6000 GPU donated to us
by NVIDIA.

## References

  * [1] M. Akian, R. Bapat, and S. Gaubert, Max-plus algebra, Handbook of linear algebra, 39 (2006).
  * [2] M. Akian, S. Gaubert, and A. Lakhoua, The max-plus finite element method for solving deterministic optimal control problems: basic properties and convergence analysis, SIAM Journal on Control and Optimization, 47 (2008), pp. 817–848.
  * [3] A. Alla, M. Falcone, and L. Saluzzi, An efficient DP algorithm on a tree-structure for finite horizon optimal control problems, SIAM Journal on Scientific Computing, 41 (2019), pp. A2384–A2406.
  * [4] A. Alla, M. Falcone, and S. Volkwein, Error analysis for POD approximations of infinite horizon problems via the dynamic programming approach, SIAM Journal on Control and Optimization, 55 (2017), pp. 3091–3115.
  * [5] A. Bachouch, C. Huré, N. Langrené, and H. Pham, Deep neural networks algorithms for stochastic control problems on finite horizon: numerical applications, arXiv preprint arXiv:1812.05916, (2018).
  * [6] S. Bansal and C. J. Tomlin, Deepreach: A deep learning approach to high-dimensional reachability, in 2021 IEEE International Conference on Robotics and Automation (ICRA), IEEE, 2021, pp. 1817–1824.
  * [7] R. E. Bellman, Adaptive control processes: a guided tour, Princeton university press, 1961.
  * [8] D. P. Bertsekas, Reinforcement learning and optimal control, Athena Scientific, Belmont, Massachusetts, (2019).
  * [9] O. Bokanowski, N. Gammoudi, and H. Zidani, Optimistic Planning Algorithms For State-Constrained Optimal Control Problems.  working paper or preprint, July 2021, https://hal.archives-ouvertes.fr/hal-03283075.
  * [10] O. Bokanowski, J. Garcke, M. Griebel, and I. Klompmaker, An adaptive sparse grid semi-Lagrangian scheme for first order Hamilton-Jacobi Bellman equations, Journal of Scientific Computing, 55 (2013), pp. 575–605.
  * [11] M. Chen, J. F. Fisac, S. Sastry, and C. J. Tomlin, Safe sequential path planning of multi-vehicle systems via double-obstacle Hamilton-Jacobi-Isaacs variational inequality, in 2015 European Control Conference (ECC), IEEE, 2015, pp. 3304–3309.
  * [12] M. Chen, Q. Hu, J. F. Fisac, K. Akametalu, C. Mackin, and C. J. Tomlin, Reachability-based safety and goal satisfaction of unmanned aerial platoons on air highways, Journal of Guidance, Control, and Dynamics, 40 (2017), pp. 1360–1373.
  * [13] M. Chen and C. J. Tomlin, Exact and efficient Hamilton-Jacobi reachability for decoupled systems, in 2015 54th IEEE Conference on Decision and Control (CDC), IEEE, 2015, pp. 1297–1303.
  * [14] P. Chen, J. Darbon, and T. Meng, Hopf-type representation formulas and efficient algorithms for certain high-dimensional optimal control problems, arXiv preprint arXiv:2110.02541, (2021).
  * [15] P. Chen, J. Darbon, and T. Meng, Lax-Oleinik-type formulas and efficient algorithms for certain high-dimensional optimal control problems, arXiv preprint arXiv:2109.14849, (2021).
  * [16] A. R. Conn, N. I. Gould, and P. Toint, A globally convergent augmented Lagrangian algorithm for optimization with general constraints and simple bounds, SIAM Journal on Numerical Analysis, 28 (1991), pp. 545–572.
  * [17] M. Coupechoux, J. Darbon, J.-M. Kélif, and M. Sigelle, Optimal trajectories of a UAV base station using Lagrangian mechanics, in IEEE INFOCOM 2019-IEEE Conference on Computer Communications Workshops (INFOCOM WKSHPS), IEEE, 2019, pp. 626–631.
  * [18] J. Darbon, On convex finite-dimensional variational methods in imaging sciences and Hamilton–Jacobi equations, SIAM Journal on Imaging Sciences, 8 (2015), pp. 2268–2293, https://doi.org/10.1137/130944163, https://arxiv.org/abs/https://doi.org/10.1137/130944163.
  * [19] J. Darbon, P. M. Dower, and T. Meng, Neural network architectures using min plus algebra for solving certain high dimensional optimal control problems and Hamilton-Jacobi PDEs, arXiv preprint arXiv:2105.03336, (2021).
  * [20] J. Darbon, G. P. Langlois, and T. Meng, Overcoming the curse of dimensionality for some Hamilton-Jacobi partial differential equations via neural network architectures, Res. Math. Sci., 7 (2020), p. 20, https://doi.org/10.1007/s40687-020-00215-6, https://doi.org/10.1007/s40687-020-00215-6.
  * [21] J. Darbon and T. Meng, On decomposition models in imaging sciences and multi-time Hamilton–Jacobi partial differential equations, SIAM Journal on Imaging Sciences, 13 (2020), pp. 971–1014, https://doi.org/10.1137/19M1266332, https://doi.org/10.1137/19M1266332, https://arxiv.org/abs/https://doi.org/10.1137/19M1266332.
  * [22] J. Darbon and T. Meng, On some neural network architectures that can represent viscosity solutions of certain high dimensional Hamilton–Jacobi partial differential equations, Journal of Computational Physics, 425 (2021), p. 109907, https://doi.org/https://doi.org/10.1016/j.jcp.2020.109907, http://www.sciencedirect.com/science/article/pii/S0021999120306811.
  * [23] J. Darbon, T. Meng, and E. Resmerita, On Hamilton-Jacobi PDEs and image denoising models with certain non-additive noise, arXiv preprint arXiv:2105.13997, (2021).
  * [24] J. Darbon and S. Osher, Algorithms for overcoming the curse of dimensionality for certain Hamilton-Jacobi equations arising in control theory and elsewhere, Res Math Sci Research in the Mathematical Sciences, 3 (2016), pp. 1–26, https://doi.org/10.1186/s40687-016-0068-7, https://doi.org/10.1186/s40687-016-0068-7.
  * [25] D. Delahaye, S. Puechmorel, P. Tsiotras, and E. Féron, Mathematical models for aircraft trajectory design: A survey, in Air Traffic Management and Systems, Springer, 2014, pp. 205–247.
  * [26] J. Denk and G. Schmidt, Synthesis of a walking primitive database for a humanoid robot using optimal control techniques, in Proceedings of IEEE-RAS International Conference on Humanoid Robots, 2001, pp. 319–326.
  * [27] B. Djeridane and J. Lygeros, Neural approximation of PDE solutions: An application to reachability computations, in Proceedings of the 45th IEEE Conference on Decision and Control, Dec 2006, pp. 3034–3039, https://doi.org/10.1109/CDC.2006.377184.
  * [28] S. Dolgov, D. Kalise, and K. K. Kunisch, Tensor decomposition methods for high-dimensional Hamilton–Jacobi–Bellman equations, SIAM Journal on Scientific Computing, 43 (2021), pp. A1625–A1650, https://doi.org/10.1137/19M1305136, https://doi.org/10.1137/19M1305136, https://arxiv.org/abs/https://doi.org/10.1137/19M1305136.
  * [29] P. M. Dower, W. M. McEneaney, and H. Zhang, Max-plus fundamental solution semigroups for optimal control problems, in 2015 Proceedings of the Conference on Control and its Applications, SIAM, 2015, pp. 368–375.
  * [30] A. El Khoury, F. Lamiraux, and M. Taix, Optimal motion planning for humanoid robots, in 2013 IEEE International Conference on Robotics and Automation, IEEE, 2013, pp. 3136–3141.
  * [31] F. Fahroo and I. M. Ross, Costate estimation by a Legendre pseudospectral method, Journal of Guidance, Control, and Dynamics, 24 (2001), pp. 270–277.
  * [32] F. Fahroo and I. M. Ross, Direct trajectory optimization by a Chebyshev pseudospectral method, Journal of Guidance, Control, and Dynamics, 25 (2002), pp. 160–166.
  * [33] M. Fallon, S. Kuindersma, S. Karumanchi, M. Antone, T. Schneider, H. Dai, C. P. D’Arpino, R. Deits, M. DiCicco, D. Fourie, T. Koolen, P. Marion, M. Posa, A. Valenzuela, K.-T. Yu, J. Shah, K. Iagnemma, R. Tedrake, and S. Teller, An architecture for online affordance-based perception and whole-body planning, Journal of Field Robotics, 32 (2015), pp. 229–254.
  * [34] S. Feng, E. Whitman, X. Xinjilefu, and C. G. Atkeson, Optimization based full body control for the atlas robot, in 2014 IEEE-RAS International Conference on Humanoid Robots, IEEE, 2014, pp. 120–127.
  * [35] W. Fleming and W. McEneaney, A max-plus-based algorithm for a Hamilton–Jacobi–Bellman equation of nonlinear filtering, SIAM Journal on Control and Optimization, 38 (2000), pp. 683–710, https://doi.org/10.1137/S0363012998332433.
  * [36] K. Fujiwara, S. Kajita, K. Harada, K. Kaneko, M. Morisawa, F. Kanehiro, S. Nakaoka, and H. Hirukawa, An optimal planning of falling motions of a humanoid robot, in 2007 IEEE/RSJ International Conference on Intelligent Robots and Systems, IEEE, 2007, pp. 456–462.
  * [37] J. Garcke and A. Kröner, Suboptimal feedback control of PDEs by solving HJB equations on adaptive sparse grids, Journal of Scientific Computing, 70 (2017), pp. 1–28.
  * [38] D. Garg, Advances in global pseudospectral methods for optimal control, PhD thesis, University of Florida Gainesville, FL, 2011.
  * [39] S. Gaubert, W. McEneaney, and Z. Qu, Curse of dimensionality reduction in max-plus based approximation methods: Theoretical estimates and improved pruning algorithms, in 2011 50th IEEE Conference on Decision and Control and European Control Conference, IEEE, 2011, pp. 1054–1061.
  * [40] S. Greydanus, M. Dzamba, and J. Yosinski, Hamiltonian neural networks, in Advances in Neural Information Processing Systems, H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, eds., vol. 32, Curran Associates, Inc., 2019, https://proceedings.neurips.cc/paper/2019/file/26cd8ecadce0d4efd6cc8a8725cbd1f8-Paper.pdf.
  * [41] E. Hairer, M. Hochbruck, A. Iserles, and C. Lubich, Geometric numerical integration, Oberwolfach Reports, 3 (2006), pp. 805–882.
  * [42] J. Han, A. Jentzen, and W. E, Solving high-dimensional partial differential equations using deep learning, Proceedings of the National Academy of Sciences, 115 (2018), pp. 8505–8510, https://doi.org/10.1073/pnas.1718942115.
  * [43] M. Hofer, M. Muehlebach, and R. D’Andrea, Application of an approximate model predictive control scheme on an unmanned aerial vehicle, in 2016 IEEE International Conference on Robotics and Automation (ICRA), IEEE, 2016, pp. 2952–2957.
  * [44] M. B. Horowitz, A. Damle, and J. W. Burdick, Linear Hamilton Jacobi Bellman equations in high dimensions, in 53rd IEEE Conference on Decision and Control, IEEE, 2014, pp. 5880–5887.
  * [45] C. Huré, H. Pham, A. Bachouch, and N. Langrené, Deep neural networks algorithms for stochastic control problems on finite horizon: Convergence analysis, SIAM Journal on Numerical Analysis, 59 (2021), pp. 525–557, https://doi.org/10.1137/20M1316640, https://doi.org/10.1137/20M1316640, https://arxiv.org/abs/https://doi.org/10.1137/20M1316640.
  * [46] C. Huré, H. Pham, and X. Warin, Some machine learning schemes for high-dimensional nonlinear PDEs, arXiv preprint arXiv:1902.01599, (2019).
  * [47] F. Jiang, G. Chou, M. Chen, and C. J. Tomlin, Using neural networks to compute approximate and guaranteed feasible Hamilton-Jacobi-Bellman PDE solutions, arXiv preprint arXiv:1611.03158, (2016).
  * [48] L. Jin, S. Li, J. Yu, and J. He, Robot manipulator control using neural networks: A survey, Neurocomputing, 285 (2018), pp. 23 – 34, https://doi.org/https://doi.org/10.1016/j.neucom.2018.01.002, http://www.sciencedirect.com/science/article/pii/S0925231218300158.
  * [49] P. Jin, Z. Zhang, I. G. Kevrekidis, and G. E. Karniadakis, Learning Poisson systems and trajectories of autonomous systems via Poisson neural networks, arXiv preprint arXiv:2012.03133, (2020).
  * [50] P. Jin, Z. Zhang, A. Zhu, Y. Tang, and G. E. Karniadakis, SympNets: Intrinsic structure-preserving symplectic networks for identifying Hamiltonian systems, Neural Networks, 132 (2020), pp. 166–179.
  * [51] D. Kalise, S. Kundu, and K. Kunisch, Robust feedback control of nonlinear PDEs by numerical approximation of high-dimensional Hamilton-Jacobi-Isaacs equations, arXiv preprint arXiv:1905.06276, (2019).
  * [52] D. Kalise and K. Kunisch, Polynomial approximation of high-dimensional Hamilton–Jacobi–Bellman equations and applications to feedback control of semilinear parabolic PDEs, SIAM Journal on Scientific Computing, 40 (2018), pp. A629–A652.
  * [53] W. Kang and L. C. Wilcox, Mitigating the curse of dimensionality: sparse grid characteristics method for optimal feedback control and HJB equations, Computational Optimization and Applications, 68 (2017), pp. 289–315.
  * [54] Y. H. Kim, F. L. Lewis, and D. M. Dawson, Intelligent optimal control of robotic manipulators using neural networks, Automatica, 36 (2000), pp. 1355 – 1364, https://doi.org/https://doi.org/10.1016/S0005-1098(00)00045-5, http://www.sciencedirect.com/science/article/pii/S0005109800000455.
  * [55] M. R. Kirchner, M. J. Debord, and J. P. Hespanha, A Hamilton–Jacobi formulation for optimal coordination of heterogeneous multiple vehicle systems, in 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), IEEE, 2020, pp. 11623–11630.
  * [56] S. Kuindersma, R. Deits, M. Fallon, A. Valenzuela, H. Dai, F. Permenter, T. Koolen, P. Marion, and R. Tedrake, Optimization-based locomotion planning, estimation, and control design for the atlas humanoid robot, Autonomous robots, 40 (2016), pp. 429–455.
  * [57] K. Kunisch, S. Volkwein, and L. Xie, HJB-POD-based feedback design for the optimal control of evolution problems, SIAM Journal on Applied Dynamical Systems, 3 (2004), pp. 701–722.
  * [58] P. Lambrianides, Q. Gong, and D. Venturi, A new scalable algorithm for computational optimal control under uncertainty, Journal of Computational Physics, 420 (2020), p. 109710, https://doi.org/https://doi.org/10.1016/j.jcp.2020.109710, https://www.sciencedirect.com/science/article/pii/S0021999120304848.
  * [59] D. Lee and C. J. Tomlin, A Hopf-Lax formula in Hamilton–Jacobi analysis of reach-avoid problems, IEEE Control Systems Letters, 5 (2020), pp. 1055–1060.
  * [60] D. Lee and C. J. Tomlin, A Computationally Efficient Hamilton-Jacobi-based Formula for State-Constrained Optimal Control Problems, arXiv e-prints, (2021), arXiv:2106.13440, p. arXiv:2106.13440, https://arxiv.org/abs/2106.13440.
  * [61] F. Lewis, D. Dawson, and C. Abdallah, Robot Manipulator Control: Theory and Practice, Control engineering, Marcel Dekker, 2004, https://books.google.com/books?id=BDS_PQAACAAJ.
  * [62] A. Li, S. Bansal, G. Giovanis, V. Tolani, C. Tomlin, and M. Chen, Generating robust supervision for learning-based visual navigation using Hamilton-Jacobi reachability, in Learning for Dynamics and Control, PMLR, 2020, pp. 500–510.
  * [63] F. Lin and R. D. Brandt, An optimal control approach to robust control of robot manipulators, IEEE Transactions on robotics and automation, 14 (1998), pp. 69–77.
  * [64] L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis, Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators, Nature Machine Intelligence, 3 (2021), pp. 218–229.
  * [65] V. R. Makkapati, J. Ridderhof, P. Tsiotras, J. Hart, and B. van Bloemen Waanders, Desensitized trajectory optimization for hypersonic vehicles, in 2021 IEEE Aerospace Conference (50100), 2021, pp. 1–10, https://doi.org/10.1109/AERO50100.2021.9438511.
  * [66] W. McEneaney, A curse-of-dimensionality-free numerical method for solution of certain HJB PDEs, SIAM Journal on Control and Optimization, 46 (2007), pp. 1239–1276, https://doi.org/10.1137/040610830.
  * [67] W. M. McEneaney, Max-plus methods for nonlinear control and estimation, Systems & Control: Foundations & Applications, Birkhäuser Boston, Inc., Boston, MA, 2006.
  * [68] W. M. McEneaney, A. Deshpande, and S. Gaubert, Curse-of-complexity attenuation in the curse-of-dimensionality-free method for HJB PDEs, in 2008 American Control Conference, IEEE, 2008, pp. 4684–4690.
  * [69] W. M. McEneaney and L. J. Kluberg, Convergence rate for a curse-of-dimensionality-free method for a class of HJB PDEs, SIAM Journal on Control and Optimization, 48 (2009), pp. 3052–3079.
  * [70] X. Meng and G. E. Karniadakis, A composite neural network that learns from multi-fidelity data: Application to function approximation and inverse PDE problems, Journal of Computational Physics, 401 (2020), p. 109020.
  * [71] T. Nakamura-Zimmerer, Q. Gong, and W. Kang, Adaptive deep learning for high-dimensional Hamilton–Jacobi–Bellman equations, SIAM Journal on Scientific Computing, 43 (2021), pp. A1221–A1247, https://doi.org/10.1137/19M1288802, https://doi.org/10.1137/19M1288802, https://arxiv.org/abs/https://doi.org/10.1137/19M1288802.
  * [72] T. Nakamura-Zimmerer, Q. Gong, and W. Kang, QRnet: Optimal regulator design with LQR-augmented neural networks, IEEE Control Systems Letters, 5 (2021), pp. 1303–1308, https://doi.org/10.1109/LCSYS.2020.3034415.
  * [73] K. N. Niarchos and J. Lygeros, A neural approximation to continuous time reachability computations, in Proceedings of the 45th IEEE Conference on Decision and Control, Dec 2006, pp. 6313–6318, https://doi.org/10.1109/CDC.2006.377358.
  * [74] D. Onken, L. Nurbekyan, X. Li, S. W. Fung, S. Osher, and L. Ruthotto, A neural network approach for high-dimensional optimal control, arXiv preprint arXiv:2104.03270, (2021).
  * [75] C. Parzani and S. Puechmorel, On a Hamilton-Jacobi-Bellman approach for coordinated optimal aircraft trajectories planning, Optimal Control Applications and Methods, 39 (2018), pp. 933–948.
  * [76] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational Physics, 378 (2019), pp. 686–707.
  * [77] A. V. Rao, A survey of numerical methods for optimal control, Advances in the Astronautical Sciences, 135 (2009), pp. 497–528.
  * [78] C. Reisinger and Y. Zhang, Rectified deep neural networks overcome the curse of dimensionality for nonsmooth value functions in zero-sum games of nonlinear stiff systems, Analysis and Applications, 18 (2020), pp. 951–999, https://doi.org/10.1142/S0219530520500116, https://doi.org/10.1142/S0219530520500116, https://arxiv.org/abs/https://doi.org/10.1142/S0219530520500116.
  * [79] D. R. Robinson, R. T. Mar, K. Estabridis, and G. Hewer, An efficient algorithm for optimal trajectory generation for heterogeneous multi-agent systems in non-convex environments, IEEE Robotics and Automation Letters, 3 (2018), pp. 1215–1222.
  * [80] I. M. Ross and M. Karpenko, A review of pseudospectral optimal control: From theory to flight, Annual Reviews in Control, 36 (2012), pp. 182–197, https://doi.org/https://doi.org/10.1016/j.arcontrol.2012.09.002, https://www.sciencedirect.com/science/article/pii/S1367578812000375.
  * [81] V. R. Royo and C. Tomlin, Recursive regression with neural networks: Approximating the HJI PDE solution, arXiv preprint arXiv:1611.02739, (2016).
  * [82] A. Rucco, G. Notarstefano, and J. Hauser, An efficient minimum-time trajectory generation strategy for two-track car vehicles, IEEE Transactions on Control Systems Technology, 23 (2015), pp. 1505–1519.
  * [83] A. Rucco, P. Sujit, A. P. Aguiar, J. B. De Sousa, and F. L. Pereira, Optimal rendezvous trajectory for unmanned aerial-ground vehicles, IEEE Transactions on Aerospace and Electronic Systems, 54 (2017), pp. 834–847.
  * [84] J. Sirignano and K. Spiliopoulos, DGM: A deep learning algorithm for solving partial differential equations, Journal of Computational Physics, 375 (2018), pp. 1339 – 1364, https://doi.org/10.1016/j.jcp.2018.08.029.
  * [85] S. Spedicato and G. Notarstefano, Minimum-time trajectory generation for quadrotors in constrained environments, IEEE Transactions on Control Systems Technology, 26 (2017), pp. 1335–1344.
  * [86] E. Todorov, Efficient computation of optimal actions, Proceedings of the national academy of sciences, 106 (2009), pp. 11478–11483.
  * [87] I. Yegorov and P. M. Dower, Perspectives on characteristics based curse-of-dimensionality-free numerical approaches for solving Hamilton–Jacobi equations, Applied Mathematics & Optimization, (2017), pp. 1–49.
  * [88] Z. Zhang, Y. Shin, and G. E. Karniadakis, GFINNs: GENERIC formalism informed neural networks for deterministic and stochastic dynamical systems, arXiv preprint arXiv:2109.00092, (2021).
  * [89] M. Zhou, J. Han, and J. Lu, Actor-critic method for high dimensional static Hamilton–Jacobi–Bellman partial differential equations based on neural networks, arXiv preprint arXiv:2102.11379, (2021).