\(}{\left(}
\renewcommand{\)}{\right)}
\newcommand{\q}{\mathbf{q}}
\newcommand{\x}{\mathbf{x}}
\newcommand{\f}{\mathbf{f}}
\renewcommand{\u}{\mathbf{u}}
\newcommand{\w}{\mathbf{w}}
\newcommand{\nchoosek}[2]{\left(\begin{matrix}
\end{matrix}\right)}

\maketitle

\blfootnote{Sandia National Laboratories is a multimission laboratory managed and operated by National Technology \& Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525.}

\blfootnote{This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.}

\begin{abstract}
This paper presents a technique for trajectory planning based on continuously parameterized high-level actions (motion primitives) of variable duration.  This technique leverages deep reinforcement learning (Deep RL) to formulate a policy which is suitable for real-time implementation.  There is no separation of motion primitive generation and trajectory planning: each individual short-horizon motion is formed during the Deep RL training to achieve the full-horizon objective.  Effectiveness of the technique is demonstrated numerically on a well-studied trajectory generation problem and a planning problem on a known obstacle-rich map.  This paper also develops a new loss function term for policy-gradient-based Deep RL, which is analogous to an anti-windup mechanism in feedback control.  We demonstrate the inclusion of this new term in the underlying optimization increases the average policy return in our numerical example.
\end{abstract}

\section{Introduction}
Autonomous vehicles have become an active area of research over the past several decades.  One of the biggest challenges for these vehicles is trajectory planning, in which a complex optimization problem must be solved in real-time, reliably, and typically on a system with limited computational resources.  Many advances have been made in trajectory planning and optimization.  Examples include nonlinear programming methods [Schlossman et~al., 2021, Patterson and Rao, 2014], sampling based methods [LaValle and Kuffner~Jr, 2001, Karaman and Frazzoli, 2011, Kuwata et~al., 2009], and discrete motion planning such as A* [Hart et~al., 1968], Hybrid A* [Dolgov et~al., 2008], and Dijkstra's Algorithm [Dijkstra et~al., 1959].  However, existing methods still suffer from drawbacks.  Solving optimization problems in real-time can still be computationally challenging.  In addition, sampling-based methods and discrete motion planners must be modified to account for kinodynamic constraints.  This can be done by performing closed loop forward simulations (CLRRT) [Kuwata et~al., 2009] or using stereotyped behaviors such as motion primitives [Frazzoli et~al., 2000].  Motion primitives (MPs) are an attractive option because they are a type of high-level action that does not need to be specified at each time step of the horizon.  However, generating a motion primitive library can be time-consuming and non-intuitive.

Deep Reinforcement learning (RL) [Boyan and Moore, 1995, Bertsekas and Tsitsiklis, 1996, Sutton et~al., 1998, Mnih et~al., 2015]  offers the potential to improve motion planning performance [Faust et~al., 2018, Everett et~al., 2018, Aradi, 2020, G{\'o}mez et~al., 2012, Semnani et~al., 2020, Qureshi} et~al., 2021].  One method of using Deep RL in motion planning is to generate end-to-end trajectories for systems with complex dynamics and constraints.  This leverages a key benefit of Deep RL, which is its ability to perform well for problems where a nonlinear programming formulation is not suitably defined.  Similarly, Deep RL can provide benefit when the underlying problem is stochastic [Bellemare et~al., 2020].  Deep RL has recently been applied in aerospace domains.  Deep RL is applied to spacecraft orbit guidance in [LaFarge et~al., 2020] and to planetary descent in [Gaudet et~al., 2020] on a six degree of freedom system model.  In both works, rewards are carefully shaped to produce desirable behavior.  In [Hovell and Ulrich, 2020] Deep RL is applied to a spacecraft docking problem.  The recent work of [Sun and van Kampen, 2021] applies Deep RL for adaptive optimal flight control.

Despite much promise, Deep RL does not always provide effective end-to-end planning performance.  Specifically, RL can struggle when there are sparse rewards [Faust et~al., 2018] or on long horizon tasks [Qureshi} et~al., 2021].  In addition, achieving convergence often requires significant hyper-parameter tuning and reward shaping [Berner et~al., 2019].   As a result, many past works have used Deep RL to solve a smaller part of the planning problem.  For example, [Chiang} et~al., 2019] used Deep RL to estimate reachability, [Faust et~al., 2018] used RL for local planning, and [Goddard et~al., 2021, Li et~al., 2020, Camci and Kayacan, 2019]  used RL to learn high-level actions (primitives).

In this paper we also examine how navigation problems can be simplified so that Deep RL can be utilized effectively.  We propose a technique that uses a high-level action space. These actions are sub-trajectories of variable duration and shape.  High-level actions have been examined in a range of contexts for planning.  For example, the maneuver automaton consists of trims and maneuvers that can be combined to create end-to-end kinodynamic plans [Frazzoli et~al., 2000, Frazzoli, 2001, Frazzoli et~al., 2002, Gray et~al., 2012].  Similarly, motions can be created using dynamic movement primitives (DMPs) [Stulp and Schaal, 2011, Ijspeert et~al., 2002].  These primitives are represented using dynamical systems theory and can be modulated and combined to achieve complex behaviors.  Each DMP is characterized by an intermediate goal parameter and shape parameters which decide how the intermediate goal is reached.  DMPs have been applied to robotic flight in [Perk and Slotine, 2006] where individual primitives are learned directly from training trajectory demonstrations.

In our approach, we draw inspiration from dynamic movement primitives.  The action shape and duration are designed at each segment along the horizon, and each action can represent a variety of quantities including the profile of a high-level ``desired dynamics'' function which is provided to a lower-level tracking controller.  Separating the trajectory planning problem into high-level and low-level tasks potentially reduces the learning burden on the neural network, as the Deep RL agent is not required to memorize the full dynamics of the environment.  We show that our method is highly compatible with Deep RL, encouraging exploration within the environment, and we illustrate how only simple reward shaping is needed.  Rapid planning for quasi-static or even static maps remains a challenging problem for systems with highly complex dynamics and constraints [Grant and Antony, 2016].  Traditional approaches for solving these types of problems in real-time often rely on receding horizon control where an optimization problem is performed on-line utilizing modeling simplifications [Kuwata et~al., 2005, Schouwenaars et~al., 2004].  Alternatively, we utilize our Deep RL approach for trajectory planning, where a state-feedback policy is first learned through off-line training.  Afterwards, the trained policy can be evaluated rapidly to produce system inputs in response to arbitrary states.  Trajectory planning is demonstrated on two problems involving a shuttle reentry vehicle: the first problem is a free-space optimization problem, the second problem focuses on feasibility and incorporates a quasi-static map with regions of denied flight.  On the first problem we demonstrate that our planner can achieve near-optimality compared to an optimal control solution.  On the second problem we demonstrate that our planner can rapidly and robustly produce feasible solutions under the quasi-static map.  Additionally, we present a new loss function term for policy gradient Deep RL.  This new loss function term prevents the policy from producing actions which, on average, are outside action limits.  We describe how this method is analogous to the classical ``integrator windup'' problem in feedback control systems in which actuator limits are present.

\begin{figure}[h!]
\centering
\begin{minipage}{.99\textwidth}
\includegraphics[page=1, trim = 5mm 105mm 32mm 0mm, clip,width=.99\textwidth]{OverviewFig.pdf}
\end{minipage}
\captionsetup{format=hang}
\captionsetup{width=.9\textwidth}
\caption{Overview of the approach.}
\label{fig:OverviewFig}
\end{figure}

The core contributions of this work are 1) the description of a high-level action space for motion planning, 2) the illustration of this method promotes exploration and minimizes reward shaping, 3) the development of a new loss function term for the policy gradient, and 4) numerical validation of the proposed methods using a shuttle re-entry example.

This paper is organized as follows: A mathematical background is provided in Section II.  Section III describes the problem formulation where we describe our primary technique and frame the problem in the context of Deep RL.  We conclude the section with a short analysis describing the quality of solutions, and the influence the technique has on exploration in Deep RL problems.  Section IV describes a new mechanism to prevent ``action windup'' in the context of Deep RL policy gradient methods.  Numerical experiments are performed on a shuttle reentry problem in Section V.  We present our conclusions in Section VI.

\section{Mathematical Background}
In this section we provide a brief background on the mathematical concepts that are relevant for trajectory generation and Deep RL.
\subsection{Optimal Control}\label{section:Optimal Control}
Optimal control examines trajectory generation by considering the system dynamics and a reward function.  Optimal control seeks the time-series trajectory or feedback law that maximizes the reward function [Bryson and Ho, 1975] (or equivalently, minimizes a cost function).  The system dynamics are formulated as differential equations.  Consider a system with the following dynamics,
\begin{equation}\label{eq:dynamics}
\dot x = f(x(t),u(t))
\end{equation}
where $x \in \mathbb{R}^n$ is the state and $u \in \mathbb{R}^m$ is the control.  With this system we associate the following performance reward function [Bertsekas et~al., 1995, Lewis et~al., 2012]
\begin{align} \label{eq: J}
J &= \Phi(x(t_H)) + \int_{t_0}^{t_H} L(x(t),u(t)) dt
\end{align}
where typically the initial time $t_0$ is fixed and horizon time $t_H$ is a variable to be optimized.  The function $L(x(t), u(t)) \in \mathbb{R}$ measures the performance of the system along the horizon, while horizon reward function $\Phi(x(t_H)) \in \mathbb{R}$ associates a final reward received at the end of horizon.  The goal is to design a control trajectory, $u(t):[t_0,t_H]\rightarrow \mathbb{R}^m$, which maximizes (\ref{eq: J}) subject to dynamics (\ref{eq:dynamics}) while simultaneously satisfying the following constraint conditions:
\begin{enumerate}
\item $u(t) \in [u_{min}, u_{max}] \subset \mathbb{R}^m ~\forall~ t\in[t_0,t_H]$
\item $x(t)\in X \subset \mathbb{R}^n ~\forall~ t\in[t_0,t_H]$
\item $g_j(x(t), u(t), t) \leq 0, ~j=1,\dots,J, ~\forall~ t\in[t_0,t_H]$
\item $\Psi(x(t_H)) = 0, \Psi\in\mathbb{R}^q$
\end{enumerate}
The first two conditions ensures the control and state trajectories are \textit{admissible}, i.e., they do not violate any control or state constraints.  In this work we consider box form control constraints.  The third condition is a set of joint state-control path inequality constraints at each point along the horizon.  The last condition is a set of terminal constraints, which differs from the horizon reward function $\Phi(x(t_H))$.  In most cases the optimal control problem cannot be solved analytically.  Rather, so-called direct transcription methods [Kelly, 2017, Schlossman et~al., 2021] or indirect methods [Grant et~al., 2011] can be applied to solve the problem numerically, in which case the dynamics (\ref{eq:dynamics}) can be approximated with a numerical discretization of the form
\begin{equation} \label{eq:dynamics_discretetime}
x_{k+1} = F_k(x_k,u_k)
\end{equation}
Here $F_k$ is some time-discretized approximation to the dynamics function $f$, and $x_k=x(k\Delta t)$, $u_k=u(k\Delta t)$ for some time step $\Delta t$ and $k\in \mathbb{Z}$. A distinct feature which many of these methods share is they are \textit{model-based}, whereby gradients of the dynamics and objective function are computed along trajectories based on analytical models.

\subsection{Deep Reinforcement Learning}\label{section: deep rl}
Rather than designing an optimal and admissible control trajectory, RL [Sutton and Barto, 2018] seeks to design a feedback policy to maximize an associated performance index.  Unlike the model-based approach described in \ref{section:Optimal Control}, RL formulations can be \textit{model-free}, in which an \textit{agent} interacts directly with the \textit{environment}, periodically receiving rewards and state observations from the environment.  This is particularly valuable for environments that cannot be easily described by differential equations.  In this work, we refer to each step of the agent-environment interaction as an \textit{action step}, $h$, and distinguish it from the time step, $\Delta t$, associated with the discrete time approximation Eq. (\ref{eq:dynamics_discretetime}) of the dynamics.  It is through this direct interaction with the environment that the agent learns to adjust its policy to maximize the performance index.  In this work, we utilize RL to perform trajectory generation for nonlinear aerospace systems subject to complex constraints.

\subsubsection{Markov Decision Process}
In general, the RL problem is formalized by a discrete-time Markov Decision Process (MDP) consisting of the following components:
\begin{enumerate}
\item The state observation space, $\mathcal{S}$, where each state $s\in \mathcal{S}\subset \mathbb{R}^{\dim(s)}$.
\item The set of feasible actions the environment can accept, called the \textit{action space}, $\mathcal{A}$.  We restrict our attention to continuous actions spaces, where $a\in\mathcal{A} \subset \mathbb{R}^{\dim(a)}$ is suitable for continuous control problems.
\item The state transition probability density function $p(s'|s,a)$ which describes the likelihood of the agent observing state $s'$ at action step $h+1$ given that the agent took action $a$ from state $s$ at action step $h$.
\item The scalar function $r_h(s,a,s')$ is the expected reward received by the agent when taking action $a$ from state $s$ at action step $h$ and arriving at state $s'$ at action step $h+1$.
\end{enumerate}
Associated with the MDP is the discounted return,
\begin{equation}\label{eq:discounted return}
R_h = \sum_{k=h}^{H} \gamma^{k-h}r_{k}(s_k, a_k, s_{k+1})
\end{equation}
Here $0 < \gamma \leq 1$ is referred to as the discount factor and controls how much effect future rewards have on near-term optimal decisions (in [Schulman et~al., 2015], $\gamma$ is also interpreted as a variance reduction parameter for value function estimation).  We will show later that the discount factor brings additional benefits to our framework by encouraging exploration.  The agent-environment interactions may naturally separate into subsequences referred to as \textit{episodes}, in which case horizon length $H$ is a random variable which is determined by the agent reaching a terminal state.  This episodic case, in which there is a clear terminal objective, is our focus in this work.

Actions are generated according to some policy $\pi$, which may be deterministic, $a = \pi(s)$, or stochastic, $a \sim \pi(a|s)$ in which case $\pi(a|s)$ represents the probability of choosing action $a$ when the state is $s$.  For any MDP there is an optimal policy which is deterministic [Littman, 1994].  Nonetheless, defining the agent with a stochastic policy during training promotes exploration within the environment [Sutton and Barto, 2018].  The goal of an MDP is to design an optimal policy $\pi^* = \arg\max_\pi \E \left[ R_0 | \pi \right]$ which maximizes Eq. (\ref{eq:discounted return}) when actions are generated from that policy.

\subsubsection{Policy Gradient Methods}
With \textit{policy gradient methods} [Sutton et~al., 1999, Peters and Schaal, 2006, Silver et~al., 2014, Lillicrap et~al., 2015, Mnih et~al., 2016, Duan et~al., 2016], the policy is explicitly parameterized, often as a deep feed-forward neural network, $\pi=\pi_\theta$. Here $\theta$ is a vector of adjustable network parameters.  A simple (and widely used) stochastic policy is constructed as a diagonal Gaussian distribution, where the actions are uncorrelated as all non-diagonal entries of the covariance matrix are zero
\begin{equation}\label{eq:GaussianPolicy}
\begin{aligned}
a &= \mu_\theta(s) + \sigma_\theta(s) \odot z \\
z &\sim \mathcal{N}(0,I)
\end{aligned}
\end{equation}
Here $z \in\mathbb{R}^{\dim(a)}$ is sampled from the standard normal distribution, while $\mu_\theta\in\mathbb{R}^{\dim(a)}$ and $\sigma_\theta \in\mathbb{R}^{\dim(a)}$ are the parameterized state-dependent mean and standard deviation, respectively, so that $\pi_\theta(a|s)$ represents a diagonal Gaussian density. As an example, in the one-dimensional case we have $\pi_\theta(a|s)=\frac{1}{\sqrt{2\pi}\sigma_\theta(s)} \exp\left(-\frac{(a-\mu_\theta(s))^2}{2\sigma_\theta(s)^2}\right)$.  Policy gradient methods seek to maximize the expected return over actions and states induced from the policy and state transition distributions, respectively,
\begin{equation}
\begin{aligned}
\min_\theta& ~~L(\theta) = \E_{\substack{a\sim\pi_\theta(\cdot|s) \\ s'\sim p(\cdot|s,a)}} \left[R_0(\pi_\theta)\right]
\end{aligned}
\end{equation}
The policy parameters are adjusted through stochastic gradient ascent of the form
\begin{equation}\label{eq:SGD}
\theta_{k+1} = \theta_k + \alpha_{LR} \hat{g}(\theta)|_{\theta_k}
\end{equation}
where $\alpha_{LR}$ is the learning rate.  Here $\hat g(\theta) \triangleq \nabla_\theta L^{PG}(\theta) \approx \nabla_\theta L(\theta)$ is a sample average approximation of the policy gradient, obtained by differentiating a surrogate objective function $L^{PG}(\theta)$.  This surrogate objective function is created over a finite batch of samples in an algorithm that alternates between sampling and optimization.  Several forms of the surrogate objective are available [Schulman et~al., 2015], one of the simplest being\footnote{For taking the derivative of $L^{PG}$, $\log \pi_\theta(a|s)$ can be analtyically expressed in terms of $\mu_\theta(s)$ and $\sigma_\theta(s)$ and the dependence of these terms on $\theta$ is computed through automatic differentiation.}
\begin{equation}\label{eq:L_pg}
\begin{aligned}
&=  \frac{1}{N}\sum_{i=1}^{N}\left(\sum_{h=0}^H  R_h^{(i)} \log \pi_\theta(a_h^{(i)}|s_h^{(i)}) \right)
\end{aligned}
\end{equation}
where the $s_h^{(i)}, a_h^{(i)}, R_h^{(i)}$ notation indicates states, actions and returns  from the $i^{th}$ of $N$ trajectory samples.  In section \ref{section:anti-windup} we formulate an inequality constrained optimization problem using a modified surrogate objective function.

\section{Problem Formulation}\label{section: ProblemFormulation}
In this section we formulate trajectory generation as a Deep RL problem endowed with a high level action space.  Trajectory generation is often formulated as an optimal control problem and then converted into a parameter optimization problem through direct transcription methods [Kelly, 2017, Schlossman et~al., 2021].  This produces a problem formulation that is compatible with modern nonlinear programming techniques [Biegler and Zavala, 2009]. For many problems this method works well [Schlossman et~al., 2021, Patterson and Rao, 2014], but long solution times and convergence stability can limit real-time application.

In this work we develop an alternative method based on the notion of a continuously parameterized motion primitive [Ijspeert et~al., 2002, Dever et~al., 2006, Paranjape et~al., 2015, Camci and Kayacan, 2019] and show the method is compatible with Deep RL.  We emphasize the novel use of a variable duration action, where each action-duration is uniquely chosen at every decision point along the horizon.  Our proposed method improves planning speed and stability by moving the computational burden off-line.  Difficult scenarios are worked out off-line through randomized exploration of the state-space, and the best trajectory formulations are memorized by the policy network.

\subsection{High-Level Action Space (HLAS) Method}\label{section:HLAS}
Consider a system with dynamics as described by (\ref{eq:dynamics}) and a subset of the state $y=Cx\in\mathbb{R}^\ell$ corresponding to the system output.  We define actions as sub-trajectories of variable duration.  Rather than requiring an action decision to be made at every time step, a constrained input sequence executes over multiple time steps and action decisions are only required between sub-trajectories.  The sub-trajectory at action-step $h$, denoted as $z_h$, can represent a range of quantities such as:
\begin{enumerate}
\item A control input function, $z_h(t)= u(t)$
\item A desired dynamics function, $z_h(t)=\dot y^{des}(t)$.  Here it is assumed the control can be designed so the output dynamics follow the desired dynamics function.  A tracking controller for this task is discussed in Section \ref{section:NumericalExp_TrackingController}.
\item A desired system output function, $z_h(t)=y^{des}(t)$.  Here it is assumed the control can be designed with closed-loop feedback so the system output follows the desired output.
\end{enumerate}
In each case above $t_{h-1} \leq t \leq t_{h}, ~~t_{h} = \sum_{i=1}^{h}\tau_i$, and $\tau_h$ is the duration of action-step $h$ for $h=1,\dots, H$.  Here $H$ is the terminal step, a random variable identical to that described in Section \ref{section: deep rl}.  Each sub-trajectory is a time-varying function parameterized as a polynomial of degree $p \geq 0$
\begin{subequations} \label{eq:z}
\begin{align}
z(t) &= z_h\left(\frac{t-t_{h-1}}{\tau_h}\right), ~~t_{h-1} \leq t \leq t_{h}, ~~t_{h} = \sum_{i=1}^{h}\tau_i\\
z_h(t') &= \sum_{k=0}^{p} c_{h,k+1} t'^k, ~~0\leq t' \leq 1
%z_{h+1}(0)&=z_{h}(1) \label{eq:z_partc}
\end{align}
\end{subequations}
for $h=1,\dots,H$ and $t_0=0$.  The polynomial degree, $p \geq 0$, is a design choice made before the problem is solved.  $p+1$ nodes are required to uniquely specify a $p^{th}$ order polynomial.  Nodes $n_{h,1}, \dots, n_{h,p+1}$ are evenly placed along each sub-horizon $\tau_h$ with node $n_{h,1}$ placed at the start of the sub-horizon and $n_{h,p+1}$ placed at the end of the sub-horizon.  The $p+1$ coefficients $c_{h,1},\dots, c_{h,p+1}$ for each step $h$ along the horizon are designed through an appropriate curve fitting interpolation.
Decision variables are:
\begin{enumerate}
\item the $p+1$ nodes $n_{h,1},\dots, n_{h,p+1}$ for each step $h=1,\dots,H$ along the horizon, and
\item the \textit{action duration}, $\tau_h$, over which the polynomial profile is to be applied
\end{enumerate}
Figure \ref{fig:z visual} visually describes the formation of $z(t)$ with variable length segments $z_h$.
\begin{figure}[h!]
\centering
\begin{minipage}{0.65\textwidth}
\includegraphics [trim = 0mm 0mm 0mm 0mm, clip,width=.99\textwidth]{z_profile.png}
\end{minipage}
\captionsetup{format=hang}
\captionsetup{width=.9\textwidth}
\caption{Example piecewise polynomial construction ($p=3$) for $z(t)$ for a horizon with three segments, $H$=3.  $t_{h} = \sum_{i=1}^{h}\tau_i$.}
\label{fig:z visual}
\end{figure}
As evident in Fig. \ref{fig:z visual} discontinuities in $z$ are permitted in the most general form.  However, during implementation an additional constraint of the form $z_{h+1}(0)=z_{h}(1)$ can be added to ensure continuity across sub-trajectories.

In the case where $z = \dot y^{des}$, for a $p^{th}$ order desired dynamics function $z(t)$ we are effectively defining piecewise polynomial desired system output profiles, $y^{des}(t)$, of degree $p+1$,
\begin{equation}\label{eq:y^des}
\begin{aligned}
&= y_{h-1}^{des}(t_{h-1}) + \int_{0}^{(t-t_{h-1}) / \tau_h} z_h(s) ds \\
&= y_{h-1}^{des}(t_{h-1}) + \sum_{k=0}^{p} \frac{1}{(k+1)} c_{h,k+1} \left(\frac{t-t_{h-1}}{\tau_h} \right)^{k+1} \\
t_{h-1} &\leq t \leq t_h ~~,~~ t_{h} = \sum_{i=1}^{h}\tau_i
\end{aligned}
\end{equation}

\subsection{Formulation as a Deep RL Problem} \label{section: HLAS Deep RL Formulation}
\subsubsection{Rewards}
The reward signal of the discounted return (\ref{eq:discounted return}) is defined as
\begin{equation}\label{eq:HLAS Deep RL reward}
r_h =
\begin{cases}
\Phi(x_h) + L_h + c_0\hat\Psi(x_h)^{-1} &\text{ if $h=H$} \\
L_h &\text{ otherwise}
\end{cases}
\end{equation}
\begin{equation} \label{eq:HLAS Deep RL reward performance index}
L_h = \int_{t_{h-1}}^{t_h} L(x(t),u(t))dt,~~t_h=\sum_{i=1}^{h}\tau_i
\end{equation}
is the performance index accumulated along action step $h$, and $\Phi(x_h)$ is the same horizon reward function in (\ref{eq: J}).  Here we have relaxed the terminal constraint function $\Psi=0$ from section \ref{section:Optimal Control} into a reward of the form $c_0\hat\Psi(x_h)^{-1}$ where $\hat\Psi(x_h)$ is a measure of distance from $\Psi(x_h)=0$.
\subsubsection{Action Space}\label{section: ProblemFormulation_actionspace}
The action space of the RL agent is defined as $\mathcal{A}=(\tau, n_{1},\dots,n_{p+1})$ as described in section \ref{section:HLAS}.  At each action step the agent chooses the time duration $\tau$ of the step and the location of the $p+1$ evenly spaced node points $n_i$ as shown in Fig \ref{fig:z visual}.  Lower and upper limits are set on the action duration: $\tau_{min} \leq \tau \leq \tau_{max}$.  Additionally, lower and upper limits are placed on each node point so that each $z_h$ remains bounded:  $z_{min} \leq n_{i} \leq z_{max}, ~i={1,\dots,p+1}$.
\subsubsection{Constraints}
The box form control input constraints, $u\in U$, are enforced directly in the environment through clamping. The control inputs are computed based on the actions and the tracking controller.  For the state constraint, $x\in X$, and path constraints, $g_j(x,u)\leq 0, j=1,\dots,J$, we employ a simple strategy in which the episode is ended and the agent receives no further reward if the agent violates any of these constraints.  Over the course of learning, the agent learns to avoid regions of infeasibility.

\subsubsection{HLAS Diagram}
The figure below shows the HLAS Deep RL agent-environment interaction loop. At every action-step $h$ the agent generates an action $a_h=(\tau_h, \{n_{h,i}\})$ which specifies the profile of $z_h(t)$ along $t_{h-1} \leq t \leq t_h$ where $t_h=\sum_{i=1}^{h}\tau_i$.  The number of timesteps $dt$ between action-steps $h$ is determined by $\tau_h$.
The tracking controller converts $y^{des}$ or $\dot y^{des}$ into a control input $u$ depending on the definition of $z$ (when $z_h(t)=u(t)$ the tracking controller is the identity function).
The environment (which is comprised of the $z_h$ profile generator, the tracking controller, the dynamics, and the reward definition) provides the agent with an updated reward and state observation every action-step.

\begin{figure}[h!]
\centering
\begin{minipage}{0.99\textwidth}
\includegraphics [trim = 0mm 0mm 0mm 0mm, clip,width=.99\textwidth]{HLAS_diagram.png}
\end{minipage}
\captionsetup{format=hang}
\captionsetup{width=.9\textwidth}
\caption{Implementation of the HLAS Deep RL formulation.}
\label{fig:HLAS_diagram}
\end{figure}

\subsection{Analysis}\label{section:Analysis}
\subsubsection{Accuracy of Approximation}
The method described in Section \ref{section:HLAS} is a way to parameterize a trajectory generation problem for Deep RL.  It can also be viewed as a way to approximate a solution to the optimal control problem in Section \ref{section:Optimal Control}.  In Section \ref{section:DebrisAvoidance} we will use the methods described in Section III to perform trajectory planning on a problem which is challenging for traditional optimal control solvers. The following shows that the best approximation error grows at a bounded linear rate along the horizon in the case where $z=\dot y^{des}$.

\textbf{Lemma} (Accuracy of HLAS approximation, for $z=\dot y^{des}$).
Define $x^*:[t_0,t_H]\rightarrow\mathbb{R}^n$, $u^*:[t_0,t_H]\rightarrow\mathbb{R}^m$ as a solution to (\ref{eq:dynamics}) which is optimal with respect to (\ref{eq: J}).  Define $\dot {\hat x}(t)=z^*(t)$, where $z^*(t)$ is a piecewise polynomial as described in (\ref{eq:z}) over $t\in [t_0,t_H]$ and best approximates $f(x^*(t),u^*(t))$ over $t\in [t_0,t_H]$ in some sense.  Assume $\hat x(t_0) = x^*(t_0)$ is a specified initial condition.  Furthermore assume $\Vert z^*(t) - f(x^*(t), u^*(t))\Vert_\infty \leq m_1 ~\forall t \in [t_0,t_H]$, for some $m_1>0$.  Then $\Vert \hat x(t) - x^*(t) \Vert_\infty \leq m_1[t-t_0] \forall t \in[t_0,t_H]$.

\textit{Proof}.  $\hat x(t) - x^*(t) = \int_{t_0}^{t} z^*(\tau) d\tau + \hat x(t_0) - \int_{t_0}^{t} f(x^*(\tau),u^*(\tau)) d\tau - x^*(t_0) = \int_{t_0}^{t} z^*(\tau) - f(x^*(\tau),u^*(\tau)) d\tau $.  We then have $\Vert \hat x(t) - x^*(t) \Vert_\infty = \int_{t_0}^{t} \Vert z^*(\tau) - f(x^*(\tau),u^*(\tau)) \Vert_\infty d\tau \leq \int_{t_0}^{t} m_1 d\tau = m_1 [t-t_0]$.

\textbf{Remark} (Bounds on $m_1$).  Bound each component of $f(x^*(t),u^*(t))$ as $L_i(t_h) = \inf\limits_{t_h \leq t \leq t_h+\tau_h} f_i(x^*(t),u^*(t))$,
$U_i(t_h) =  \sup\limits_{t_h \leq t \leq t_h+\tau_h} f_i(x^*(t),u^*(t))$,
where $\tau_h$ is the action duration as described in Section \ref{section:HLAS}.  Set $z_i(t) = \frac{U_i(t_h)-L_i(t_h)}{2}$ as a suboptimal approximation to $f_i(x^*(t),u^*(t))$ along $t \in [t_h,t_h+\tau_h]$.  We then have $\vert z_i(t)-f_i^*(t) \vert  \leq \frac{\Delta_i}{2}$, where $\Delta_i=\sup_{h} (U_i(t_h) - L_i(t_h))$.  We can therefore take $m_1=\sup_i\frac{\Delta_i}{2}$.

The bound $\Delta_i$ will generally decrease as $\tau_h$ is decreased by the definitions of $L_i(t_h)$ and $U_i(t_h)$.  Furthermore, by the Weierstrass Approximation Theorem, $f_i(x^*, u^*)$ can be approximated to arbitrary accuracy by a polynomial of sufficient degree.  Therefore, $\Delta_i$ will generally decrease as the polynomial degree $p$, as described in Section \ref{section:HLAS}, is increased.

\subsubsection{Effect on Environment Exploration}
The variable action duration $\tau_h$ has interesting implications in terms of environment exploration during the RL training process.  We now examine the effect of discounting on the action duration for a specific reward structure, and show longer action durations are preferred. Longer action durations induce larger movements through the state space between gradient updates, promoting environment exploration and ultimately producing improved policies.  This effect is empirically demonstrated in Section \ref{section:shuttle_ablation_studies}.

\textbf{Definition} (Strictly Episodic Reward).  A reward that is guaranteed to occur every episode, but only at a terminal state: $r_h=0 ~\forall~ h \neq H$, $r_H > 0$, where $H$ is the horizon length described in Section \ref{section: deep rl}.

When positive and strictly episodic rewards are used in conjunction with the HLAS Deep RL formulation, a discount factor $0 < \gamma < 1$ encourages the agent to favor larger action durations $\tau_h$.  The discounted episode return is $R=\sum_{h=1}^{H} \gamma^{h-1}r_h = \gamma^{H-1} r_H$, where the term $\gamma^{H-1}$ is the amount by which the episodic reward $r_H$ is discounted.  Consider two horizon lengths $H_1, H_2$ where $H_1<H_2$.  If $r_{H_1} = r_{H_2}$, horizon length $H_1$ will be preferred by the agent as this produces a larger discounted episode return.  As a result, the agent will attempt to keep $H$ small by choosing larger action durations $\tau_h$.

\subsection{Deep RL Training Setup}\label{section:RLTrainingSetup}
A range of algorithms can be used to determine an optimal policy through Deep RL training, including off-policy algorithms such as Soft Actor Critic (SAC) [Haarnoja et~al., 2018], Deep Deterministic Policy Gradient (DDPG) [Silver et~al., 2014], and Twin Delayed Deep Deterministic Policy (TD3) Gradients [Fujimoto et~al., 2018], and on-policy algorithms such as Trust Region Policy Optimization (TRPO) [Schulman et~al., 2015], Asynchronous Advantage Actor Critic (A3C) [Mnih et~al., 2016] and Proximal Policy Optimization (PPO) [Schulman et~al., 2017].  While off-policy methods tend to be more sample efficient (since they can reuse old data through the use of a replay buffer), on-policy methods directly optimize the objective and favor stability over sample efficiency [Achiam, 2018].  TRPO is a second-order algorithm which prevents large, sometimes destructive updates to the policy parameters by restricting the amount by which an update is allowed to change the policy [Gangwani and Peng, 2018].  PPO is a first order approximation to TRPO and restricts the size of policy update via a KL divergence penalty or clipping mechanism [Gangwani and Peng, 2018].  PPO is one of the most popular Deep RL methods, achieving state-of-the-art performance across a wide range of challenging tasks [Wang et~al., 2019, Schulman et~al., 2017].  In this work we use PPO for all Deep RL training.

Our implementation of PPO uses the following surrogate objective function,
\begin{equation}\label{eq:PPOSurrogateObj}
\begin{aligned}
L^{PPO}(\theta) &= \hat\E_{h}\left[ L_h^{\text{clip}}(\theta) - c_1 L_h^{VF}(\theta) + c_2 S[\pi_\theta](s_h) + c_3 \mu_\theta^+(s_h) \right]
\end{aligned}
\end{equation}
Here $\hat \E_h[f(s_h,a_h)]$ indicates the sample average approximation to $\E_{\substack{s\sim \rho^{\pi_\theta}\\a\sim\pi_\theta(\cdot|s)}} \left[f(s,a)\right]$ where $\rho^{\pi_\theta}$ is the distribution of states induced by the policy.  The terms $L_h^{\text{clip}}(\theta)$, $L_h^{VF}(\theta)$ and $S[\pi_\theta](s_h)$ are the standard PPO-clip, value function and entropy exploration bonus terms, respectively, taken directly from [Schulman et~al., 2017]. Only the first two terms are required for PPO, we include the entropy term as it has been shown to improve exploration by discouraging premature convergence to suboptimal policies [Williams and Peng, 1991, Mnih et~al., 2016].  The last term in the objective is a new contribution in this work for bounding the learned action distribution, given by
\begin{equation}\label{eq:AntiWindupLossTerm}
\mu_\theta^+(s_h) = \hat\E_{h}\left[\sum_{j=1}^{\dim(a)} \max\left(\left|\mu_{j,\theta}(s_h)\right|-(1-\epsilon), 0\right)^2 \right]
\end{equation}
as described in Section \ref{section:anti-windup}.  Coefficients $c_1$ and $c_2$ are hyperparameters chosen heuristically, while coefficient $c_3$ is automatically adjusted according to the update scheme described in Section \ref{section:anti-windup}.

The PPO algorithm alternates between data collection and policy optimization.  During data collection the policy $\pi_\theta$ is executed on multiple instances of the environment in parallel [Mnih et~al., 2016, Schulman et~al., 2017], performing a fixed number of action-steps on each parallel environment.  During optimization stochastic gradient descent is performed by cycling through minibatches of the collected data to obtain sampled estimates of $\nabla_\theta L^{PPO}(\theta)$.  The final output of DeepRL training is a fully trained policy network $\pi_\theta(a|s)$.  In this work we use the Stable Baselines3 codebase [Raffin et~al., 2019] for implementation of the PPO algorithm.

After training is complete, the trained policy can be used to produce actions $a$ from any state $s$.  Referring to Eq. (\ref{eq:GaussianPolicy}), rather than sampling from the full stochastic policy we extract actions deterministicly from the action-mean portion of the policy network: $a=\mu_\theta(s)$.

\section{Preventing Action Windup in Stochastic Policy Gradients}
Almost all physical systems have action limits due to hardware limitations.  These can be limits on control surface deflections or actuator force magnitudes.  When an agent exceeds these bounds during learning, problems can arise.  In this section we introduce a novel method for bounding the learned action distribution in policy gradient methods.  We focus on settings where the action the environment can accept is limited with box type bounds, producing a ``clipped'' action.  We demonstrate the effectiveness of this method in Section \ref{section:TraversedLatitude}.

The issues associated with bounded actions in Deep RL policy gradients has been explored previously.  One approach is to use the clipped action, referred to here as $\hat a_h$, directly in the policy gradient estimate in place of the unclipped action $a_h$.  Unfortunately this has been shown to introduce significant bias into the policy gradient estimation process [Chou et~al., 2017] by effectively corrupting the empirical average (\ref{eq:L_pg}).  The clipped action policy gradient method [Fujita and Maeda, 2018] attempts to alleviate this bias issue by replacing the gradient with calculations based on the cumulative distribution function (CDF).  However, computing the CDF can be non-trivial (such as with Gaussian distributions where no closed form exists).  The work of [Haarnoja et~al., 2018] employs the use of a $\tanh$ \textit{squashing function} to transform an unbounded Gaussian distribution to a bounded distribution as $\hat a_h = \tanh(a_h)$, where $a_h$ is the action produced by the (Gaussian) policy (\ref{eq:GaussianPolicy}).  However, the $\tanh$ transformation does not solve the ``windup issue,'' as the input to the $\tanh$ function can grow without bound with no further effect on the output.  The work of [Chou et~al., 2017] investigates the applicability of the Beta distribution, which has bounded support, to represent the stochastic policy.  One drawback of the Beta distribution is that the policy must be nearly deterministic for near-limit action values to be chosen with high probability [Fujita and Maeda, 2018].

\subsection{The Windup Problem}
Consider the typical agent-environment reinforcement learning training process [Sutton and Barto, 2018], shown in Fig. \ref{fig:agent-env-interaction} for box-type action constraints.
\begin{figure}[h!]
\centering
\begin{minipage}{0.65\textwidth}
\includegraphics [trim = 0mm 0mm 0mm 0mm, clip,width=.99\textwidth]{agent-env.png}
\end{minipage}
\captionsetup{format=hang}
\captionsetup{width=.9\textwidth}
\caption{The agent–environment interaction in reinforcement learning, with action limits explicitly shown separated from the environment.}
\label{fig:agent-env-interaction}
\end{figure}
Here we have separated the action limits from the environment for clarity in our explanation.  This interaction can be viewed as a typical closed-loop feedback control process [\AA}str{\"o}m and Murray, 2021] where the agent serves the purpose of a typical controller and the environment consists of the system being controlled, including all measurements.  The presence of the action limit in the forward path of the loop can be problematic.  For systems which operate over a wide range of conditions, the action delivered by the agent may very well reach a physical bound represented here by the action limit.  When this happens,  the feedback loop is broken and the process runs in open-loop, independent of the feedback from the environment [\AA}str{\"o}m and Murray, 2021].  For feedback controllers which possess an integrator process (i.e. for controllers which have a form of memory) the consequence can be long and undesirable transient behavior, which is sometimes destabilizing.  In the feedback control literature this well-studied phenomenon is referred to as \textit{integrator windup}.  \textcolor{black}{In the context of Deep RL, windup can occur in the action mean $\mu_\theta$ as it approaches the action limit.  Actions are sampled from the Gaussian policy, $a_h \sim \pi_\theta(\cdot|s_h)$, and all samples outside the limit will have the same effect on the environment as actions at the limit.  If a favorable reward $r_{h+1}$ is consistently received at the action limit at some point during training, then the aggregate effect of a policy gradient method will push the action mean $\mu_\theta$ towards and even well beyond the limit.  This can be problematic if the probability of sampling an action within the action limit becomes sufficiently low (due to the drifting $\mu_\theta$), since at this point the agent has effectively stopped exploring the environment.}

In feedback control literature, a well-established method of \textit{anti-windup} involves making an adjustment to the input of the feedback controller based on the difference between the pre-limited and post-limited actions [\AA}str{\"o}m and Murray, 2021].  We now introduce an analogous method for policy gradient methods in Deep RL.
\subsection{Inequality Constrained Policy Optimization for Anti-Windup} \label{section:anti-windup}
\subsubsection{KKT Conditions}
The policy optimization problem of minimizing $L^{PG}(\theta)$ in (\ref{eq:L_pg}) is reformulated as the following inequality constrained optimization problem
\begin{subequations}\label{eq:ConstrainedSurrogateOpt}
\begin{align}
\min_\theta&~~~ L^{PG}(\theta) \\
\text{subject to}&~~~ \hat\E_{h}\left[\mu_{j,\theta}^2(s_h)-1\right] \leq 0, ~~j=1,\dots,\dim(a) \label{eq:Inequality}
\end{align}
\end{subequations}
where $\hat\E_h[\dots]$ indicates the sample average with respect to quantities indexed by $h$, and $\mu_{j,\theta}$ is the mean of action $j$.  The inequality constrained optimization (\ref{eq:ConstrainedSurrogateOpt}) seeks to minimize the standard surrogate objective function $L^{PG}(\theta)$ while maintaining the mean of each action dimension bounded between $[-1,1]$.  We construct the following Lagrangian function and associated Karush–Kuhn–Tucker (KKT) first order necessary conditions,
\begin{subequations}
\begin{align}
L^{CONS}(\theta) &= L^{PG}(\theta) + \sum_{j}c_j  \hat\E_{h}\left[\mu_{j,\theta}^2(s_h)-1\right] \label{eq:KKT Lagrangian} \\
\left.\nabla_\theta L^{CONS}(\theta)\right|_{\theta^*,c_j^*} &= 0 \label{eq:KKT stationary point} \\
c_j^* \hat\E_{h}\left[\mu_{j,\theta^*}^2(s_h)-1\right] &=0 \label{eq:KKT complimentary} \\
c_j^* & \geq 0 \label{eq:KKT positive cj}
\end{align}
\end{subequations}
for $j=1,\dots,\dim(a)$. We next devise a scheme to approximately solve for the conditions (\ref{eq:KKT stationary point})-(\ref{eq:KKT positive cj}) under stochastic gradient ascent optimization.

\subsubsection{Approximately Solving for the KKT Conditions} \label{section:anti-windup_approxKKT}
We adopt a penalty function approach from constrained optimization [Bertsekas, 2016, Chong and Zak, 2008] in which case the KKT conditions are approximately satisfied under stochastic gradient ascent optimization.  The main idea is the inequality constraint (\ref{eq:Inequality}) is replaced with a positive penalty function which resembles the inequality, and penalty coefficients are adapted sequentially.  In the surrogate objective function below, we introduce a penalty function of the form $P(\mu)=\max\left(\left|\mu\right|-1, 0\right)^2$,
\begin{equation}\label{eq:KKT_Approx Lagrangian}
L^{CONS}(\theta) = L^{PG}(\theta) + \sum_{j}c_j  \hat\E_{h}\left[ \max\left(\left|\mu_{j,\theta}(s_h)\right|-1, 0\right)^2 \right]
\end{equation}
This penalty function has several desirable properties: it has continuous derivatives with respect to $\mu$ at $\mu=+/-1$, $P(\mu)>= 0$ for all $\mu$, and $P(\mu)=0$ if and only if $\mu\in[-1,1]$.  Condition (\ref{eq:KKT complimentary}) is the so-called complimentary slackness condition, stating that $c_j^*$ must be zero whenever the inequality in (\ref{eq:Inequality}) is inactive.  We relax (\ref{eq:KKT complimentary}) as
\begin{equation}\label{eq: KKT Approx complimentary}
c_j^* \hat\E_{h}\left[ \max\left(\left|\mu_{j,\theta^*}(s_h)\right|-1, 0\right)^2 \right] = 0, ~~j=1,\dots,\dim(a)
\end{equation}
noting that any ($\theta^*, c_j^*$) which satisfy (\ref{eq:KKT complimentary}) automatically satisfy (\ref{eq: KKT Approx complimentary}).  The conditions specified in (\ref{eq:KKT stationary point}) and (\ref{eq:KKT positive cj}) place further restrictions on the choice of $\theta^*$ and $c_j^*$.  Specifically (\ref{eq:KKT stationary point}) states that the optimal values must be concurrently chosen to produce a stationary point in the surrogate objective.  For this we employ two processes.  First, standard stochastic gradient ascent is applied to the surrogate objective (\ref{eq:KKT_Approx Lagrangian}) to adjust parameters $\theta$.  Second, at each step of stochastic gradient ascent, the coefficients $c_j$ are adapted according to a method developed in [Schulman et~al., 2017]\footnote{The method developed in [Schulman et~al., 2017] is for adapting a coefficient in KL divergence targeting.  We adopt the method nearly verbatim only replacing the definition of $d_j$ as appropriate.  The parameters 1.5 and 2 were taken directly.}:
\begin{itemize}
\item Define some target value $d_{\text{tar}}$
\item Compute $d_j = \hat\E_{h}\left[ \max\left(\left|\mu_{j,\theta}(s_h)\right|-1, 0\right)^2 \right]$
\subitem - If $d_j < d_{\text{tar}}/1.5$, ~$c_j \leftarrow c_j/2$
\subitem - If $d_j > d_{\text{tar}}*1.5$, ~$c_j \leftarrow c_j*2$
\end{itemize}
Note that this adaptation process automatically satisfies (\ref{eq:KKT positive cj}), assuming each $c_j$ is initialized with some positive value.  Examining our composite surrogate objective (\ref{eq:KKT_Approx Lagrangian}), the second term places external pressure on the network optimization process to find values $\theta$ which, on average, keep the action means bounded within the feasible range of the action space.

\subsubsection{Implementation Note}
In actual implementation we slightly modify the summation term in (\ref{eq:KKT_Approx Lagrangian}) as
$c  \hat\E_{h}\left[\sum_{j} \max\left(\left|\mu_{j,\theta}(s_h)\right|-(1-\epsilon), 0\right)^2 \right]$.
For one, we have made the simplification of summing over all action dimensions inside the empirical average.  As a result, we only need to adapt one coefficient $c$\footnote{For some applications there may be benefit in adapting all coefficients $c_j$, $j=1,\dots,\dim(a)$, but in this work we found the simplification to a single coefficient acting upon a summation of terms to be sufficient.} using $d=\hat\E_{h}\left[\sum_{j} \max\left(\left|\mu_{j,\theta}(s_h)\right|-(1-\epsilon), 0\right)^2 \right]$ in place of $d_j$ in the update scheme above.  Second, we use $(1-\epsilon)$ with $\epsilon\approx 0.1$ inside the max function and set $d_\text{tar}=\epsilon^2$.  This allows us to target a positive value of $d_\text{tar}$ in the update scheme above while seeking $\hat\E_h\left[\sum_{j=1}^{\dim(a)}\left|\mu_{j,\theta}(s_h)\right|\right] \leq 1$.

\section{Numerical Experiments with the Space Shuttle Reentry Problem}\label{section:NumericalExp}
In this section we evaluate the performance of our approach using the well-studied space shuttle reentry problem [Betts, 2010].  We benchmark our approach against a state-of-the-art nonlinear programming solution, and a straightforward implementation of Deep RL in which action decisions are made by the RL agent every simulation time step.  A trajectory planning problem is also studied where the reentry must be planned over a wide range of initial conditions, avoiding a quasi-static group of obstacles within the map.

Deep RL training is performed using the PPO algorithm as described in Section \ref{section:RLTrainingSetup}.  We use six\footnote{Training was performed on a machine with a six core CPU.} parallel instances of the environment with 4096 action-steps per environment, and minibatch sizes of 128 action-steps during the optimization phase unless specified otherwise.  We use a network with two shared hidden layers of size 256, followed by an additional hidden layer of size 256 which is evenly split between the policy and value function approximation [Raffin et~al., 2019, Mnih et~al., 2016].  The rectified linear unit (ReLU) activation function was used in each layer, various hyperparameters used throughout training are listed in the following sections.  The shuttle dynamics Eq. (\ref{eq:sphericalDynamics}) are propagated along the horizon using a fourth order Runge-Kutta numerical integration scheme with a time step of two seconds.  The Deep RL training is performed on a laptop equipped with a six-core 2.60GHz Intel i7 CPU and an NVIDIA Quadro P3200 GPU.

\subsection{Shuttle Vehicle Description}\label{section:NumericalExp_shuttle_descr}
We briefly describe the problem setup, the full details can be found in [Betts, 2010].  The problem is concerned with a non-thrusted space reentry vehicle starting at a specified initial condition, gliding towards a specified terminal condition.  The motion of the vehicle is described by the following nonlinear dynamics,
\begin{equation}
\label{eq:sphericalDynamics}
\begin{aligned}
&\dot h =v \sin\gamma , & \dot v &= -\frac{D}{m}-g\sin\gamma \\
& \dot \theta = \frac{v}{h+Re}\cos\gamma\sin\psi/\cos\phi , & \dot \phi &= \frac{v}{h+Re}\cos\gamma\cos\psi \\
& \dot\gamma = \frac{L\cos\sigma}{mv} + \(\frac{v}{h+Re}-\frac{g}{v}\)\cos\gamma , &
\dot\psi &= \frac{L\sin\sigma }{mv\cos\gamma}+\frac{v}{h+Re}\cos\gamma\sin\psi \frac{\sin\phi}{\cos\phi} \\
& \dot \alpha  = \frac{1}{\tau_{\alpha}}\left[\alpha_{cmd}-\alpha \right],
& \dot \sigma  &= \frac{1}{\tau_{\sigma}}\left[\sigma_{cmd}-\sigma \right]
\end{aligned}
\end{equation}
where $R_e$ is radius of the Earth [m], $m$ is vehicle mass [kg], $h$ is altitude [$m$], $\theta$ is longitude [rad], $\phi$ is latitude [rad], $v$ is velocity [m/s], $\gamma$ and $\psi$ are vertical and horizontal flight path angles [rad], respectively, $\alpha$ is angle of attack [rad], and $\sigma$ is bank angle [rad].  The last two equations represent a first order dynamics model of the angle of attack and bank angle control loops, where $\tau_{\alpha}=\tau_{\sigma}=1$ second.  The following path constraints are present
\begin{equation} \label{eq:shuttle state path constraints}
\begin{aligned}
&h \geq 20000 \text{ [m]}, ~~
v \geq 600 \text{ [m/s]}, ~~
-20 \text{ [deg]} \leq ~\gamma \leq 20\text{ [deg]}
\end{aligned}
\end{equation}
The aerodynamic and gravitational forces are computed as
\begin{equation}
\label{eq:aeroForces}
L=\frac{1}{2}C_L S \rho v^2~,~D=\frac{1}{2}C_D S \rho v^2~,~\rho=\rho_0\exp(-h/H_0)~,~g=\frac{\mu}{(h+R_e)^2}
\end{equation}
where $S$ is surface area [m$^2$] and $H_0$, $\rho_0$, $\mu$ are Earth-specific constants. The variables $C_D=b_0+b_1\hat\alpha+b_2\hat\alpha^2$ and $C_L=a_0+a_1\hat\alpha$ are the aerodynamic drag and lift coefficients, respectively.  The state and control vectors are defined as $x=[h,v,\theta,\phi,\gamma,\psi, \alpha, \sigma]\Tr$ and $u=[\alpha_{cmd}, \sigma_{cmd}]$, respectively.  The following control constraints are imposed with clipping inside the environment: $-45 \text{ [deg]} \leq \alpha_{cmd} \leq 45\text{ [deg]}$, $-89\text{ [deg]} \leq \sigma_{cmd} \leq 89\text{ [deg]}$.

\subsection{Tracking Controller} \label{section:NumericalExp_TrackingController}
As discussed in Section \ref{section:HLAS}, a tracking controller is required in the case where $z_h(t)$ represents $\dot y^{des}(t)$ or $y^{des}(t)$.  Here we describe such a controller based on the dynamic inversion (DI) principle.  In typical DI controller synthesis the system is first linearized so the control input appears as an affine term in the system dynamics [Tip{\`a}n et~al., 2020].  Here we take an alternative approach which does not require any linearization.  We define the system output $y=\left[\gamma, \psi\right]\Tr$ and desired dynamics $\dot y^{des} = \left[\dot\gamma^{des}, \dot\psi^{des}\right]\Tr$.
\def \fpaDot {\dot y^{des}}
\def \bankCmd {\sigma_{cmd}}
\def \alphaCmd {\alpha_{cmd}}
As a gliding vehicle, the Space Shuttle must achieve its high level goals (i.e. the flight path angle rate commands) through adjustment of its attitude with respect to the ``wind vector''.
Since no wind is modeled, the ``wind vector'' is the velocity vector.
The wind-relative attitude is parameterized by angle of attack $ \alphaCmd $ and bank angle $ \bankCmd $.  We will first construct an acceleration command that is a function of $ \fpaDot $.
Assuming the vehicle has an onboard nominal aerodynamic model, the acceleration command can be used to determine the bank angle and angle of attack necessary to achieve the desired acceleration.

\def \lat {\phi}
\def \lng {\theta}

\def \vpa {\gamma}
\def \hpa {\psi}

\def \Cev {\dcm[E][V]}

To develop this method, we first define the so-called ``velocity frame''.
We define this frame by describing the direction cosine matrix relating the earth-fixed-inertial (ECI) frame to this velocity frame assuming a non-rotating earth ([Regan and Anandakrishnan, 1993] Appendix H.2 and H.3), $\Cev=\dcm[N][V](\vpa,\hpa)\dcm[E][N](\lat,\lng)$.
\def \vE {\stdVec{v}[E][][E]}
By definition, the inertial velocity $ \stdVec{v}[][][E] $ of the vehicle lies completely along the x-axis of the velocity frame such that
\begin{equation}\label{eq:velocityInVFrame}
\rowVec{v,0,0}\trans = \Cev \vE
\end{equation}
Taking the derivative of both sides of \eqref{eq:velocityInVFrame} and substituting \eqref{eq:sphericalDynamics}, we can solve for the inertial acceleration of the vehicle as a function of $ \fpaDot $, $ \dot{v} $, and the current state.
\def \accelV {\stdVec{a}[V][][E]}
\def \Re {R_e}
\begin{equation}\label{key}
\accelV = \colVec{
\dot{v},
\dfrac{v\cos(\gamma)\big(h~\dot{\psi}\cos(\phi) + \Re\dot{\psi}\cos(\phi) - v\cos(\gamma)\sin(\phi)\sin(\psi)\big)}{\cos(\phi)(\Re + h)},
- \dfrac{v\big(\Re\dot{\gamma} + \dot{\gamma}h - v\cos(\gamma)\big)}{\Re + h}
\end{equation}
Only the force due to gravity and the aerodynamic force act on the vehicle.
Applying Newton's second law we find that
\def \FaeroV {\stdVec{f}[V][a]}
\def \gravV {\stdVec{g}[V]}
\begin{equation}\label{eq:NewtonsSecond}
m \accelV	= \FaeroV + m \gravV
\end{equation}
Assuming that gravity acts in the local ``down'' direction, we can define the gravity vector as $\gravV = [{-g \sin(\gamma), 0, g \cos(\gamma)}]\trans$.
Solving \eqref{eq:NewtonsSecond} for $ \FaeroV $ yields the desired aerodynamic force that would achieve our desired acceleration.
\begin{equation}\label{key}
\FaeroV = m(\accelV-\gravV) = \rowVec{f_1, f_2, f_3}\trans
\end{equation}
Using $ \FaeroV $ we can implement a ``bank to turn'' system.
In ``bank to turn'', a bank angle is applied in order to line up the lift vector in the direction of desired acceleration (see [Regan and Anandakrishnan, 1993] chapter 9.4).
No aerodynamic side force is experienced.
This is in contrast to a ``skid to turn'' system which uses both lift and side forces to set the aerodynamic force direction.

The ``wind frame'' is the result of a right-handed rotation about the velocity x-axis through the bank angle.
The lift force acts in the negative direction of the wind-frame z-axis.
We can thus find the bank angle that ``lines up'' the lift force in the desired direction using $\bankCmd = \arctan \bigg(\frac{f_2}{-f_3} \bigg)$.
\def \liftCmd {L_{cmd}}
Next the desired lift force magnitude $ \liftCmd $ is determined $\liftCmd=\sqrt{f_2^2 + f_3^2}$.
\def \cLift {C_{L}}
\def \sRef {S}
As shown in \eqref{eq:aeroForces}, the lift force is a function of the lift coefficient $ \cLift(\alpha) $, air density $ \rho $, speed $ v $, and reference surface area $ \sRef $.
We can use \eqref{eq:aeroForces} to solve $ \alphaCmd  $ such that $\liftCmd = \frac{1}{2}\cLift(\alphaCmd) \sRef \rho v^2$.  In general, the solution to this $L_{cmd}$ may be analytical or may require numerical solutions.
In this work, it can be solved analytically.

\subsection{Problem 1: Maximizing Traversed Latitude}\label{section:TraversedLatitude}

Consider the following well-studied trajectory generation problem with terminal condition as defined in [Betts, 2010]
\begin{equation}\label{eq: shuttle terminal condition}
\begin{aligned}
h_f&=24384\text{ [m]}, & v_f &=762\text{ [m/s]}, & \gamma_f&= -5\text{ [deg]}
\end{aligned}
\end{equation}
The goal is to maximize traversed latitude $\phi(t_f)$ along the horizon from the following nominal initial conditions
\begin{equation}\label{eq:shuttle_ics}
\begin{aligned}
h_0&=79248\text{ [m]}, & v_0 &=7802\text{ [m/s]}, & \gamma_0&= -1\text{ [deg]}, & \alpha_0&=0 \text{ [deg]} \\
\theta_0&=0 \text{ [deg]}, & \phi_0&=0 \text{ [deg]}, & \psi_0 &= 90\text{ [deg]}, & \sigma_0&=0 \text{ [deg]}
\end{aligned}
\end{equation}
In addition to the path constraints listed in Eq. (\ref{eq:shuttle state path constraints}), we also enforce the following heating constraint taken from [Betts, 2010] at each time step along the horizon
\begin{equation} \label{eq:shuttle_heating_constraint}
\begin{aligned}
&q = 779.67 (c_0 + c_1\hat\alpha + c_2\hat\alpha^2 + c_3\hat\alpha^3) \sqrt{\rho}(3.28084\times 10^{-4} v)^{3.07} \\
&q \leq 80 \text{ [BTU/ft$^2$ - s]}
\end{aligned}
\end{equation}
where $q$ is the aerodynamic heating on the vehicle wing leading edge [$\text{BTU}/\text{ft}^2-s$] and $\hat \alpha = \frac{180}{\pi} \alpha$.  Referring to Eq. (\ref{eq: J}), in this problem there is no running performance measure and the horizon reward, $\Phi$, is simply the latitude achieved at episode termination.  Following the reward formulation Eq. (\ref{eq:HLAS Deep RL reward}), we set
\begin{subequations}
\begin{align}
L &= 0 && \text{(running performance reward)}\\
\Phi &= e^{\phi}\times\mathbf{1}_{\phi<0} + (1+\phi)\times\mathbf{1}_{\phi\geq0} && \text{(horizon reward)} \label{eq:shuttle_perf_reward}\\
\Psi &= \left(\frac{h-h_f}{\bar h}\right)^2
+ \left(\frac{v-v_f}{\bar v}\right)^2
+ \left(\frac{\gamma-\gamma_f}{\bar \gamma}\right)^2 && \text{(terminal ellipsoid)}\\
\hat\Psi^{-1} &= \min\left(1, \Psi^{-1}\right) && \text{(terminal reward)} \label{eq:shuttle_terminal_reward} \\
c_0 &= 5 && \text{(terminal reward weight)}
\end{align}
\end{subequations}
where $\bar h=250$ [m], $\bar v=8$ [m/s] and $\bar\gamma=0.1$ [deg] are scale factors.  Here the performance reward Eq. (\ref{eq:shuttle_perf_reward}) is always positive and increases monotonically with latitude, and the terminal condition Eq. (\ref{eq: shuttle terminal condition}) is cast as a bounded inverse distance reward Eq. (\ref{eq:shuttle_terminal_reward}).  We additionally define the \textit{terminal tolerance} criteria,
\begin{equation}\label{eq:shuttle_terminal_tolerance}
\begin{aligned}
|h-h_f|\leq h_{tol}, ~~~~ |v-v_f| \leq v_{tol}, ~~~~|\gamma-\gamma_f| \leq \gamma_{tol}
\end{aligned}
\end{equation}
where $h_{tol}=500$ [m], $v_{tol}=16$ [m/s] and $\gamma_{tol}=0.5$ [deg].  The training episode is ended at the end of action-step $h$ if: (i) any of the path constraints from Eq. (\ref{eq:shuttle state path constraints}) or Eq. (\ref{eq:shuttle_heating_constraint}) are violated at any time $t$ for $t_{h-1} < t \leq t_h ~~,~~ t_{h} = \sum_{i=1}^{h}\tau_i$, or (ii) if each of the terminal tolerance criteria Eq. (\ref{eq:shuttle_terminal_tolerance}) are satisfied at the end of the action-step.  The Deep RL agent receives no further positive reward when the training episode ends, so the agent is therefore intrinsically incentivized to satisfy all path constraints while simultaneously maximize latitude, and ending the episode inside the ellipsoid $\Psi=1$ to maximize the terminal reward.  To promote robustness to initial condition variation, during training the initial conditions are perturbed uniformly as $x_0=x_0^* \pm \Delta x$ where $x_0^*= [h_0,v_0,\theta_0,\phi_0,\gamma_0,\psi_0, \alpha_0, \sigma_0]\Tr$ is the nominal initial condition Eq. (\ref{eq:shuttle_ics}) and $\Delta x$ is a uniformly distributed random variable with limits [4km, 390m/s, 2deg, 2deg, 2deg, 2deg, 0, 0]$\Tr$.

\subsubsection{Deep RL Comparison to NLP}\label{section:NumericalExp_DeepRLComparetoNLP}
We consider two Deep RL setups:
\begin{enumerate}
\item The \textit{HLAS-Dynamics} setup is the HLAS formulation described in Section \ref{section: HLAS Deep RL Formulation}, with the new ``anti-windup'' term Eq. (\ref{eq:AntiWindupLossTerm}) incorporated into the PPO surrogate objective Eq. (\ref{eq:PPOSurrogateObj}).  Each action sub- trajectory represents the desired dynamics function $z_h(t)=\dot y^{des}(t)$.  Referring to Section \ref{section: ProblemFormulation_actionspace}, the action space is $\mathcal{A}=(\tau, \{n^{z^{(1)}}_{i}\}_{i=1}^{p+1}, \{n^{z^{(2)}}_{i}\}_{i=1}^{p+1})$, where $z^{(1)}$ corresponds to $\dot \gamma^{des}$ and $z^{(2)}$ corresponds to $\dot \psi^{des}$.  Referring to Eq. (\ref{eq:z}) we choose to represent $z$ with a linear profile, $p=1$, corresponding to a second order polynomial representation of $y^{des}$ in Eq. (\ref{eq:y^des}).
\item The \textit{HLAS-Control} setup is similar to \textit{HLAS-Dynamics}, except each action sub-trajectory represents the control input $z_h(t)=u(t)$.  Referring to Section \ref{section: ProblemFormulation_actionspace}, the action space is $\mathcal{A}=(\tau, \{n^{z^{(1)}}_{i}\}_{i=1}^{p+1}, \{n^{z^{(2)}}_{i}\}_{i=1}^{p+1})$, where $z^{(1)}$ corresponds to $\alpha_{cmd}$ and $z^{(2)}$ corresponds to $\sigma_{cmd}$, with $p=1$.
\end{enumerate}

The Deep RL results are compared against a benchmark solution obtained through nonlinear programming (NLP) techniques.  To generate the NLP benchmark we utilize Pyomo [Hart et~al., 2017], a general, open-source, Python-based algebraic modeling language developed at Sandia National Labs.  The extension \verb|pyomo.dae| [Nicholson et~al., 2018] integrates with Pyomo to accommodate optimization problems with differential and algebraic equations (DAEs).  In this work a Backward Euler finite difference scheme with a fixed timestep of 2 seconds is used to approximate the shuttle DAEs.  The optimal time horizon length is determined by the solution of the NLP.  The final output from solving the NLP is a full trajectory along the horizon. The underlying optimization is solved through the interior-point solver \verb|IPOPT| [Biegler and Zavala, 2009].  Details of using \verb|pyomo.dae| in an aerospace optimal control context can be found in [Schlossman et~al., 2021].

The resulting trajectory from the NLP benchmark is shown in Fig. \ref{fig:shuttle_3d}. Also shown are the Deep RL-induced trajectories, which are produced by sequentially executing the trained policy and simulated environment (see Fig. \ref{fig:HLAS_diagram}).  The Deep RL trained policy executes in approximately 0.001 seconds for arbitrary inputs to the policy network.  Each trajectory is generated from the same fixed initial condition shown in Eq. (\ref{eq:shuttle_ics}).  The NLP solution achieves a final latitude of 32.19 [deg], and exactly satisfies the terminal condition.  Solution times for the NLP trajectory are approximately 30 seconds.  The \textit{HLAS-Control}-trained Deep RL policy achieves a final latitude of 31.05 [deg] (96\% of the NLP solution) with terminal error $h_{\text{err}}=128.8$ [m], $v_{\text{err}}=2.14$ [m/s], $\gamma_{\text{err}}=0.035$ [deg]. The \textit{HLAS-Dynamics}-trained Deep RL policy performs slightly worse, achieving a final latitude of 29.96 [deg].
\begin{figure}[h!]
\centering
\begin{minipage}{0.7\textwidth}
\includegraphics [trim = 5mm 5mm 5mm 10mm, clip, width=.99\textwidth]{shuttle_3d.pdf}
\end{minipage}
\captionsetup{format=hang}
\captionsetup{width=0.8\textwidth}
\caption{Trajectories generated from two Deep RL trained policies and the NLP solution.  Deep RL solutions are in blue (HLAS-Dynamics) and red (HLAS-Control), NLP solution in dashed grey.}
\label{fig:shuttle_3d}
\end{figure}
The vertical flight path and heading angles are shown in Fig. \ref{fig:shuttle_fpas}, along with the heating profile and Deep RL action duration distribution along the horizon.  The policy tends to prefer longer action-durations throughout the horizon, except where the policy is finely adjusting the trajectory (for example, in order to maximize terminal accuracy near horizon end).
\begin{figure}[h!]
\centering
\begin{minipage}{.65\textwidth}
\includegraphics [trim = 0mm 0mm 5mm 0mm, clip, width=.99\textwidth]{shuttle_fpas.pdf}
\end{minipage} \\
\begin{minipage}{.49\textwidth}
\includegraphics [trim = 0mm 0mm 5mm 0mm, clip, width=.99\textwidth]{shuttle_heating.pdf}
\end{minipage} ~
\begin{minipage}{.49\textwidth}
\includegraphics [trim = 0mm 0mm 5mm 0mm, clip, width=.99\textwidth]{shuttle_tau.pdf}
\end{minipage}
\captionsetup{format=hang}
\captionsetup{width=.85\textwidth}
\caption{Top: flight path angle, $\gamma$, and heading angle, $\psi$.  Left: Heating profiles, heating must not exceed a value of 80 BTU/ft$^2$-s along the horizon. Right: Action duration from Deep RL policy.  Deep RL solutions in blue (HLAS-Dynamics) and red (HLAS-Control), NLP solution in dashed grey.}
\label{fig:shuttle_fpas}
\end{figure}

\subsubsection{Deep RL Variants}\label{section:shuttle_ablation_studies}
We now compare the performance of several Deep RL variants against the \textit{HLAS-Dynamics} and \textit{HLAS-Control} methods.  The same network architecture (described in Section \ref{section:RLTrainingSetup}) was used throughout,  and algorithm hyperparamters were optimized when necessary.  The variants are described below:
\begin{enumerate}
\item The \textit{HLAS-Dynamics-Fixed $\tau$} variant is identical to \textit{HLAS-Dynamics} except the action duration is fixed as $\tau_h=4$ seconds.
\item The \textit{HLAS-Dynamics-NoAntiWindUp} variant is idential to \textit{HLAS-Dynamics} except the ``anti-windup'' loss term Eq. (\ref{eq:AntiWindupLossTerm}) is removed from the PPO surrogate objective Eq. (\ref{eq:PPOSurrogateObj}).
\item The \textit{Baseline} is a straightforward Deep RL setup where the attitude (angle of attack and bank angle) is commanded every timestep (two seconds).  The new ``anti-windup'' loss term Eq. (\ref{eq:AntiWindupLossTerm}) is not included in the PPO surrogate objective Eq. (\ref{eq:PPOSurrogateObj}).
\end{enumerate}
The training curves for each variant are shown in Fig. \ref{fig:shuttle_baselines_training}.  During training a 100-episode moving average of the undiscounted episode return is maintained, and the network parameters are saved during training whenever this average return reaches a new best value.
\begin{figure}[h!]
\centering
\begin{minipage}{0.7\textwidth}
\includegraphics [trim = 20mm 0mm 20mm 0mm, clip,width=.99\textwidth]{shuttle_baselines_training.pdf}
\end{minipage}
\captionsetup{format=hang}
\captionsetup{width=0.8\textwidth}
\caption{Training curves for Deep RL variants.  Each training is run until a reasonable steady state is reached, up to a maximum of roughly 35 hours.  Each data point (lighter curves) is the average return of the previous 100 training episodes.  The lighter curves are passed through a smoothing first order filter so the trending behavior can be easily observed.}
\label{fig:shuttle_baselines_training}
\end{figure}
Hyperparameters used for training each variant are shown in Table \ref{table:shuttle_baselines_hyperparameters}.  The entropy-based exploration bonus was included in the \textit{Baseline} and \textit{HLAS-Dynamics-Fixed-$\tau$} variants to encourage exploration.  Various combinations of hyperparameters were used for the Baseline training, and the combination producing the best training event is shown.  Although not shown, the HLAS formulations using the PPO training algorithm was found to perform consistently well under a reasonable range of hyperparameters.
\begin{table}[h!]
\renewcommand{\arraystretch}{1.2}
\caption{PPO hyperparameters. $\gamma=0.9999$, LR $=5\times 10^{-5}$ and PPO Clip $= 0.2$ in all methods.} \label{table:shuttle_baselines_hyperparameters}
\vspace{-4mm}
\centering % used for centering table
\begin{tabular}{lcccccc}
\hline {Method} & $\tau_{max}, \tau_{min}$ [s] & {VF Coef} & {Ent Coef} & AntiWindup & {Steps / Proc.} & Batch Size \\
\hline\vspace{-4mm} &&&&&& \\
HLAS-Dynamics & 30, 2 & 0.5 & 0 & Yes & 4096 & 128\\
HLAS-Control & 30, 2 & 0.5 & 0 & Yes & 4096 & 128\\
HLAS-Dyn-NoAntiWindUp & 30, 2 & 0.5 & 0 & No & 4096 & 128\\
HLAS-Dyn-Fixed-$\tau$ & 4 [fixed] & 0.5 & 0.001 & Yes & 4096 & 128\\
Baseline & 2 [fixed] & 100 & 0.001 & No & 8192 & 256\\
\hline
\end{tabular}
\end{table}

All five Deep RL methods were evaluated over 1000 randomized initial conditions  $x_0=x_0^* \pm 0.5\times\Delta x$, with $\Delta x$ as described above.  Post-training metrics are listed in Table \ref{table:shuttle_baselines_metrics}, which includes performance from the nominal initial condition, average episode return (given by (\ref{eq:discounted return}), evaluated with the discount factor $\gamma=1$), and the number of terminal misses where the state lands outside the bounds defined in (\ref{eq:shuttle_terminal_tolerance}).  Interestingly, \textit{HLAS-Control} slightly outperforms \textit{HLAS-Dynamics}, suggesting the policy network was sufficiently capable of representing the control dependencies within the system dynamics in this problem.  The \textit{HLAS-Dynamics-Fixed-$\tau$} variant performs poorly, a possible explanation is the shorter action duration will see fewer states due to the shorter flight time between action samples.  The \textit{HLAS-Dynamics-NoAntiWindUp} variant also performed poorly. Although it is not shown here, with \textit{HLAS-Dynamics-NoAntiWindUp} the action mean drifts beyond the action limit in several parts of the trajectory during training and never falls back within limits.

\begin{table}[h!]
\renewcommand{\arraystretch}{1.2}
\caption{Post training metrics.  Averages were taken over 1000 initial condition perturbations.} \label{table:shuttle_baselines_metrics}
\vspace{-4mm}
\centering % used for centering table
\begin{tabular}{lp{25mm}p{27mm}p{25mm}p{25mm}}
\hline {Method} & Latitude Achieved\newline $\Phi$ (nominal) & Terminal Miss Dist. \newline $\Psi^{-1}$ (nominal) &  Episode Return \newline $R_0$ (average) & \# Terminal Misses \newline (out of 1000) \\
\hline\vspace{-4mm} & & & & \\
HLAS-Dynamics & 29.96 deg & 0.305 & 34.67 & 0 \\
HLAS-Control & 31.05 deg & 0.679 & 35.26 & 0 \\
HLAS-Fixed-$\tau$ & 27.96 deg & 27.78 & 28.94 & 955 \\
HLAS-NoAntiWindUp & 24.72 deg & 2.08 & 28.65 & 7\\
Baseline & 26.39 deg & 1.67 & 29.38 & 21\\
\hline
\end{tabular}
\end{table}

In this section we have shown that Deep RL using the HLAS method can achieve similar optimality to NLP for a simple trajectory generation problem.  In contrast to the NLP solution which takes approximately 30 seconds to obtain, the trained Deep RL policy can be evaluated quickly for a given input (approximately 0.001 seconds) making it a strong candidate for real-time application.  In the next section we study a problem for which NLP struggles.

\subsection{Problem 2: Trajectory Planning for Debris Avoidance}\label{section:DebrisAvoidance}
In the previous section we showed our Deep RL approach can solve a difficult trajectory generation problem which is well-suited for traditional optimal control solvers (e.g. nonlinear programming).  We now demonstrate the effectiveness of our approach as a trajectory planner on a quasi-static map, and show that traditional optimal control solvers struggle here.  In this problem we employ the \textit{HLAS-Control} variant and again consider the shuttle vehicle described in Section \ref{section:NumericalExp_shuttle_descr}.  The objective is to guide the shuttle path to reach a specific reentry location $h_f=24384\pm1000$[m], $\theta_f=0\pm 1$ [deg], $\phi_f=50\pm 1$ [deg] while avoiding floating space debris represented as obstacles at fixed locations.  In Fig. \ref{fig:shuttle_obstacle_field} below the terminal location is shown as the blue dot and the debris field is shown as the collection of light red ellipses.
\begin{figure}[h!]
\centering
\begin{minipage}{0.7\textwidth}
\includegraphics [trim = 15mm 0mm 15mm 0mm, clip,width=.99\textwidth]{shuttle_obstacle_field.pdf}
\end{minipage}
\captionsetup{format=hang}
\captionsetup{width=.85\textwidth}
\caption{Floating debris field (light red ellipses) and reentry target location (blue dot).}
\label{fig:shuttle_obstacle_field}
\end{figure}
The initial shuttle location is randomly set along a semi-circle in the longitude-latitude plane of radius = $50 \pm 5$ [deg] surrounding the terminal location.  To ensure the problem is feasible the initial heading, $\psi_0$, points the shuttle vehicle towards the terminal location with a randomly set $\pm 5$ degree error.  Other initial conditions are $h_0=79248 \pm 2000$ [m], $v_0=7802 \pm 100$ [m/s], $\gamma=-1\pm 1$ [deg].

For Deep RL training we use a reward structure similar to that described in Section \ref{section:TraversedLatitude} which again follows the reward formulation Eq. (\ref{eq:HLAS Deep RL reward}): $L=0$ (running performance reward), $\Phi = 0$ (horizon reward), $\Psi = \left(\frac{h-h_f}{\bar h}\right)^2
+ \left(\frac{\theta-\theta_f}{\bar \theta}\right)^2
+ \left(\frac{\phi-\phi_f}{\bar \phi}\right)^2$ (terminal ellipsoid), $\hat\Psi^{-1} = \min\left(1, \Psi^{-1}\right)$ (terminal reward), $c_0 = 5$ (terminal reward weight). Here $\bar h=1000$ [m], $\bar \theta=1$ [deg] and $\bar\phi=1$ [deg] are scale factors.  Each training episode (successfully) ends when three simultaneous terminal conditions are satisfied: $|h-h_f|\leq 500$ [m], $|\theta-\theta_f|\leq 1$ [deg], $|\phi-\phi_f| \leq 1$ [deg].  The training episode ends prematurely if the agent intersects any of the obstacles shown in Fig. \ref{fig:shuttle_obstacle_field}.  Over approximately 16 hours of off-line training the agent learns to avoid the obstacles in order to maximize the terminal reward.  We train with the \textit{HLAS-Control} variant described in the previous section.  Referring to Section \ref{section: ProblemFormulation_actionspace}, the action space is $\mathcal{A}=(\tau, \{n^{z^{(1)}}_{i}\}_{i=1}^{p+1}, \{n^{z^{(2)}}_{i}\}_{i=1}^{p+1})$, where $z^{(1)}$ corresponds to $\alpha_{cmd}$ and $z^{(2)}$ corresponds to $\sigma_{cmd}$.  We choose $p=1$ corresponding to linear profiles for $\alpha_{cmd}$ and $\sigma_{cmd}$.

\subsubsection{Trajectory Planning Results}
Planning with a trained policy is straight-forward, which is a significant benefit of using Deep RL: trajectories are produced by sequentially executing the trained policy and then evaluating in the simulated environment, shown in Algorithm 1 (also Fig. \ref{fig:HLAS_diagram}).  An action is computed by evaluating the policy network's action-mean, $a=\mu_\theta(s)$, which takes approximately 0.001 seconds using the laptop hardware specified in Section \ref{section:NumericalExp}.  Profiles for angle of attack and bank angle are extracted from the action and passed to a simulator\footnote{In real applications the simulator would be the actual vehicle.}.  The simulator numerically propagates the shuttle vehicle dynamics Eq. (\ref{eq:sphericalDynamics}) along a sub-horizon of length $\tau$ determined by the action.  At the end of sub-horizon $\tau$, the updated state $s$ is passed to the policy network's action-mean and a new action is computed.  The planning process ends when the shuttle reaches its goal state or intersects a debris obstacle.  In this work the debris field and goal state are quasi-static, so these pieces of information are not provided to the policy (the policy network has essentially memorized this particular map configuration).  Given that end-to-end training completes in 16 hours on our laptop hardware, this is not an unreasonable assumption depending on the application.  However, in future work we seek to perform planning in response to dynamic obstacle fields and goal states, in which obstacles and goal states are provided as inputs to the policy, so that the policy can learn to generalize to arbitrary maps.

\begin{algorithm}[] %[H]
\setstretch{1.0}
\SetAlgoLined
\KwIn{Trained policy $\mu_\theta$, initial state $s$, simulator timestep $dt$}
\While{done is False}{
$a = \mu_\theta(s)$  \\
$\tau=a(1)$ \\
$p_{\alpha_{cmd}}=$ ProfileExtract($a(2)$) \\
$p_{\sigma_{cmd}}=$ ProfileExtract($a(3)$)\\
initialize: $t=0$\\
\While{$t\leq \tau$}{
$\alpha_{cmd}=p_{\alpha_{cmd}}(t)$, $\sigma_{cmd}=p_{\sigma_{cmd}}(t)$ \\
$s\leftarrow$ Simulator($\alpha_{cmd}, \sigma_{cmd}, s$) \\
$t\leftarrow t+dt$\\
\eIf{state limit violation \Or hit debris obstacle \Or reach goal}{
done = True\\
done = False
\caption{Trajectory planning with Deep RL policy}
\end{algorithm}

Trajectory planning results for the Deep RL planner for 100 initial conditions distributed over the semi-circle radius are shown in Fig. \ref{fig:shuttle_obstacle_field_deeprl_results}.  The initial conditions have been perturbed about nominal values as discussed above.
\begin{figure}[h!]
\centering
\begin{minipage}{0.8\textwidth}
\includegraphics [trim = 15mm 0mm 15mm 0mm, clip,width=.99\textwidth]{shuttle_obstacle_field_deeprl_results.pdf}
\end{minipage}
\captionsetup{format=hang}
\captionsetup{width=.85\textwidth}
\caption{Deep RL planner results from 100 randomly chosen initial conditions.}
\label{fig:shuttle_obstacle_field_deeprl_results}
\end{figure}
In each case the Deep RL planner is able to guide the shuttle to its reentry location within tolerance without intersecting any of the obstacles.  This problem was also attempted with the NLP methods discussed in Section \ref{section:NumericalExp_DeepRLComparetoNLP}.  Our implementation of NLP had difficulty solving this problem consistently, and required careful trial-and-error calibration of the numerical discretization scheme.  If the number of discretization points was too large or too small the solver would not converge.  Additionally, a set of discretization points which worked for one initial condition did not necessarily work for another.  This can be problematic for real-time planning as each NLP solution can take several minutes or longer if the solver gets ``stuck" while trying to find a solution.  In an attempt to promote convergence robustness with the NLP solver we also attempted a ``warm-start" procedure where the problem was first solved without obstacles present, and then re-solved with obstacles present using the first solution as an initial guess for the solver.  This often worked, but at times the trajectory became ``trapped" between two obstacles as shown in Fig. \ref{fig:shuttle_obstacle_field_results_trappedsoln}.

\begin{figure}[h!]
\centering
\begin{minipage}{0.7\textwidth}
\includegraphics [trim = 15mm 0mm 15mm 0mm, clip,width=.99\textwidth]{shuttle_obstacle_field_results_trappedsoln.pdf}
\end{minipage}
\captionsetup{format=hang}
\captionsetup{width=.9\textwidth}
\caption{Deep RL (red curve) comparison to NLP (grey curves) on a challenging initial condition.  The first solution generated with NLP does not consider the debris field (dashed grey).  This solution is provided as an initial guess to the NLP solver in a second iteration where the debris field is present (solid grey).  On the second iteration, the NLP trajectory becomes ``trapped" between two obstacles and cannot satisfy the boundary conditions.}
\label{fig:shuttle_obstacle_field_results_trappedsoln}
\end{figure}

\section{Conclusions}
This paper presents an approach for trajectory planning based on Deep RL, where actions are sub-trajectories of variable duration and shape.  The developed approach is referred to as the high level action space (HLAS) approach.  The HLAS approach appears to promote exploration within the environment as shown through a brief analysis and demonstrated empirically.  The HLAS approach is shown to significantly improve the Deep RL training process on a long-range trajectory generation problem, and also outperforms NLP on a trajectory planning problem involving obstacle avoidance.

We first assessed our HLAS approach on a well-studied long-range trajectory generation problem using the shuttle reentry vehicle.  First, off-line training is performed over a range of initial conditions to produce a state-feedback policy (policy training).  Second, after training is complete, the policy is rapidly evaluated in response to state inputs to induce a trajectory in a sequential manner (trajectory generation).  During the trajectory generation phase of Deep RL, the trained policy achieves optimality similar to NLP solutions, achieving 96\% of the NLP solution objective.  At trajectory generation time, the trained Deep RL policy can be evaluated rapidly making it a good candidate for real-time implementation.  The Deep RL policy can be evaluated in approximately 0.001 seconds, whereas the full NLP solution takes approximately 30 seconds to obtain on the same hardware.

On the same problem our approach outperforms a straight-forward implementation of Deep RL where actions are low level inputs, and action decisions are made every time step (in contrast to the HLAS approach where a constrained input sequence executes over multiple time steps).  The trained policy from HLAS Deep RL produces an 18\% improvement in the objective function, more efficiently and with less sensitivity to hyperparameters, compared to the straight-forward implementation of Deep RL.  We also develop a method to prevent action wind-up in the learned policy, which can be implemented directly with existing policy gradient algorithms such as PPO. We demonstrate the action anti-windup method produces a higher quality policy than when this method is absent.

Next, we investigated shuttle trajectory planning around a fixed obstacle field, with a fixed goal state.  Training is performed off-line using the HLAS Deep RL approach, and the trained policy is used directly as the planner from randomized initial conditions to the fixed goal state.  Our planning approach is shown to produce solutions more robustly and more rapidly than NLP.  The HLAS Deep RL policy can be evaluated in approximately 0.001 seconds on our computational hardware, whereas the NLP solution can take minutes to obtain (if the solution converges).  Our proposed method improves planning speed and stability by moving the computational burden off-line, where difficult scenarios are worked out through randomized exploration of the state-space.  Total training time is approximately 16 hours on our hardware for this problem.  The HLAS Deep RL trained policy is capable of memorizing the configuration of the environment over the course of training, so that planning can be performed from randomized initial conditions to the goal state.  In future work we seek to augment the HLAS Deep RL training, allowing the policy to plan in response to arbitrarily placed obstacles and goal states.

\section{Funding Statement}
This work was supported by the Laboratory Directed Research and Development program at Sandia National Laboratories, a multi-mission laboratory managed and operated by the National Technology and Engineering Solutions of Sandia LLC, a wholly owned subsidiary of Honeywell International Inc. for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

\begin{thebibliography}{74}
\newcommand{\enquote}[1]{``#1''}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\providecommand{\urlprefix}{URL }
\expandafter\ifx\csname urlstyle\endcsname\relax
\providecommand{\doi}[1]{\discretionary{}{}{}https://doi.org/#1}\else
\providecommand{\doi}[1]{\discretionary{}{}{}\urlstyle{rm}\url{https://doi.org/#1}}\fi

[Schlossman et~al., 2021]
Schlossman, R., Williams, K., Kozlowski, D., and Parish, J.,
\enquote{Open-Source, Object-Oriented, Multi-Phase Pseudospectral
Optimization Using Pyomo,} \emph{AIAA Scitech 2019 Forum}, 2021, p. 0430.

[Patterson and Rao, 2014]
Patterson, M.~A., and Rao, A.~V., \enquote{GPOPS-II: A MATLAB software for
solving multiple-phase optimal control problems using hp-adaptive Gaussian
quadrature collocation methods and sparse nonlinear programming,} \emph{ACM
Transactions on Mathematical Software (TOMS)}, Vol.~41, No.~1, 2014, pp.

[LaValle and Kuffner~Jr, 2001]
LaValle, S.~M., and Kuffner~Jr, J.~J., \enquote{Randomized kinodynamic
planning,} \emph{The international journal of robotics research}, Vol.~20,
No.~5, 2001, pp. 378--400.

[Karaman and Frazzoli, 2011]
Karaman, S., and Frazzoli, E., \enquote{Sampling-based algorithms for optimal
motion planning,} \emph{The international journal of robotics research},
Vol.~30, No.~7, 2011, pp. 846--894.

[Kuwata et~al., 2009]
Kuwata, Y., Teo, J., Fiore, G., Karaman, S., Frazzoli, E., and How, J.~P.,
\enquote{Real-time motion planning with applications to autonomous urban
driving,} \emph{IEEE Transactions on control systems technology}, Vol.~17,
No.~5, 2009, pp. 1105--1118.

[Hart et~al., 1968]
Hart, P.~E., Nilsson, N.~J., and Raphael, B., \enquote{A formal basis for the
heuristic determination of minimum cost paths,} \emph{IEEE transactions on
Systems Science and Cybernetics}, Vol.~4, No.~2, 1968, pp. 100--107.

[Dolgov et~al., 2008]
Dolgov, D., Thrun, S., Montemerlo, M., and Diebel, J., \enquote{Practical
search techniques in path planning for autonomous driving,} \emph{Ann Arbor},
Vol. 1001, No. 48105, 2008, pp. 18--80.

[Dijkstra et~al., 1959]
Dijkstra, E.~W., et~al., \enquote{A note on two problems in connexion with
graphs,} \emph{Numerische mathematik}, Vol.~1, No.~1, 1959, pp. 269--271.

[Frazzoli et~al., 2000]
Frazzoli, E., Dahleh, M.~A., and Feron, E., \enquote{Robust hybrid control for
autonomous vehicle motion planning,} \emph{Proceedings of the 39th IEEE
Conference on Decision and Control (Cat. No. 00CH37187)}, Vol.~1, IEEE, 2000,
pp. 821--826.

[Boyan and Moore, 1995]
Boyan, J., and Moore, A.~W., \enquote{Generalization in reinforcement learning:
Safely approximating the value function,} \emph{Advances in neural
information processing systems}, 1995, pp. 369--376.

[Bertsekas and Tsitsiklis, 1996]
Bertsekas, D.~P., and Tsitsiklis, J.~N., \emph{Neuro-Dynamic Programming},
Athena Scientific, 1996.

[Sutton et~al., 1998]
Sutton, R.~S., Barto, A.~G., et~al., \emph{Introduction to reinforcement
learning}, Vol. 135, MIT press Cambridge, 1998.

[Mnih et~al., 2015]
Mnih, V., Kavukcuoglu, K., Silver, D., Rusu, A.~A., Veness, J., Bellemare,
M.~G., Graves, A., Riedmiller, M., Fidjeland, A.~K., Ostrovski, G., et~al.,
\enquote{Human-level control through deep reinforcement learning,}
\emph{nature}, Vol. 518, No. 7540, 2015, pp. 529--533.

[Faust et~al., 2018]
Faust, A., Oslund, K., Ramirez, O., Francis, A., Tapia, L., Fiser, M., and
Davidson, J., \enquote{PRM-RL: Long-range robotic navigation tasks by
combining reinforcement learning and sampling-based planning,} \emph{2018
IEEE International Conference on Robotics and Automation (ICRA)}, IEEE, 2018,
pp. 5113--5120.

[Everett et~al., 2018]
Everett, M., Chen, Y.~F., and How, J.~P., \enquote{Motion planning among
dynamic, decision-making agents with deep reinforcement learning,} \emph{2018
IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)},
IEEE, 2018, pp. 3052--3059.

[Aradi, 2020]
Aradi, S., \enquote{Survey of deep reinforcement learning for motion planning
of autonomous vehicles,} \emph{IEEE Transactions on Intelligent
Transportation Systems}, 2020.

[G{\'o}mez et~al., 2012]
G{\'o}mez, M., Gonz{\'a}lez, R., Mart{\'\i}nez-Mar{\'\i}n, T., Meziat, D., and
S{\'a}nchez, S., \enquote{Optimal motion planning by reinforcement learning
in autonomous mobile vehicles,} \emph{Robotica}, Vol.~30, No.~2, 2012, p.

[Semnani et~al., 2020]
Semnani, S.~H., Liu, H., Everett, M., de~Ruiter, A., and How, J.~P.,
\enquote{Multi-agent motion planning for dense and dynamic environments via
deep reinforcement learning,} \emph{IEEE Robotics and Automation Letters},
Vol.~5, No.~2, 2020, pp. 3221--3226.

[Qureshi} et~al., 2021]
{Qureshi}, A.~H., {Miao}, Y., {Simeonov}, A., and {Yip}, M.~C., \enquote{Motion
Planning Networks: Bridging the Gap Between Learning-Based and Classical
Motion Planners,} \emph{IEEE Transactions on Robotics}, Vol.~37, No.~1, 2021,
pp. 48--66.
\newblock \doi{10.1109/TRO.2020.3006716}.

[Bellemare et~al., 2020]
Bellemare, M.~G., Candido, S., Castro, P.~S., Gong, J., Machado, M.~C., Moitra,
S., Ponda, S.~S., and Wang, Z., \enquote{Autonomous navigation of
stratospheric balloons using reinforcement learning,} \emph{Nature}, Vol.
588, No. 7836, 2020, pp. 77--82.

[LaFarge et~al., 2020]
LaFarge, N.~B., Miller, D., Howell, K.~C., and Linares, R., \enquote{Guidance
for closed-loop transfers using reinforcement learning with application to
libration point orbits,} \emph{AIAA Scitech 2020 Forum}, 2020, p. 0458.

[Gaudet et~al., 2020]
Gaudet, B., Linares, R., and Furfaro, R., \enquote{Deep reinforcement learning
for six degree-of-freedom planetary landing,} \emph{Advances in Space
Research}, Vol.~65, No.~7, 2020, pp. 1723--1741.

[Hovell and Ulrich, 2020]
Hovell, K., and Ulrich, S., \enquote{On Deep Reinforcement Learning for
Spacecraft Guidance,} \emph{AIAA Scitech 2020 Forum}, 2020, p. 1600.

[Sun and van Kampen, 2021]
Sun, B., and van Kampen, E.-J., \enquote{Reinforcement-Learning-Based Adaptive
Optimal Flight Control with Output Feedback and Input Constraints,}
\emph{Journal of Guidance, Control, and Dynamics}, 2021, pp. 1--8.

[Berner et~al., 2019]
Berner, C., Brockman, G., Chan, B., Cheung, V., D{\k{e}}biak, P., Dennison, C.,
Farhi, D., Fischer, Q., Hashme, S., Hesse, C., et~al., \enquote{Dota 2 with
large scale deep reinforcement learning,} \emph{arXiv preprint
arXiv:1912.06680}, 2019.

[Chiang} et~al., 2019]
{Chiang}, H.~L., {Hsu}, J., {Fiser}, M., {Tapia}, L., and {Faust}, A.,
\enquote{RL-RRT: Kinodynamic Motion Planning via Learning Reachability
Estimators From RL Policies,} \emph{IEEE Robotics and Automation Letters},
Vol.~4, No.~4, 2019, pp. 4298--4305.
\newblock \doi{10.1109/LRA.2019.2931199}.

[Goddard et~al., 2021]
Goddard, Z.~C., Wardlaw, K., Krishnan, R., Tsiotras, P., Smith, M.~R., Sena,
M.~R., Parish, J.~J., and Mazumdar, A., \enquote{Utilizing Reinforcement
Learning to Continuously Improve a Primitive-Based Motion Planner,}
\emph{AIAA Scitech 2021 Forum}, 2021, p. 1752.

[Li et~al., 2020]
Li, T., Lambert, N., Calandra, R., Meier, F., and Rai, A., \enquote{Learning
generalizable locomotion skills with hierarchical reinforcement learning,}
\emph{2020 IEEE International Conference on Robotics and Automation (ICRA)},
IEEE, 2020, pp. 413--419.

[Camci and Kayacan, 2019]
Camci, E., and Kayacan, E., \enquote{Learning motion primitives for planning
swift maneuvers of quadrotor,} \emph{Autonomous Robots}, Vol.~43, No.~7,
2019, pp. 1733--1745.

[Frazzoli, 2001]
Frazzoli, E., \enquote{Robust hybrid control for autonomous vehicle motion
planning,} Ph.D. thesis, Massachusetts Institute of Technology, 2001.

[Frazzoli et~al., 2002]
Frazzoli, E., Dahleh, M.~A., and Feron, E., \enquote{Real-time motion planning
for agile autonomous vehicles,} \emph{Journal of guidance, control, and
dynamics}, Vol.~25, No.~1, 2002, pp. 116--129.

[Gray et~al., 2012]
Gray, A., Gao, Y., Lin, T., Hedrick, J.~K., Tseng, H.~E., and Borrelli, F.,
\enquote{Predictive control for agile semi-autonomous ground vehicles using
motion primitives,} \emph{2012 American Control Conference (ACC)}, IEEE,
2012, pp. 4239--4244.

[Stulp and Schaal, 2011]
Stulp, F., and Schaal, S., \enquote{Hierarchical reinforcement learning with
movement primitives,} \emph{2011 11th IEEE-RAS International Conference on
Humanoid Robots}, IEEE, 2011, pp. 231--238.

[Ijspeert et~al., 2002]
Ijspeert, A.~J., Nakanishi, J., and Schaal, S., \enquote{Movement imitation
with nonlinear dynamical systems in humanoid robots,} \emph{Proceedings 2002
IEEE International Conference on Robotics and Automation (Cat. No.
02CH37292)}, Vol.~2, IEEE, 2002, pp. 1398--1403.

[Perk and Slotine, 2006]
Perk, B.~E., and Slotine, J.-J.~E., \enquote{Motion primitives for robotic
flight control,} \emph{arXiv preprint cs/0609140}, 2006.

[Grant and Antony, 2016]
Grant, M.~J., and Antony, T., \enquote{Rapid indirect trajectory optimization
of a hypothetical long range weapon system,} \emph{AIAA Atmospheric Flight
Mechanics Conference}, 2016, p. 0276.

[Kuwata et~al., 2005]
Kuwata, Y., Schouwenaars, T., Richards, A., and How, J., \enquote{Robust
constrained receding horizon control for trajectory planning,} \emph{AIAA
Guidance, Navigation, and Control Conference and Exhibit}, 2005, p. 6079.

[Schouwenaars et~al., 2004]
Schouwenaars, T., How, J., and Feron, E., \enquote{Receding horizon path
planning with implicit safety guarantees,} \emph{Proceedings of the 2004
American control conference}, Vol.~6, IEEE, 2004, pp. 5576--5581.

[Bryson and Ho, 1975]
Bryson, A.~E., and Ho, Y.-C., \emph{Applied Optimal Control}, Taylor
Scientific, 1975.

[Bertsekas et~al., 1995]
Bertsekas, D.~P., Bertsekas, D.~P., Bertsekas, D.~P., and Bertsekas, D.~P.,
\emph{Dynamic programming and optimal control}, Vol.~1, Athena scientific
Belmont, MA, 1995.

[Lewis et~al., 2012]
Lewis, F.~L., Vrabie, D., and Syrmos, V.~L., \emph{Optimal control}, John Wiley
\& Sons, 2012.

[Kelly, 2017]
Kelly, M., \enquote{An introduction to trajectory optimization: How to do your
own direct collocation,} \emph{SIAM Review}, Vol.~59, No.~4, 2017, pp.

[Grant et~al., 2011]
Grant, M., Clark, I., and Braun, R., \enquote{Rapid simultaneous hypersonic
aerodynamic and trajectory optimization using variational methods,}
\emph{AIAA Atmospheric Flight Mechanics Conference}, 2011, p. 6640.

[Sutton and Barto, 2018]
Sutton, R.~S., and Barto, A.~G., \emph{Reinforcement learning: An
introduction}, MIT press, 2018.

[Schulman et~al., 2015]
Schulman, J., Levine, S., Abbeel, P., Jordan, M., and Moritz, P.,
\enquote{Trust region policy optimization,} \emph{International conference on
machine learning}, PMLR, 2015{\natexlab{a}}, pp. 1889--1897.

[Littman, 1994]
Littman, M.~L., \enquote{Markov games as a framework for multi-agent
reinforcement learning,} \emph{Machine learning proceedings 1994}, Elsevier,
1994, pp. 157--163.

[Sutton et~al., 1999]
Sutton, R.~S., McAllester, D.~A., Singh, S.~P., Mansour, Y., et~al.,
\enquote{Policy gradient methods for reinforcement learning with function
approximation.} \emph{NIPs}, Vol.~99, Citeseer, 1999, pp. 1057--1063.

[Peters and Schaal, 2006]
Peters, J., and Schaal, S., \enquote{Policy gradient methods for robotics,}
\emph{2006 IEEE/RSJ International Conference on Intelligent Robots and
Systems}, IEEE, 2006, pp. 2219--2225.

[Silver et~al., 2014]
Silver, D., Lever, G., Heess, N., Degris, T., Wierstra, D., and Riedmiller, M.,
\enquote{Deterministic policy gradient algorithms,} \emph{International
conference on machine learning}, PMLR, 2014, pp. 387--395.

[Lillicrap et~al., 2015]
Lillicrap, T.~P., Hunt, J.~J., Pritzel, A., Heess, N., Erez, T., Tassa, Y.,
Silver, D., and Wierstra, D., \enquote{Continuous control with deep
reinforcement learning,} \emph{arXiv preprint arXiv:1509.02971}, 2015.

[Mnih et~al., 2016]
Mnih, V., Badia, A.~P., Mirza, M., Graves, A., Lillicrap, T., Harley, T.,
Silver, D., and Kavukcuoglu, K., \enquote{Asynchronous methods for deep
reinforcement learning,} \emph{International conference on machine learning},
PMLR, 2016, pp. 1928--1937.

[Duan et~al., 2016]
Duan, Y., Chen, X., Houthooft, R., Schulman, J., and Abbeel, P.,
\enquote{Benchmarking deep reinforcement learning for continuous control,}
\emph{International conference on machine learning}, PMLR, 2016, pp.

[Schulman et~al., 2015]
Schulman, J., Moritz, P., Levine, S., Jordan, M., and Abbeel, P.,
\enquote{High-dimensional continuous control using generalized advantage
estimation,} \emph{arXiv preprint arXiv:1506.02438}, 2015{\natexlab{b}}.

[Biegler and Zavala, 2009]
Biegler, L.~T., and Zavala, V.~M., \enquote{Large-scale nonlinear programming
using IPOPT: An integrating framework for enterprise-wide dynamic
optimization,} \emph{Computers \& Chemical Engineering}, Vol.~33, No.~3,
2009, pp. 575--582.

[Dever et~al., 2006]
Dever, C., Mettler, B., Feron, E., Popovic, J., and McConley, M.,
\enquote{Nonlinear trajectory generation for autonomous vehicles via
parametrized maneuver classes,} \emph{Journal of guidance, control, and
dynamics}, Vol.~29, No.~2, 2006, pp. 289--302.

[Paranjape et~al., 2015]
Paranjape, A.~A., Meier, K.~C., Shi, X., Chung, S.-J., and Hutchinson, S.,
\enquote{Motion primitives and 3D path planning for fast flight through a
forest,} \emph{The International Journal of Robotics Research}, Vol.~34,
No.~3, 2015, pp. 357--377.

[Haarnoja et~al., 2018]
Haarnoja, T., Zhou, A., Hartikainen, K., Tucker, G., Ha, S., Tan, J., Kumar,
V., Zhu, H., Gupta, A., Abbeel, P., et~al., \enquote{Soft actor-critic
algorithms and applications,} \emph{arXiv preprint arXiv:1812.05905}, 2018.

[Fujimoto et~al., 2018]
Fujimoto, S., Hoof, H., and Meger, D., \enquote{Addressing function
approximation error in actor-critic methods,} \emph{International Conference
on Machine Learning}, PMLR, 2018, pp. 1587--1596.

[Schulman et~al., 2017]
Schulman, J., Wolski, F., Dhariwal, P., Radford, A., and Klimov, O.,
\enquote{Proximal policy optimization algorithms,} \emph{arXiv preprint
arXiv:1707.06347}, 2017.

[Achiam, 2018]
Achiam, J., \enquote{{Spinning Up in Deep Reinforcement Learning},} 2018.

[Gangwani and Peng, 2018]
Gangwani, T., and Peng, J., \enquote{Policy optimization by genetic
distillation,} \emph{6th International Conference on Learning
Representations, ICLR 2018}, 2018.

[Wang et~al., 2019]
Wang, Y., He, H., Tan, X., and Gan, Y., \enquote{Trust region-guided proximal
policy optimization,} \emph{arXiv preprint arXiv:1901.10314}, 2019.

[Williams and Peng, 1991]
Williams, R.~J., and Peng, J., \enquote{Function optimization using
connectionist reinforcement learning algorithms,} \emph{Connection Science},
Vol.~3, No.~3, 1991, pp. 241--268.

[Raffin et~al., 2019]
Raffin, A., Hill, A., Ernestus, M., Gleave, A., Kanervisto, A., and Dormann,
N., \enquote{Stable Baselines3,}
\url{https://github.com/DLR-RM/stable-baselines3}, 2019.

[Chou et~al., 2017]
Chou, P.-W., Maturana, D., and Scherer, S., \enquote{Improving stochastic
policy gradients in continuous control with deep reinforcement learning using
the beta distribution,} \emph{International conference on machine learning},
PMLR, 2017, pp. 834--843.

[Fujita and Maeda, 2018]
Fujita, Y., and Maeda, S.-i., \enquote{Clipped action policy gradient,}
\emph{International Conference on Machine Learning}, PMLR, 2018, pp.

[\AA}str{\"o}m and Murray, 2021]
{\AA}str{\"o}m, K.~J., and Murray, R.~M., \emph{Feedback systems: an
introduction for scientists and engineers}, Princeton university press, 2021.

[Bertsekas, 2016]
Bertsekas, D.~P., \emph{Nonlinear Programming}, Athena Scientific, 2016.

[Chong and Zak, 2008]
Chong, E.~K., and Zak, S.~H., \emph{An introduction to optimization}, John
Wiley \& Sons, 2008.

[Betts, 2010]
Betts, J.~T., \emph{Practical methods for optimal control and estimation using
nonlinear programming}, SIAM, 2010.

[Tip{\`a}n et~al., 2020]
Tip{\`a}n, S., Theodoulis, S., Thai, S., and Proff, M., \enquote{Nonlinear
Dynamic Inversion Flight Control Design for Guided Projectiles,}
\emph{Journal of Guidance, Control, and Dynamics}, Vol.~43, No.~5, 2020, pp.

[Regan and Anandakrishnan, 1993]
Regan, F.~J., and Anandakrishnan, S.~M., \emph{Dynamics of Atmospheric
Re-Entry}, American Institute of Aeronautics and Astronautics, Inc., 1993.

[Hart et~al., 2017]
Hart, W.~E., Laird, C.~D., Watson, J.-P., Woodruff, D.~L., Hackebeil, G.~A.,
Nicholson, B.~L., and Siirola, J.~D., \emph{Pyomo-optimization modeling in
python}, Vol.~67, Springer, 2017.

[Nicholson et~al., 2018]
Nicholson, B., Siirola, J.~D., Watson, J.-P., Zavala, V.~M., and Biegler,
L.~T., \enquote{pyomo. dae: a modeling and automatic discretization framework
for optimization with differential and algebraic equations,}
\emph{Mathematical Programming Computation}, Vol.~10, No.~2, 2018, pp.

\end{thebibliography}

\end{document}