text:
\section{Appendix}

\input{background_related_work}

\subsection{Algorithm for Model-Based RL}
\label{sec:modelRL}
RL follows the formal framework of a Markov decision process (MDP) which is represented using a tuple $M := ({\mathbb{S}}, {\mathbb{A}}, f, R)$ where ${\mathbb{S}}$, ${\mathbb{A}}$ and $f$ denote the state space, action space and true system dynamics, as defined in Section 2.
The last element in the MDP tuple is the reward function $R: {\mathbb{S}} \times {\mathbb{A}} \rightarrow {\mathbb{R}}$ which quantifies the \textit{goodness} of taking an action ${\bm{a}}_t$ in state ${\bm{s}}_t$, by computing the scalar value $r_t := R({\bm{s}}_t, {\bm{a}}_t)$.
Then the goal of the RL agent is to develop a policy (controller) $\pi$ which maximizes the expectation of the total return $G = \sum_{t=1}^T \gamma^{t-1} r_t$ over time $T$, where $\gamma \in [0, \; 1)$ is the discount factor.
In the model-based RL framework, the policy $\pi$ can simply be a planning algorithm or it can be a function approximator.

Under the model-based RL framework, the true dynamics function $f$ is unknown and needs to be learnt from observed data samples collected from the actual environment.
The learnt model is represented using a function approximator as $\hat{f}_{\bm{\theta}}$.
A generic algorithm for model-based RL is provided in Algorithm 1.
To train the model, an initial training data set $D^{(0)}$ is first generated by sampling trajectories $\xi$ using a randomly initialized policy $\pi$.
In any iteration $k$, the data tuples $D^{(k)} = \{ ({\bm{s}}_t, {\bm{a}}_t, \Delta {\bm{s}}_{t+1})_i \}_{i=1}^{N_k}$ are used to train $\hat{f}_{\bm{\theta}}$ by minimizing the prediction error $E$
using a  supervised learning technique.
This trained model $\hat{f}_{\bm{\theta}}$ is then used to update the policy $\pi$ or used by a planning algorithm to generate new trajectories and improve the return $G$ obtained by the agent.
The data tuples obtained from the new trajectories are appended to the existing training data set to create $D^{(k+1)} = D^{(k)} \cup \{ \xi_i \}$ and the process continues until the performance of the controller converges.

\begin{algorithm}
\SetAlgoLined
\KwInput{$\pi^{(0)}$, $\hat{f}^{(0)}_{\bm{\theta}}$}
 Generate trajectories $\xi$ using base policy $\pi^{(0)}$\;
 Collect and construct initial data set $D^{(0)} = \{ ({\bm{s}}_t, {\bm{a}}_t, \Delta {\bm{s}}_{t+1})_i \}_{i=1}^{N_0}$\;
 k = 0\;
 \While{$\Delta E > \epsilon$}{
  Learn dynamics model $\hat{f}^{(k)}_{{\bm{\theta}}}({\bm{s}}, {\bm{a}})$ over $D^{(k)}$ by minimizing error $E(D^{(k)})$\;
  \While{$\Delta G > \delta$}{
  At time $t$ plan through $\hat{f}^{(k)}_{\bm{\theta}}({\bm{s}},{\bm{a}})$ to select action ${\bm{a}}_t$\;
  Execute action ${\bm{a}}_t$ and observe reward $r_t$ and next state ${\bm{s}}_{t+1}$ (MPC)\;
  Obtain updated policy $\pi^{(k+1)}$ using backpropagation to maximize return $G = \sum_{t=1}^T\gamma^{t-1} r_t$\;
  Create data set $D^{(k+1)}$ by appending new observed data points\;
  }
  k = k+1\;
 }
 \caption{Generic model-based RL algorithm}
 \label{algo1}
\end{algorithm}

\subsection{Implementation of Proposed Safety Verification Framework}
\label{sec:implementation}
To implement the proposed safety verification framework in Fig. \ref{fig:flowchart}, the MATLAB helperOC \cite{helperOC_tool} developed and maintained by the HJ Reachability Group on Github was used.
This toolbox is dependent on the level set toolbox developed by \citet{LSM_toolbox}.
In the helperOC toolbox, the \textit{level set method} is used to provide a numerical solution to the PDE in (\ref{eqn:BRT_final}).
Level set methods have been used to track and model the evolution of dynamic surfaces (here the reachable tube) by representing the dynamic surface implicitly using a level set function at any instant of time \citep{osher1988fronts, osher2004level, sethian1999level}.
There are several advantages to using the level set algorithms, including the fact that they can represent non-linear, non-convex surfaces even if the surface merges or divides over time.
Additionally, it is easy to extend the theory to the higher dimensional state space.
However, the computational time increases exponentially with a dimension of the state space.
The helperOC \citet{helperOC_tool} has demonstrated results on problems with 10 dimensional space and since most physical systems can be represented within this limit, it is sufficient for the problems of interest in this work.

The helperOC \citet{helperOC_tool} comes with inbuilt numerical schemes to compute the Hamiltonian $H$, the spacial derivatives $\nabla V_R$ and the time derivative $\frac{dV_R}{dt}$ in (\ref{eqn:BRT_final}).
This toolbox is modified to solve the problem in (\ref{eqn:BRT_final}) by substituting the true system dynamics $f$ with the augmented learnt dynamics $\hat{f}^{(r)}_{\bm{\theta}}({\bm{s}}, {\bm{a}}, {\bm{d}})$ and computing the action at any state ${\bm{s}}$ as ${\bm{a}} = \pi({\bm{s}})$.
Subsequently, to determine the level set function at any time instant, $\frac{dV}{dt}$ is computed using the Runge-Kutta integration scheme while the gradient of the level set function $\nabla V_R$ is computed using an upwind finite difference scheme.
Additionally, a Lax-Friedrichs approximation is used to stably evaluate the Hamiltoniaan \citep{Lax_fried}.

\subsection{Experiment Results: Safe Land Navigation}
\label{sec:results_land}
In the land navigation problem, a ground robot is moving in an indoor environment with obstacles.
The ground vehicle has the state vector ${\bm{s}} = [x, y]^T$ representing the $x$ and $y$ coordinates of the car.
The control vector ${\bm{a}} = [v, \psi]$ comprises of the velocity $v$ and the direction $\psi$ of the vehicle.
The problem setting in shown in Fig. \ref{fig:2D_init}.
The set of initial and goal states are represented by circles centered at coordinates $(0.0, 0.0)$ and $(6.0, 4.5)$, respectively, and with radii of $0.7$ and $0.5$ units, respectively.
There are two rectangular obstacles in the environment with centers at $(1.5, 4.5)$ and $(4.0, 1.5)$.

To demonstrate the working of the proposed safety verification framework, different models and controllers were trained for the 2D land navigation problem setting.
Three different models were trained over training data sets of different sizes.
Model 1 was trained over 300 data samples, model 2 over 600 data samples, and model 3 over 1000 data samples.
To estimate the modeling error for all the three models, a $\pm 3\sigma$ bound was computed by first computing the standard deviation (s.d.) over the validation data sets (Table \ref{tab:sd_2d}).
\begin{table}
\begin{center}
\caption{Standard deviations (s.d.) computed for different models.} \label{tab:sd_2d}
\begin{tabular}{| c | c | c | c | }
  \hline
  & Model 1 & Model 2 & Model 3 \\
  \hline
  s.d. for state $x$ & 0.00114 & 0.00045 & 0.00031 \\
 \hline
  s.d. for state $y$ & 0.00120 & 0.00053 & 0.00042 \\
 \hline
\end{tabular}
\end{center}
\end{table}
Then, for each model, three different policies were sampled after training them for 100, 300 and 600 episodes, respectively.
We first present the FRT results for all controllers correspondong to every model, followed by the BRT results.

\begin{figure}[h]
\begin{center}
\includegraphics[width=0.32\textwidth]{images/FRT_300_100_v2.png}
\includegraphics[width=0.32\textwidth]{images/FRT_300_300_v2.png}
\includegraphics[width=0.32\textwidth]{images/FRT_300_600_v2.png}
\\
\includegraphics[width=0.32\textwidth]{images/FRT_600_100_rr.png}
\includegraphics[width=0.32\textwidth]{images/FRT_600_300_rr.png}
\includegraphics[width=0.32\textwidth]{images/FRT_600_600_rr.png}
\\
\includegraphics[width=0.32\textwidth]{images/FRT_1000_100_rr.png}
\includegraphics[width=0.32\textwidth]{images/FRT_1000_300_rr.png}
\includegraphics[width=0.32\textwidth]{images/FRT_1000_600_rr.png}
\end{center}
\caption{The FRTs computed over the augmented learnt dynamics and the true system dynamics are presented. (Top Row) The FRTs corresponding to policy 1 (left), policy 2 (middle) and policy 3 (right) trained over learnt model 1. (Middle Row) The FRTs corresponding to policy 1 (left), policy 2 (middle) and policy 3 (right) trained over learnt model 2. (Bottom Row) The FRTs corresponding to policy 1 (left), policy 2 (middle) and policy 3 (right) trained over learnt model 3.}
\label{fig:FRT_2D}
\end{figure}

The FRT results for all the different models and controllers are shown in Fig. \ref{fig:FRT_2D}.
Each plot in Fig. \ref{fig:FRT_2D} shows two FRTs - the first corresponding to the augmented learnt dynamics and the second corresponding to the true system dynamics.
This is done to validate the result of the proposed safety verification framework with the ground truth classification of $\pi$ as safe or unsafe.
In general, it is observed that the FRT computed by the proposed framework, over the augmented learnt dynamics, closely resembles the FRT computed over the true system dynamics.
Thus, in every case, classification of $\pi$ determined by the proposed framework is consistent with the ground truth result.
All the above policies are found to be unsafe because the FRT intersects with obstacle 2.
In all the cases, the policies are safe with respect to obstacle 1.
Since all the above policies are unsafe, the BRT analysis is performed to determine ${\mathbb{S}}_{safe}$ for each policy.

The results corresponding to the determination of ${\mathbb{S}}_{safe}$, using the proposed BRT formulation, for each of the 3 controllers trained over learnt models 1, 2 and 3 are presented in Fig. \ref{fig:BRT_land_300}.
${\mathbb{S}}_{safe}$ determined by the proposed method is compared with the ground truth (G.T.) data, which was generated by using Monte-Carlo samples of possible trajectories from ${\mathbb{S}}_0$.
It is observed that the most accurate results are obtained for model 1.
For model 2, a small error due to modeling uncertainty is observed while determining ${\mathbb{S}}_{safe}$.
This error due to modeling uncertainty is related to how conservative the bound on the modeling error is.
For model 3, the proposed BRT formulation correctly identifies policy 1 as completely unsafe, since ${\mathbb{S}}_{safe} = \emptyset$.
For the remaining policies, a small error is observed while determining ${\mathbb{S}}_{safe}$.

It can be seen from Fig. \ref{fig:BRT_land_300} that there is no consistent trend for the fraction of safe initial states across, either the models trained with increasing number of sample points (top to bottom) or, across the policies trained for increasing number of iterations (left to right).
This highlights one of the primary concerns for the lack of safety in RL based controllers: even though the cost function is designed to mitigate collisions and head towards the target, the controller learnt using it may not be safe because of various reasons stemming from state-representations learned by NN based policies and model dynamics, the hyper-parameters being used, or the optimization getting stuck in local minima, etc.
In such scenarios, the proposed framework can be used to check before deployment whether, the developed RL controller abides by the required safety constraints or not.

\begin{figure}[h]
\begin{center}
\includegraphics[width=0.32\textwidth]{images/BRT_300_100_v2.png}
\includegraphics[width=0.32\textwidth]{images/BRT_300_300_v2.png}
\includegraphics[width=0.32\textwidth]{images/BRT_300_600_v2.png}
\\
\includegraphics[width=0.32\textwidth]{images/BRT_600_100_rr.png}
\includegraphics[width=0.32\textwidth]{images/BRT_600_300_rr.png}
\includegraphics[width=0.32\textwidth]{images/BRT_600_600_rr.png}
\\
\includegraphics[width=0.32\textwidth]{images/BRT_1000_100_rr.png}
\includegraphics[width=0.32\textwidth]{images/BRT_1000_300_rr.png}
\includegraphics[width=0.32\textwidth]{images/BRT_1000_600_rr.png}
\end{center}
\caption{ Comparison of ${\mathbb{S}}_{safe}$ computed by the proposed BRT formulation with the ground truth (G.T.) data.  (Top Row) The BRT computed for policy 1 (left), policy 2 (middle) and policy 3 (right) trained over learnt model 1. (Middle Row) The BRT computed for policy 1 (left), policy 2 (middle) and policy 3 (right) trained over learnt model 2. (Bottom Row) The BRT computed for policy 1 (left), policy 2 (middle) and policy 3 (right) trained over learnt model 3.}
\label{fig:BRT_land_300}
\end{figure}

\subsection{Experiment Results: Safe Aerial Navigation}
\label{sec:results_aerial}
In the aerial navigation problem, an unmanned aerial vehicle (UAV) is flying in an urban environment.
To ensure a safe flight, the UAV has to avoid collision with obstacles, like buildings and trees, in it's path.
This problem is abstracted in a continuous state and action domain, where the state vector ${\bm{s}} = [x, y, z]$ gives the 3D coordinates of the UAV and the action vector ${\bm{a}} = [v, \psi, \phi]$ comprises of the velocity, heading angle and pitch angle.
The UAV has to takeoff from the ground and reach a desired set of goal states.
Both the initial set of states and the goal states are represented using cuboids centered at coordinates $(0.0, 0.0, 0.0)$ and $(3.8, 4.5, 4.5)$, respectively.
There are two cylindrical obstacles in the environment, each of height $6$ units, which are centered at coordinates $(2.0, 4.0, 0.0)$ and $(4.0, 3.0, 0.0)$, respectively, with a radius of $0.5$ units.
This problem setting is presented in Fig. \ref{fig:3D_init}.

A learnt model represented by an NN is trained over 1000 data points and the policy $\pi$ trained over this model is evaluated.
The FRT computed by the proposed framework, over the augmented learnt dynamics, is presented from different perspectives in Fig. \ref{fig:FRT_3D}.
While the computed FRT doesn't intersect with obstacle 1, it clearly intersects with obstacle 2.
Therefore, $\pi$ is deemed unsafe.
The next step is to compute the BRT and determine ${\mathbb{S}}_{safe}$ for the unsafe policy $\pi$ as shown in Fig. \ref{fig:FRT_3D}.
It is observed that the proposed method under approximates ${\mathbb{S}}_{safe}$ when compared with the G.T. data.

\begin{figure}[h]
\begin{center}
\includegraphics[width=0.32\textwidth]{images/FRT_3D_pic2.png}
\includegraphics[width=0.32\textwidth]{images/FRT_3D_pic4.png}
\includegraphics[width=0.32\textwidth]{images/FRT_3D_pic6.png}
\\
\includegraphics[width=0.32\textwidth]{images/BRT_3D_pic2.png}
\includegraphics[width=0.32\textwidth]{images/BRT_3D_pic5.png}
\includegraphics[width=0.32\textwidth]{images/BRT_3D_pic4.png}
\end{center}
\caption{(Top Row) Views from different angles for the FRT computed by the proposed method over the augmented learnt system dynamics.
It can be noticed that FRT intersects with obstacle 2, therefore, the policy is deemed unsafe.
(Bottom Row) Views from different angles of the BRT computed by the proposed method over the augmented learnt system dynamics.
The proposed method under approximates ${\mathbb{S}}_{safe}$when compared with the G.T. data, marked in green. }
\label{fig:FRT_3D}
\end{figure}
\subsection{Extended Related Work} \label{sec:ERW}

NNs are one of the most commonly used function approximators for RL algorithms.
Therefore, often, verifying the safety of an RL controller reduces to verifying the safety of the underlying NN model or controller.
In this section, we briefly review the techniques used to evaluate the safety of NNs.

In recent years, NN-based supervised learning of image classifiers have found applications in autonomous driving and perception modules.
These are safety critical systems and it is necessary to verify the robustness of NN-based classifiers to minimize misclassification errors after they are deployed \citep{kwiatkowska2019safety}.
Some of the recent works have relied upon using existing optimization frameworks, such as \textit{linear programming} and \textit{satisfiability modulo theory} (SMT), to determine adversarial examples in the neighborhood of input images \citep{huang2017safety, bastani2016measuring, pulina2012challenging}.
However, these techniques usually
require an approximation of the constraints on the NN model.
\citet{katz2017reluplex} improved the scalability of SMT verification algorithms and \citet{sun2019formal} used \textit{satisfiability modulo convex} (SMC) over a partitioned workspace, to determine the set of safe initial states for a robot, but their algorithms were limited to NNs with the $\texttt{ReLU}$ activation function.
However, the problem of interest in this work is more general than the \textit{single-step} classification problems considered in these related work, and instead we look at RL problems which require \textit{multi-step} decision making in real-valued domains.

For real-valued output domains, the safety property needs to be verified for a set of states as opposed to a single label (i.e., a single point) in the output domain.
Therefore, to verify safety in the continuous problem domain, most works \citep{xiang2018reachable, tran2019parallelizable, xiang2018output, tran2019safety} first compute the forward reachable tube (FRT) of an NN, i.e, the set of output values achieved by the NN.
Then, the NN is said to be unsafe if the FRT intersects with an unsafe set.
\citet{xiang2018reachable}, \citet{tran2019parallelizable} and \citet{xiang2018output} follow the FRT algorithm by discretizing the input space into subspaces and executing forward pass over the NN to evaluate its output range.
However, the above works use convex polyhedrons to approximate non-convex reachable sets and the approximation error compounds over each layer of the NN.
\citet{tran2019safety} used a star-shaped polygon approach to generate the FRT which gives a less conservative solution than most polyhedra-based reachable tube computation.
Additionally, some of the algorithms are applicable specifically to $\texttt{ReLU}$ activation functions \citep{xiang2018reachable, tran2019parallelizable}.
In comparison, the algorithm proposed in this work computes the exact reachable tube and can handle most commonly used activation functions, including $\texttt{tanh, sigmoid, ReLu,}$ etc.

\begin{comment}
The problem of safety verification in RL controllers has gained attention in the last few years \citep{kwiatkowska2019safety}.
Since the increasing use of neural network (NN) controllers in complex applications like autonomous driving and perception modules, much of the initial work in safety verification of NNs was focused on image classification problems.
In contrast, in this work we consider a continuous output domain.
For classification problems, the most popular way to check the safety of a controller is to find an adversarial input for which the controller falsely classifies the image
and thereby, determine a safe neighborhood around input images for which the classification label is always correct \cite{huang2017safety, ruan2018reachability, wu2020game, wicker2018feature, wang2018efficient, bastani2016measuring}.
To this end, the authors of \cite{huang2017safety, pulina2012challenging, bastani2016measuring, katz2017reluplex, scheibler2015towards}
use existing solvers like \textit{satisfiability modulo theory} (SMT) to verify the properties of NNs, but their application is limited to NNs with small number of nodes.
The authors of \cite{wicker2018feature, wu2020game} formulated the problem as a two player game and analyzed the resilience and robustness of a deep neural network (DNN).
Some of these works also divide the input space into subspaces \cite{huang2017safety, wang2018efficient} to propagate the input vectors efficiently and while this method can work for classification problems, it can result in too conservative solutions for continuous systems.
However, most of these works are focused on obtaining an adversarial example rather evaluating the system safety.

Recent work for continuous domain focus on determining the output range of a NN by computing a forward reachable set/tube \cite{xiang2017reachable, xiang2018reachable, tran2019safety, tran2019parallelizable, xiang2018output}.
A common technique is to conduct a layer-by-layer analysis of NNs which ensures scalability \cite{xiang2017reachable, xiang2018reachable}, however, it also means that the set approximation error keeps growing with each layer.
Some algorithms are developed for specific activation functions. e.g., the analysis of ReLu activation function has been well studied in literature \cite{katz2017reluplex, xiang2017reachable, tran2019safety, sun2019formal, xiang2018reachable}.
Some work have also tried to abstract or linearize a NN \cite{pulina2010abstraction, wang2018efficient} to make the problem of output range determination tractable.
Additionally, using polygon approximation for computing forward reachable sets has become a popular technique \cite{tran2019safety, xiang2017reachable, xiang2018reachable, tran2019parallelizable} but can lead to conservative solutions.

As stated before, in this work we are focused on not only verifying the safety (constraint satisfaction) property of a model-based controller but to also determine the set of safe and unsafe initial states.
To this end, we propose a framework to construct a backward reachable tube (BRT) from an unsafe target set and from it's intersection with the initial set of states, determine the subset of unsafe states.
We use reachable set analysis because it is an exhaustive verification framework \cite{maler2008computing} which can be leveraged to analyze properties of non-linear, non-convex functions, such as NNs, too.
We use the Hamilton-Jacobi-Isaac (HJI) partial differential equation (PDE) formulation to solve for the BRT over a finite time interval, by using the level set methods to numerically approximate the reachable set \cite{mitchell2005time, akametalu2014reachability}.
This method of computing the reachable set gives a less conservative solution than most polygon approximation based approaches because the level set functions can represent non-convex set boundaries while handling uncertain parameters in the system.

\section{Notation and Problem Statement}

Let ${\mathbb{S}} \subset {\mathbb{R}}^n$ denote the state space and ${\mathbb{A}} \subset {\mathbb{R}}^m$ denote the action space.
Let ${\mathbb{S}}_0 \subset {\mathbb{S}}$ denote the set of initial states for the agent and let ${\mathbb{L}} \subset {\mathbb{R}}^n$ denote the target set which the agent is trying to reach.
Let the system be under some bounded disturbance ${\mathbb{D}} \subset {\mathbb{R}}^n$ such that the agent follows the deterministic system dynamics $f_r : {\mathbb{S}} \times {\mathbb{A}} \times {\mathbb{D}} \rightarrow {\mathbb{S}}$.
Then, the backward reachable tube (BRT) is defined as the set of all states ${\bm{s}} \in {\mathbb{S}}$ from which the agent can reach the given target set ${\mathbb{L}}$ while following the system dynamics $f_r({\bm{s}}, {\bm{a}}, {\bm{d}})$ under a finite time $T < \infty$, despite the disturbance ${\bm{d}} \in {\mathbb{D}}$.
\begin{equation}
    \begin{split}
    {\mathbb{R}}_B(T) := \{ {\bm{s}}_0: \exists {\bm{a}}, \; s.t. \; \forall {\bm{d}}, \; {\bm{s}}(\cdot) \; \text{ satisfies }\; \dot{{\bm{s}}}=f_r({\bm{s}},{\bm{a}}, {\bm{d}}), \\
    {\bm{s}}(0)={\bm{s}}_0; \; \exists t \in [-T, 0], \; s.t. \; {\bm{s}}(t) \in {\mathbb{L}} \}
    \end{split}
    \label{eqn:backward_reach_dist}
\end{equation}
Similarly, the forward reachable tube (FRT) is defined as the set of all states which can be reached by an agent starting from an initial set ${\mathbb{L}}$ while following the dynamics $f_r({\bm{s}}, {\bm{a}}, {\bm{d}})$.

In this work, the learnt system dynamics $\hat{f}: {\mathbb{S}} \times {\mathbb{A}} \rightarrow {\mathbb{S}}$ is represented using an artificial neural network (ANN).
Let the ANN have $H$ number of hidden layers.
Then, the output from a layer is given as
\begin{equation}
    {\bm{y}}^{(h)} = g_h ({\bm{W}}^{(h)} {\bm{y}}^{(h-1)} + {\bm{\theta}}^{(h)}), \; h = 1, 2, \hdots H,
\end{equation}
where ${\bm{W}}^{(h)}$ is the weight matrix mapping the output from layer $(h-1)$ to the input of layer $h$ and ${\bm{\theta}}^{(h)}$ is the bias vector, of appropriate dimensions.
The first layer ${\bm{y}}^{(0)} = x$.
In this work, the following two assumptions need to be made for the system dynamics and thereby the ANN
\begin{assumption} (\textit{Bounded output}) \label{ass1} $\forall i = 1, 2, \hdots n, \; \hat{o}_i \in [a, \; b]$ where $-\infty < a, \; b < \infty$.
\end{assumption}
\begin{assumption} \label{ass2}
The ANN is uniformly continuous and Lipschitz continuous.
\end{assumption}
Assumption \ref{ass1} can be easily satisfied by choosing a bounded activation function for the output layer.
$tanh()$ and $sigmoid()$ are the two most commonly used activation functions to bound the outputs between $[-1, \; 1]$ and $[0, \; 1]$.
Therefore, to bound each element in $\hat{{\bm{o}}} \in \mathbb{R}^n$ between some range $[a, \; b]$ or $[0, \; b]$ the following post-processing step can be conducted
\begin{equation}
\hat{o}_i = a + \frac{(b-a)}{2}(y^{(H)}_i + 1), \; \; \; \; \; \; \hat{o}_i = b y^{(H)}_i, \; \; \; \forall i = 1, 2, \hdots n,
\end{equation}
where ${\bm{y}}^{(H)} = tanh({\bm{W}}^{(H)} {\bm{y}}^{(H-1} + {\bm{\theta}}^{(H)})$ or ${\bm{y}}^{(H)} = sigmoid({\bm{W}}^{(H)} {\bm{y}}^{(H-1} + {\bm{\theta}}^{(H)})$ respectively.

\begin{problem} \label{prob1}
Given a deterministic model-based controller $\pi(s)$, determine the subset of safe initial states $\mathcal{S}_{safe} \subset \mathcal{S}_0$ such that any finite time trajectory $\tau$ starting from an initial state $s_0 \in \mathcal{S}_{safe}$ and ending at a final state $s_f \in \mathcal{T}$, following policy $\pi(s)$, satisfies the state constraint $\mathcal{C}_s$ imposed on the system despite the disturbance $d$.
\end{problem}
\end{comment}
\section{Experiments}
\begin{figure}[t]
\begin{center}
\includegraphics[width=0.32\textwidth]{images/2D_initial_problem.png}
\includegraphics[width=0.32\textwidth]{images/FRT_300_100_v2.png}
\includegraphics[width=0.32\textwidth]{images/BRT_300_100_v2.png}
\end{center}
\caption{(Left) The environment for the safe land navigation problem. (Middle) The FRTs for both the augmented learnt model dynamics and true system dynamics classify the controller as unsafe. (Right) The BRT computed from obstacles 1 and 2 for the given controller. ${\mathbb{S}}_{safe}$ computed by the proposed BRT algorithm is compared with the ground truth (G.T.) data for the safe initial states, marked in green.}
\label{fig:2D_init}
\end{figure}

In our experiments, we aim to answer the following questions:  \textbf{(a)} Can safety verification be done for an NN-based $\pi$ and $\hat f_\theta$ using FRT?, and \textbf{(b)} Can ${\mathbb{S}}_\text{safe}$ be identified using BRT if $\pi$ is deemed unsafe?
To answer these two questions, we demonstrate results on the following  domains inspired by real world safety-critical problems, where RL controllers developed using a learnt model can be appealing, as they can adapt to transition dynamics involving friction, air-drag, wind, etc., which might be hard to explicitly model otherwise.

\noindent \textbf{Safe land navigation: } Navigation of a ground robot in an indoor environment is a common application which requires the satisfaction of safety constraints by $\pi$ to avoid collision with obstacles.
For this setting, we simulate a ground robot which has continuous states and actions.
The initial configuration of the domain is shown in Fig. \ref{fig:2D_init}.
The set of initial and goal states are represented by circles and the obstacles with rectangles.

\noindent \textbf{Safe aerial navigation: }
This domain simulates a navigation problem in an urban environment for an unmanned aerial vehicle (UAV).
Constraints are incorporated while training $\pi$ to ensure that collision is avoided with potential obstacles in its path.
States and actions are both continuous and the initial configuration of the domain is shown in Fig \ref{fig:3D_init}.
The set of initial and goal states are represented using cuboids, and the obstacles with cylinders.

\paragraph{Analysis:} To address the questions with respect to the above mentioned domains, we first train a NN based $\hat f_\theta$ to estimate the dynamics using sampled transitions.
This $\hat f_\theta$ is then also used to learn a NN based controller $\pi$ which is trained with a cost function designed to mitigate collisions.
For brevity, only the representative results for this $\pi$ are discussed here; implementation details and more experimental results are available in Appendix \ref{sec:modelRL}, \ref{sec:implementation}, \ref{sec:results_land} and \ref{sec:results_aerial}.

To address the first question, the FRT is computed for both the domains over the augmented learnt dynamics $\hat{f}^{(r)}_{\bm{\theta}}$ as shown in Fig. \ref{fig:2D_init} and Fig. \ref{fig:3D_init}.
Additionally, for land navigation we also compute the FRT over the true system dynamics $f$, which serves as a way to validate the safety verification result of $\pi$ from the proposed framework.
It is observed that for both the domains, FRTs deem the given policy $\pi$ as \textit{unsafe}, since the FRTs intersect with one of the obstacles.
Even when $\pi$ is learnt using a cost function designed to avoid collisions, the proposed safety verification framework successfully brings out the limitations of $\pi$, which may have resulted due to the use of function approximations, ill-specified hyper-parameters, convergence to a local optimum, etc.

For the second question, the BRT is computed from both the obstacles for the given controller $\pi$, as shown in Fig. \ref{fig:2D_init} and Fig. \ref{fig:3D_init}.
To estimate the accuracy of the BRT computation, we compare the computed ${\mathbb{S}}_{unsafe}$ and ${\mathbb{S}}_{safe}$ sets with the ground truth (G.T.) data generated using random samples of possible trajectories.
It is observed that the BRT from obstacle 1 does not intersect with ${\mathbb{S}}_0$, implying that all trajectories are safe w.r.t. obstacle 1.
However, the BRT from obstacle 2 intersects with ${\mathbb{S}}_0$ and identifies the subset of initial states which are unsafe.
Such an information can be critical to safely deploy even an unsafe controller just by restricting its starting conditions.

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.32\textwidth]{images/3D_init.png}
\includegraphics[width=0.32\textwidth]{images/FRT_3D_pic7.png}
\includegraphics[width=0.32\textwidth]{images/BRT_3D_pic1.png}
\end{center}
\caption{(Left) The safe aerial navigation domain. (Middle) The FRT is computed for the augmented learnt dynamics. For clarity in 3D, FRT with true dynamics is not plotted; instead the approximation quality can be better visualized in BRT. (Right) Comparison of ${\mathbb{S}}_{safe}$ computed by the proposed BRT algorithm with the ground truth data of the safe initial states, marked in green  (and unsafe states marked in red), it is observed that the BRT conservatively approximates  ${\mathbb{S}}_{safe}$.}
\label{fig:3D_init}
\end{figure}

\begin{comment}
\subsection{Experiment 2:}
The second problem is a 3D navigation problem for an aerial vehicle.
We consider an abstraction of the problem of flying an unmanned aerial vehicle (UAV) in an urban environment.
Simple point mass dynamics is considered for the aerial vehicle with the state vector consisting of the 3 dimensional position of the UAV, ${\bm{s}} = [x, \; y, \; z]^T$.
The action vector ${\bm{a}} = [v, \; \theta, \; \gamma]^T \in R^3$ consists of the velocity, the pitch angle and the heading angle.
The proposed algorithm can handle complex dynamics for a UAV, however, for the purpose of illustration we consider 3 dimensional point mass dynamics.
In the urban environment, the buildings form the obstacles for the UAV and the policy $\pi(s)$ is trained to reach a goal state and altitude while avoiding these obstacles.
\end{comment}
\section{Experiments}
\section{Conclusion}
In this paper, we have presented a novel framework using forward and backward reachable tubes for safety verification and determination of the subset of initial states for which a given model-based RL controller always satisfies the state constraints.
 The main contribution of this work is the formulation of the reachability problem for a neural network modeled system dynamics and the use of level set method to compute an exact reachable tube by solving the Hamilton-Jacobi partial differential equation, for the reinforcement learning framework, thereby minimizing approximation errors that other existing reachability methods suffer.
 Additionally, the proposed framework can identify the set of safe initial sets for a given policy, thereby determining the initial conditions for which even a sub-optimal, unsafe policy satisfies the safety constraints.

\begin{comment}
\section{Submission of conference papers to ICLR 2021}

ICLR requires electronic submissions, processed by
\url{https://openreview.net/}. See ICLR's website for more instructions.

If your paper is ultimately accepted, the statement {\tt
  {\textbackslash}iclrfinalcopy} should be inserted to adjust the
format to the camera ready requirements.

The format for the submissions is a variant of the NeurIPS format.
Please read carefully the instructions below, and follow them
faithfully.

\subsection{Style}

Papers to be submitted to ICLR 2021 must be prepared according to the
instructions presented here.

Authors are required to use the ICLR \LaTeX{} style files obtainable at the
ICLR website. Please make sure you use the current files and
not previous versions. Tweaking the style files may be grounds for rejection.

\subsection{Retrieval of style files}

The style files for ICLR and other conference information are available online at:
\begin{center}
   \url{http://www.iclr.cc/}
\end{center}
The file \verb+iclr2021_conference.pdf+ contains these
instructions and illustrates the
various formatting requirements your ICLR paper must satisfy.
Submissions must be made using \LaTeX{} and the style files
\verb+iclr2021_conference.sty+ and \verb+iclr2021_conference.bst+ (to be used with \LaTeX{}2e). The file
\verb+iclr2021_conference.tex+ may be used as a ``shell'' for writing your paper. All you
have to do is replace the author, title, abstract, and text of the paper with
your own.

The formatting instructions contained in these style files are summarized in
sections \ref{gen_inst}, \ref{headings}, and \ref{others} below.

\section{General formatting instructions}
\label{gen_inst}

The text must be confined within a rectangle 5.5~inches (33~picas) wide and
9~inches (54~picas) long. The left margin is 1.5~inch (9~picas).
Use 10~point type with a vertical spacing of 11~points. Times New Roman is the
preferred typeface throughout. Paragraphs are separated by 1/2~line space,
with no indentation.

Paper title is 17~point, in small caps and left-aligned.
All pages should start at 1~inch (6~picas) from the top of the page.

Authors' names are
set in boldface, and each name is placed above its corresponding
address. The lead author's name is to be listed first, and
the co-authors' names are set to follow. Authors sharing the
same address can be on the same line.

Please pay special attention to the instructions in section \ref{others}
regarding figures, tables, acknowledgments, and references.

There will be a strict upper limit of 8 pages for the main text of the initial submission, with unlimited additional pages for citations. Note that the upper page limit differs from last year!Authors may use as many pages of appendices (after the bibliography) as they wish, but reviewers are not required to read these. During the rebuttal phase and for the camera ready version, authors are allowed one additional page for the main text, for a strict upper limit of 9 pages.

\section{Headings: first level}
\label{headings}

First level headings are in small caps,
flush left and in point size 12. One line space before the first level
heading and 1/2~line space after the first level heading.

\subsection{Headings: second level}

Second level headings are in small caps,
flush left and in point size 10. One line space before the second level
heading and 1/2~line space after the second level heading.

\subsubsection{Headings: third level}

Third level headings are in small caps,
flush left and in point size 10. One line space before the third level
heading and 1/2~line space after the third level heading.

\section{Citations, figures, tables, references}
\label{others}

These instructions apply to everyone, regardless of the formatter being used.

\subsection{Citations within the text}

Citations within the text should be based on the \texttt{natbib} package
and include the authors' last names and year (with the ``et~al.'' construct
for more than two authors). When the authors or the publication are
included in the sentence, the citation should not be in parenthesis using \verb|\citet{}| (as
in ``See \citet{Hinton06} for more information.''). Otherwise, the citation
should be in parenthesis using \verb|\citep{}| (as in ``Deep learning shows promise to make progress
towards AI~\citep{Bengio+chapter2007}.'').

The corresponding references are to be listed in alphabetical order of
authors, in the \textsc{References} section. As to the format of the
references themselves, any style is acceptable as long as it is used
consistently.

\subsection{Footnotes}

Indicate footnotes with a number\footnote{Sample of the first footnote} in the
text. Place the footnotes at the bottom of the page on which they appear.
Precede the footnote with a horizontal rule of 2~inches
(12~picas).\footnote{Sample of the second footnote}

\subsection{Figures}

All artwork must be neat, clean, and legible. Lines should be dark
enough for purposes of reproduction; art work should not be
hand-drawn. The figure number and caption always appear after the
figure. Place one line space before the figure caption, and one line
space after the figure. The figure caption is lower case (except for
first word and proper nouns); figures are numbered consecutively.

Make sure the figure caption does not get separated from the figure.
Leave sufficient space to avoid splitting the figure and figure caption.

You may use color figures.
However, it is best for the
figure captions and the paper body to make sense if the paper is printed
either in black/white or in color.
\begin{figure}[h]
\begin{center}
\fbox{\rule[-.5cm]{0cm}{4cm} \rule[-.5cm]{4cm}{0cm}}
\end{center}
\caption{Sample figure caption.}
\end{figure}

\subsection{Tables}

All tables must be centered, neat, clean and legible. Do not use hand-drawn
tables. The table number and title always appear before the table. See
Table~\ref{sample-table}.

Place one line space before the table title, one line space after the table
title, and one line space after the table. The table title must be lower case
(except for first word and proper nouns); tables are numbered consecutively.

\begin{table}[t]
\caption{Sample table title}
\label{sample-table}
\begin{center}
\begin{tabular}{ll}
\multicolumn{1}{c}{\bf PART}  &\multicolumn{1}{c}{\bf DESCRIPTION}
\\ \hline \\
Dendrite         &Input terminal \\
Axon             &Output terminal \\
Soma             &Cell body (contains cell nucleus) \\
\end{tabular}
\end{center}
\end{table}

\section{Default Notation}

In an attempt to encourage standardized notation, we have included the
notation file from the textbook, \textit{Deep Learning}
\cite{goodfellow2016deep} available at
\url{https://github.com/goodfeli/dlbook_notation/}.  Use of this style
is not required and can be disabled by commenting out
\texttt{math\_commands.tex}.

\centerline{\bf Numbers and Arrays}
\bgroup
\def1.5{1.5}
\begin{tabular}{p{1in}p{3.25in}}
$\displaystyle a$ & A scalar (integer or real)\\
$\displaystyle {\bm{a}}$ & A vector\\
$\displaystyle {\bm{A}}$ & A matrix\\
$\displaystyle {\tens{A}}$ & A tensor\\
$\displaystyle {\bm{I}}_n$ & Identity matrix with $n$ rows and $n$ columns\\
$\displaystyle {\bm{I}}$ & Identity matrix with dimensionality implied by context\\
$\displaystyle {\bm{e}}^{(i)}$ & Standard basis vector $[0,\dots,0,1,0,\dots,0]$ with a 1 at position $i$\\
$\displaystyle \text{diag}({\bm{a}})$ & A square, diagonal matrix with diagonal entries given by ${\bm{a}}$\\
$\displaystyle {\textnormal{a}}$ & A scalar random variable\\
$\displaystyle {\mathbf{a}}$ & A vector-valued random variable\\
$\displaystyle {\mathbf{A}}$ & A matrix-valued random variable\\
\end{tabular}
\egroup
\vspace{0.25cm}

\centerline{\bf Sets and Graphs}
\bgroup
\def1.5{1.5}

\begin{tabular}{p{1.25in}p{3.25in}}
$\displaystyle {\mathbb{A}}$ & A set\\
$\displaystyle \mathbb{R}$ & The set of real numbers \\
$\displaystyle \{0, 1\}$ & The set containing 0 and 1 \\
$\displaystyle \{0, 1, \dots, n \}$ & The set of all integers between $0$ and $n$\\
$\displaystyle [a, b]$ & The real interval including $a$ and $b$\\
$\displaystyle (a, b]$ & The real interval excluding $a$ but including $b$\\
$\displaystyle {\mathbb{A}} \backslash {\mathbb{B}}$ & Set subtraction, i.e., the set containing the elements of ${\mathbb{A}}$ that are not in ${\mathbb{B}}$\\
$\displaystyle {\mathcal{G}}$ & A graph\\
$\displaystyle \parents_{\mathcal{G}}({\textnormal{x}}_i)$ & The parents of ${\textnormal{x}}_i$ in ${\mathcal{G}}$
\end{tabular}
\vspace{0.25cm}

\centerline{\bf Indexing}
\bgroup
\def1.5{1.5}

\begin{tabular}{p{1.25in}p{3.25in}}
$\displaystyle {a}_i$ & Element $i$ of vector ${\bm{a}}$, with indexing starting at 1 \\
$\displaystyle {a}_{-i}$ & All elements of vector ${\bm{a}}$ except for element $i$ \\
$\displaystyle {A}_{i,j}$ & Element $i, j$ of matrix ${\bm{A}}$ \\
$\displaystyle {\bm{A}}_{i, :}$ & Row $i$ of matrix ${\bm{A}}$ \\
$\displaystyle {\bm{A}}_{:, i}$ & Column $i$ of matrix ${\bm{A}}$ \\
$\displaystyle {\etens{A}}_{i, j, k}$ & Element $(i, j, k)$ of a 3-D tensor ${\tens{A}}$\\
$\displaystyle {\tens{A}}_{:, :, i}$ & 2-D slice of a 3-D tensor\\
$\displaystyle {\textnormal{a}}_i$ & Element $i$ of the random vector ${\mathbf{a}}$ \\
\end{tabular}
\egroup
\vspace{0.25cm}

\centerline{\bf Calculus}
\bgroup
\def1.5{1.5}
\begin{tabular}{p{1.25in}p{3.25in}}
$\displaystyle\frac{d y} {d x}$ & Derivative of $y$ with respect to $x$\\ [2ex]
$\displaystyle \frac{\partial y} {\partial x} $ & Partial derivative of $y$ with respect to $x$ \\
$\displaystyle \nabla_{\bm{x}} y $ & Gradient of $y$ with respect to ${\bm{x}}$ \\
$\displaystyle \nabla_{\bm{X}} y $ & Matrix derivatives of $y$ with respect to ${\bm{X}}$ \\
$\displaystyle \nabla_{\tens{X}} y $ & Tensor containing derivatives of $y$ with respect to ${\tens{X}}$ \\
$\displaystyle \frac{\partial f}{\partial {\bm{x}}} $ & Jacobian matrix ${\bm{J}} \in \mathbb{R}^{m\times n}$ of $f: \mathbb{R}^n \rightarrow \mathbb{R}^m$\\
$\displaystyle \nabla_{\bm{x}}^2 f({\bm{x}})\text{ or }{\bm{H}}( f)({\bm{x}})$ & The Hessian matrix of $f$ at input point ${\bm{x}}$\\
$\displaystyle \int f({\bm{x}}) d{\bm{x}} $ & Definite integral over the entire domain of ${\bm{x}}$ \\
$\displaystyle \int_{\mathbb{S}} f({\bm{x}}) d{\bm{x}}$ & Definite integral with respect to ${\bm{x}}$ over the set ${\mathbb{S}}$ \\
\end{tabular}
\egroup
\vspace{0.25cm}

\centerline{\bf Probability and Information Theory}
\bgroup
\def1.5{1.5}
\begin{tabular}{p{1.25in}p{3.25in}}
$\displaystyle P({\textnormal{a}})$ & A probability distribution over a discrete variable\\
$\displaystyle p({\textnormal{a}})$ & A probability distribution over a continuous variable, or over
a variable whose type has not been specified\\
$\displaystyle {\textnormal{a}} \sim P$ & Random variable ${\textnormal{a}}$ has distribution $P$\\% so thing on left of \sim should always be a random variable, with name beginning with \r
$\displaystyle  \mathbb{E}_{{\textnormal{x}}\sim P} [ f(x) ]\text{ or } \mathbb{E} f(x)$ & Expectation of $f(x)$ with respect to $P({\textnormal{x}})$ \\
$\displaystyle \mathrm{Var}(f(x)) $ &  Variance of $f(x)$ under $P({\textnormal{x}})$ \\
$\displaystyle \mathrm{Cov}(f(x),g(x)) $ & Covariance of $f(x)$ and $g(x)$ under $P({\textnormal{x}})$\\
$\displaystyle H({\textnormal{x}}) $ & Shannon entropy of the random variable ${\textnormal{x}}$\\
$\displaystyle D_{\mathrm{KL}} ( P \Vert Q ) $ & Kullback-Leibler divergence of P and Q \\
$\displaystyle \mathcal{N} ( {\bm{x}} ; {\bm{\mu}} , {\bm{\Sigma}})$ & Gaussian distribution %
over ${\bm{x}}$ with mean ${\bm{\mu}}$ and covariance ${\bm{\Sigma}}$ \\
\end{tabular}
\egroup
\vspace{0.25cm}

\centerline{\bf Functions}
\bgroup
\def1.5{1.5}
\begin{tabular}{p{1.25in}p{3.25in}}
$\displaystyle f: {\mathbb{A}} \rightarrow {\mathbb{B}}$ & The function $f$ with domain ${\mathbb{A}}$ and range ${\mathbb{B}}$\\
$\displaystyle f \circ g $ & Composition of the functions $f$ and $g$ \\
  $\displaystyle f({\bm{x}} ; {\bm{\theta}}) $ & A function of ${\bm{x}}$ parametrized by ${\bm{\theta}}$.
  (Sometimes we write $f({\bm{x}})$ and omit the argument ${\bm{\theta}}$ to lighten notation) \\
$\displaystyle \log x$ & Natural logarithm of $x$ \\
$\displaystyle \sigma(x)$ & Logistic sigmoid, $\displaystyle \frac{1} {1 + \exp(-x)}$ \\
$\displaystyle \zeta(x)$ & Softplus, $\log(1 + \exp(x))$ \\
$\displaystyle || {\bm{x}} ||_p $ & $L^p$ norm of ${\bm{x}}$ \\
$\displaystyle || {\bm{x}} || $ & $L^2$ norm of ${\bm{x}}$ \\
$\displaystyle x^+$ & Positive part of $x$, i.e., $\max(0,x)$\\
$\displaystyle \bm{1}_\mathrm{condition}$ & is 1 if the condition is true, 0 otherwise\\
\end{tabular}
\egroup
\vspace{0.25cm}

\section{Final instructions}
Do not change any aspects of the formatting parameters in the style files.
In particular, do not modify the width or length of the rectangle the text
should fit into, and do not change font sizes (except perhaps in the
\textsc{References} section; see below). Please note that pages should be
numbered.

\section{Preparing PostScript or PDF files}

Please prepare PostScript or PDF files with paper size ``US Letter'', and
not, for example, ``A4''. The -t
letter option on dvips will produce US Letter files.

Consider directly generating PDF files using \verb+pdflatex+
(especially if you are a MiKTeX user).
PDF figures must be substituted for EPS figures, however.

Otherwise, please generate your PostScript and PDF files with the following commands:
\begin{verbatim}
dvips mypaper.dvi -t letter -Ppdf -G0 -o mypaper.ps
ps2pdf mypaper.ps mypaper.pdf
\end{verbatim}

\subsection{Margins in LaTeX}

Most of the margin problems come from figures positioned by hand using
\verb+\special+ or other commands. We suggest using the command
\verb+\includegraphics+
from the graphicx package. Always specify the figure width as a multiple of
the line width as in the example below using .eps graphics
\begin{verbatim}
   \usepackage[dvips]{graphicx} ...
   \includegraphics[width=0.8\linewidth]{myfile.eps}
\end{verbatim}
or
\begin{verbatim}
   \usepackage[pdftex]{graphicx} ...
   \includegraphics[width=0.8\linewidth]{myfile.pdf}
\end{verbatim}
for .pdf graphics.
See section~4.4 in the graphics bundle documentation (\url{http://www.ctan.org/tex-archive/macros/latex/required/graphics/grfguide.ps})

A number of width problems arise when LaTeX cannot properly hyphenate a
line. Please give LaTeX hyphenation hints using the \verb+\-+ command.

\subsubsection*{Author Contributions}
If you'd like to, you may include  a section for author contributions as is done
in many journals. This is optional and at the discretion of the authors.

\subsubsection*{Acknowledgments}
Use unnumbered third level headings for the acknowledgments. All
acknowledgments, including those to funding agencies, go at the end of the paper.
\end{comment}

\section{Introduction}
\begin{comment}
Reinforcement learning (RL) paradigm has emerged as a powerful tool for developing optimal controllers for sequential decision making problems.
Particularly, RL algorithms don't assume \textit{a priori} knowledge of the system dynamics.
Instead, a data driven model of the system dynamics is developed either implicitly or explicitly.
This implies that any uncertainties in the system dynamics, like environment noise, friction etc., are also learnt by the modeled dynamics.
In this work, we narrow our focus on model-based RL algorithms with applications to safety critical physical systems.
\end{comment}

One of the primary reasons for the growing application of reinforcement learning (RL) algorithms in developing optimal controllers is that RL does not assume \textit{a priori} knowledge of the system dynamics.
Model-based RL explicitly learns a model of the system dynamics, from observed samples of state transitions.
This learnt model is used along with a planning algorithm to develop optimal controllers for different tasks.
Thus, any uncertainties in the system, including environment noise, friction, air-drag etc., can also be captured by the modeled dynamics.

However, the performance of the controller is directly related to how accurately the learnt model represents the true system dynamics.
Due to the discrepancy between the learnt model and the true model, the developed controller can behave unexpectedly when deployed on the real physical system, e.g., land robots, UAVs, etc. \citep{benbrahim1997biped, endo2008learning, morimoto2001acquisition}.
This unexpected behavior may result in the violation of constraints imposed on the system, thereby violating its safety requirements \citep{moldovan2012safe}.
Thus, it is necessary to have a framework which can ensure that the controller will satisfy the safety constraints \textit{before} it is deployed on a real system.
This raises the primary question of interest: \textit{Given a set of safety constraints imposed on the state space, how do we determine whether a given controller is safe or not?}

In the literature, there have been several works that focus on the problem of ensuring safety.
Some of the works aim at \textit{finding} a safe controller, under the assumption of a known baseline safe policy \citep{hans2008safe,garcia2012safe,berkenkamp2017safe,thomas2015high,laroche2019safe}, or several known safe policies \citep{perkins2002lyapunov}.
However, such safe policies may not be readily available in general.
Alternatively,  \citet{akametalu2014reachability} used reachability analysis to develop safe model-based controllers, under the assumption that the system dynamics can be modeled using Gaussian processes, i.e., an assumption which is violated by most modern RL methods that make use of neural networks (NN) instead.
While there have been several works proposed to develop safe controllers, some of the assumptions made in these works may not be possible to realize in practice.
In recent years, this limitation has drawn attention towards developing \textit{verification} frameworks for RL controllers, which is the focus of this paper.
Since verifying safety of an NN based RL controller is also related to verifying the safety of the underlying NN model \citep{xiang2018reachable, tran2019parallelizable, xiang2018output, tran2019safety}, we provide an additional review for these methods in Appendix \ref{sec:ERW}.

\paragraph{Contributions: }
In this work, we focus on the problem of determining whether a given controller is safe or not, with respect to satisfying constraints imposed on the state space.
To do so, we propose a novel safety verification algorithm for model-based RL controllers using \textit{forward reachable tube} analysis that can handle NN based learnt dynamics and controllers, while also being robust against modeling error.
The problem of determining the reachable tube is framed as an optimal control problem using the Hamilton Jacobi (HJ) partial differential equation (PDE), whose solution is computed using the level set method.
The advantage of using the level set method is the fact that it can represent sets with non-convex boundaries, thereby avoiding approximation errors that most existing methods suffer from.
Additionally, if a controller is deemed unsafe, we take a step further to identify if there are any starting conditions from which the given controller can be safely executed.
To achieve this, a \textit{backward reachable tube} is computed for the learnt model and, to the best of our knowledge, this is the first work which computes the backward reachable tube over an NN.
Finally, empirical results are presented on two domains inspired by real-world applications where safety verification is critical.

\begin{comment}
There are two problems to consider when developing a safety verification framework for a  model-based RL algorithm.
The first problem deals with the type of representation of the learnt model.
Artificial neural networks (ANNs) are a popular choice of function approximators used to represent the learnt model because of their property of learning complex non-linear mappings
However, their mathematical formulation is not favorable to using analytical tools for verifying the properties.
Secondly, since the true system dynamics is unavailable, the framework proposed in this work requires a conservative estimate of the modeling error.
This limitation is shared with other notable works on safety verification of RL controllers.
This implies that any safety verification algorithm should be developed to accommodate uncertainties in the learnt dynamics.

To achieve the above goals, we propose to use reachable set analysis
\end{comment}

\section{Safety Verification}
In this section, we first present the problem setting and specify the assumptions made in this work, followed by the proposed algorithm.

\subsection{Model-Based Reinforcement Learning}

\subsection{Reachable Set Analysis}
In this section we give a brief overview of the HJI PDE formulation required to compute the BRT for a non-linear system with uncertain parameters.

Before formally defining a BRT, we introduce the problem setting.
The goal is to determine the set of states which can reach a target set $\mathcal{L} \subset \mathbb{R}^n$ under finite time.
While computing the BRT, convention is that the target set is reached at final time $t_f = 0$ and the finite time trajectory starts at some initial time $t_0 = -T, \; 0 < T < \infty$.
The system follows the dynamics $f(s, a)$ and its initial states are represented as $\mathcal{S}_0 \in \mathcal{S}$.
Additionally, there might be a disturbance in the system or the environment which affects the system state as well.
This disturbance is represented using a finite set $\mathcal{D} \subset \mathbb{R}^n$ and is considered as a control input which can affect the progression of the system states in an undesired manner.
Let $f_r: \mathcal{S} \times \mathcal{A} \times \mathcal{D} \rightarrow \mathcal{S}$ denote the augmented system dynamics which also takes into account the disturbance input $d \in \mathcal{D}$ to determine the next state.
\begin{assumption}
The control inputs $a$ and $d$ are drawn from closed and bounded sets, i.e., finite sets $\mathcal{A} \subset \mathbb{R}^m$ and $\mathcal{D} \subset \mathbb{R}^m$ respectively.
\end{assumption}
Besides the assumptions made in this work, we explicitly state the assumptions from \cite{TheHolyOne, mitchell2005time} as requirements for Theorem \ref{thm1} stated below and later prove how our assumptions fulfill these requirements to leverage the results from previous works.
\begin{reqm} \label{req1}
The system  dynamics $f_r(s, a, d)$ is uniformly continuous, bounded and Lipschitz continuous in $s$ while the inputs $a$ and $d$ are fixed.
\end{reqm}

\begin{reqm} \label{req2}
The target set $\mathcal{L} \subset \mathbb{R}^n$ is closed and can be represented as the zero sublevel set of a bounded and Lipschitz continuous function $l: \mathbb{R}^n \rightarrow \mathbb{R}$ as
\begin{equation}
    \mathcal{L} = \{ s: l(s) \le 0\}.
    \label{eqn:target_level_func}
\end{equation}
\end{reqm}
Thus, (\ref{eqn:target_level_func}) determines a way to check whether a state lies inside or outside the target set $\mathcal{L}$.
If $\mathcal{L}$ can be represented using a regular, well-defined geometric shape (like a sphere, rectangle, cylinder etc.) then deriving the level set function $l(s)$ is straight forward.
An irregularly shaped $\mathcal{L}$ can be represented as a union of several well-defined geometric shapes to derive $l(s)$.

From the above information, we formally define the BRT for a finite time $T$, $R_B(T)$ as the set of states that can reach $\mathcal{L}$ within time $T$, despite the disturbance $d$ in the system.
Specifically, if $\mathcal{L}$ is an unsafe target set, we are interested in determining the subset of initial states $\mathcal{S}_{unsafe} \subset \mathcal{S}_0$ such that for all trajectories starting from an initial state $s_0 \in \mathcal{S}_{unsafe}$ there exists a sequence of actions $a(\cdot)$ which will eventually drive the system following dynamics $f_r(s, a, d)$ to $\mathcal{L}$ despite the disturbance.
\begin{equation}
\begin{split}
    R_B(T) := \{ s_0: \exists a, \; s.t. &\; \forall d, \; s(\cdot) \; satisfies \; \dot{s}=f_r(s,a, d), \\
    s(0)=s_0; &\; \exists t \in [-T, 0], \; s.t. \; s(t) \in \mathcal{L} \}
    \end{split}
    \label{eqn:backward_reach_dist}
\end{equation}

The above goal of computing $R_B(T)$ can be achieved by formulating the BRT as an optimal control problem.
From (\ref{eqn:target_level_func}), a trajectory is said to reach $\mathcal{L}$ if it's state at the final time step $t_f = 0$ satisfies $l(s(0)) \le 0$.
Similarly, if $l(s(0)) > 0$ then the trajectory never reaches $\mathcal{L}$.
Therefore, to formulate the BRT problem as an optimal control problem a function $\phi$ is defined as
\begin{equation}
\begin{split}
    \phi(s(t), t) &= \max_a \min_d l(s(0)), \; t \in [-T, 0] \\
    \text{s.t. } \; \; \dot{s} &= f_r(s, a, d) \\
    l&(s(0)) \le 0
\end{split}
\label{eqn:game_opt}
\end{equation}
where the control input $a \in \mathcal{A}$ is driving the system away from the unsafe target set $\mathcal{L}$ by maximizing the value of $l(s(0))$ at the final time step and the disturbance control is doing the opposite.
The above formulation is a 2 player game where player 1 is the controller and player 2 is the disturbance.
The choice of formulating the problem as a maximin optimization gives the advantage to player 2 \cite{TheHolyOne}, thereby ensuring that a worst case scenario is considered for player 1 and the BRT is over approximated to give a conservative solution.
It should be noted that depending on whether the target set is safe or unsafe, the optimization formulation in (\ref{eqn:game_opt}) can be modified according to the table \ref{tab:game_opt}.
\begin{table}
\begin{center}
\caption{Optimization objectives of players for safe and unsafe target sets.} \label{tab:game_opt}
\begin{tabular}{| c | c | c | }
  \hline
  & Unsafe target set & Safe target set \\
  \hline
 Player 1 & max $\phi$ & min $\phi$ \\
 \hline
 Player 2 & min $\phi$ & max $\phi$ \\
 \hline
\end{tabular}
\end{center}
\end{table}
\begin{reqm} \label{ass3}
The input signals are drawn from the following sets
\begin{equation}
    \begin{split}
        a(\cdot) \in \mathcal{A}(t) \overset{\Delta}{=} \{ \eta:[-T, \; 0] \rightarrow \mathcal{A}| \eta(\cdot) \text{ is measurable} \}\\
        d(\cdot) \in \mathcal{D}(t) \overset{\Delta}{=} \{ \eta:[-T, \; 0] \rightarrow \mathcal{D}| \eta(\cdot) \text{ is measurable} \}
    \end{split}
\end{equation}
where $\mathcal{A} \subset \mathbb{R}^m$ and $\mathcal{D} \subset \mathbb{R}^m$ are compact and $T > 0$.
\end{reqm}

We now state theorem 3 from \cite{TheHolyOne} or equivalently theorem 2 from \cite{mitchell2005time} because it is essential to validating the algorithm proposed in the next section.

\begin{theorem} \label{thm1} Under requirements 1-3, let $\phi: \mathbb{R}^n \times [-T, \; 0] \rightarrow \mathbb{R}$ be the viscosity solution of the terminal value HJI PDE
\begin{equation}
\begin{split}
    D_t \phi(s, t) + \min[0, H(s, \; D_s \phi(s, t))] &= 0, \\
    \phi(s, 0) &= l(s) \\
\text{where,} \; H(s, p) &= \max_{a \in \mathcal{A}} \min_{d \in \mathcal{D}} p^T f_r(s, a, d).
\end{split}
\label{eqn:PDE_HJI}
\end{equation}
Then, the zero sublevel set of $\phi$ describes $B_R(t)$
\begin{equation}
    B_R(t) = \{s \in \mathbb{R}^n | \phi(s, t) \le 0 \}.
\label{eqn:BRT}
\end{equation}
\end{theorem}
\begin{proof} See proof of theorem 3 in \cite{TheHolyOne} or theorem 2 in \cite{mitchell2005time}.
\end{proof}

The formulation in (\ref{eqn:PDE_HJI}) is key to computing the BRT and a numerical solution for $\phi$ using the level set methods gives an implicit representation of $B_R(t)$ in (\ref{eqn:BRT}).

\textbf{Forward Reachable Set (FRS):}
The HJI formulation can also be used to compute the forward reachable set by forming an initial value PDE.
Using the same notation defined above, the FRS is defined as all the states which can be reached from an initial set $\mathcal{L}$ in finite time $T, \; 0< T < \infty$ despite the disturbance $d$.
\begin{equation}
\begin{split}
F_R(T) := \{s: \exists a, \text{ s.t. } \forall d, s(\cdot) \text{ satisfies } \dot{s} = f_r(s, a, d), \\
s(0) \in \mathcal{L}; s(T) = s \}
\end{split}
\end{equation}

Using a change of variables the terminal value PDE in (\ref{eqn:PDE_HJI}) can be converted into an initial value PDE and solved for the FRS.
The FRS has been used to analyze the output range of an ANN controller and can be used to determine whether or not a policy $\pi$ is safe for a given set of initial states $\mathcal{L}$.
\begin{equation}
    F_R(T) \cap \mathcal{C}_s = \begin{cases} = \emptyset
    &\mbox{if } \pi \text{ is safe} \\
    \neq \emptyset &\mbox{if } \pi \text{ is not safe} \end{cases}
\end{equation}
However, the FRS can only give a binary answer regarding the safety of the policy.
But by computing the BRS, we provide the set of initial states for which a given policy $\pi$ is safe.
Thus, by computing the BRS we can make a stronger statement regarding the policy.

\section{Methodology}
\label{sec:methodology}
This section presents the proposed algorithm and details the mathematical formulation required to compute the BRT for a general non-linear model of the learnt system dynamics.
\subsection{Overview}
The goal of the proposed algorithm is to verify the safety (i.e. constraint satisfaction) property of the given controller $\pi(s)$ by solving problem \ref{prob1}.
For given initial states $\mathcal{S}_0$ a policy $\pi(s)$ can be classified as either completely safe, completely unsafe or partially safe.
$\pi(s)$ is said to be completely safe if all the initial states are safe states and is said to be completely unsafe if none of the initial states are safe.
Mathematically, this can be verified by the following set operations
\begin{equation}
\mathcal{S}_{safe} \cap \mathcal{S}_0 = \begin{cases} \emptyset
    & \implies \pi(s) \mbox{ is completely unsafe} \\
      \mathcal{S}_0
    & \implies \pi(s) \mbox{ is completely safe.}
\end{cases}
\label{eqn:safe_policy}
\end{equation}
The third case is the most interesting one where there is a subset of states for which $\pi(s)$ is safe, i.e., $\mathcal{S}_{safe} \subset \mathcal{S}_0$ and hence it is classified as a partially safe policy.

Thus, to classify a policy we need to determine $\mathcal{S}_{safe}$.
In the proposed algorithm this is achieved by computing the BRT from the constrained region $\mathcal{C}_s$ to first determine the subset of unsafe states $\mathcal{S}_{unsafe} \subset \mathcal{S}_0$, followed by using simple set operations to determine $\mathcal{S}_{safe}$.
The formulation in (\ref{eqn:PDE_HJI}) is used to compute the BRT.
However, the true system dynamics is unknown and only data tuples $D$ observed from the environment are known.
Therefore, the learnt model $\hat{f}(s, a)$ will be used to construct the BRT.

Since the learnt model is imperfect, we first determine the modelling error to bound the error uncertainty in the model.
By analyzing the prediction error of the learnt model $\hat{f}(s, a)$ over the validating data set $D_v$, we identify certain properties (like mean error) to determine the error bounds on the model as $\underline{d} < d < \overline{d}$.
Subsequently, this error bound is used to construct the augmented model $\hat{f}_r(s, a, d)$ which is used to compute the BRT.
To achieve this, the modification in the mathematical formulation in (\ref{eqn:PDE_HJI}) is presented in the next subsection.
Once the BRT is computed, it's intersection with the set of initial states $\mathcal{S}_0$ is used to get the desired solution to problem $\ref{prob2}$.
The overview is given in algorithm \ref{algo2}.

\begin{algorithm}
\label{algo2}
 \caption{Safety verification of model-based RL policy (INCOMPLETE)}
\KwInput{$\pi, \hat{f}(s, a), \mathcal{C}_s, \mathcal{S}_0$}
 $[ \underline{d}, \; \overline{d}] \leftarrow EstimateModellingError(\hat{f}(s, a), D_v)$\;
 $B_C(-T) \leftarrow BRT(\hat{f}(s, a, d), \pi(s), \mathcal{C}_s)$\;
 $\mathcal{S}_{unsafe} \leftarrow B_C(-T) \cap \mathcal{S}_0$\;
 $\mathcal{S}_{safe} = \overline{\mathcal{S}}_{unsafe} \cap \mathcal{S}_0$
\end{algorithm}

\subsection{Proposed Algorithm}
The proposed algorithm is divided in 3 steps:
  \paragraph{Step 1: \textit{Determining modelling error}}

  Recall that the model-based controller $\pi(s)$ and the learnt model $\hat{f}(s, a)$ are already known.
  %
  The first task is to determine the augmented system dynamics $\hat{f}_r(s, a, d)$ where the disturbance input $d \in \mathcal{D} \subset \mathbb{R}^n$ represents the modelling error.
  %
  Using the additive formulation in (\ref{eqn:model_error}) the augmented system dynamics is represented as
  \begin{equation}
      \hat{f}_r(s, a, d) = \hat{f}(s, a) + d.
      \label{eqn:aug_dynamics}
  \end{equation}
  %
  To determine the bounds on the modelling error, the feedforward pass in (\ref{eqn:feedforward}) is conducted over the validation data set $D_v$ for the learnt model $\hat{f}(s, a)$ and the squared error between the predicted state and the desired state values is computed.
  %
  This squared error is computed separately for each of the $n$ output states of the learnt model and for all the data points in the validation set.
  %
  Analyzing the individual state prediction errors for every data point in $D_v$ gives metrics like the mean, variance, s.d. etc. to estimate the bounds of the modelling error such that $\underline{d} < d < \overline{d}$.
  %
  The exact choice of the metric depends on the problem setting and how conservative a bound is required.

  \paragraph{Step 2: \textit{Determine BRT for unsafe sets}}

  Instead of directly computing the set of safe initial states $\mathcal{S}_{safe}$ required in problem 1, a complementary problem is framed below.

  \begin{problem} \label{prob2}
  For a given unsafe set $\mathcal{C}_s$ determine the subset of unsafe initial states $\mathcal{S}_{unsafe} \subset \mathcal{S}_0$ such that for any initial state $s_0 \in \mathcal{S}_{unsafe}$ there exists at least one trajectory $\tau$ which will reach $\mathcal{C}_s$ under finite time while following the policy $\pi(s)$ under the disturbance $d$.
  \end{problem}

  To determine $\mathcal{S}_{unsafe}$ in problem \ref{prob2} the BRT from the unsafe target set $\mathcal{C}_s$ is constructed using (\ref{eqn:PDE_HJI}) and it is denoted as $B_{\mathcal{C}}(T)$.
  %
  Then $\mathcal{S}_{unsafe} = B_{\mathcal{C}}(T) \cap \mathcal{S}_0$.

  Since the true system dynamics is unknown, the formulation in (\ref{eqn:PDE_HJI}) has to be modified to substitute the system dynamics $f_r(s, a, d)$ with the learnt dynamics in (\ref{eqn:aug_dynamics}).
  %
  Additionally, the control input $a$ is now determined by the policy $\pi(s)$.
  %
  These changes are incorporated in the formulation (\ref{eqn:PDE_HJI}).

  The overall HJI PDE formulation remains the same, however the maximin Hamiltonian is modified to a constrained problem as
  \begin{equation}
\begin{split}
    H_{const}(s, p) &= \max_{a \in \mathcal{A}} \min_{d \in \mathcal{D}} p^T \hat{f}_r(s, a, d) \\
    \text{s.t. } &a = \pi(s)
\end{split}
\label{eqn:HJI_appx_con}
\end{equation}
Since the action $a$ is computed by the deterministic policy $\pi(s)$, the outer maximization optimization is carried over a singleton set thereby simplifying (\ref{eqn:HJI_appx_con}) to
\begin{equation}
    H_{const}(s, p) = \min_{d \in \mathcal{D}} p^T [ \hat{f}(s, \pi(s)) + d ]
\end{equation}
where the modelled dynamics from (\ref{eqn:aug_dynamics}) are substituted above.
It is observed that the vector product $p^T \hat{f}(s, \pi(s))$ is independent of the design variable $d$ and can be taken out of the optimization.
\begin{equation}
    H_{const}(s, p) = p^T \hat{f}(s, \pi(s)) \; + \; \min_{d \in \mathcal{D}} p^T d
    \label{eqn:H_opt1}
\end{equation}
Further, it is possible to derive an analytical solution of (\ref{eqn:H_opt1}).
By expanding the objective function $\min_{d \in \mathcal{D}} p^T d = \min_{d \in \mathcal{D}} (p_1 d_1 + p_2 d_2 + \hdots + p_n d_n)$ it is evident that depending on the sign of the costate terms, an optimal minimizing disturbance input can be computed as
\begin{equation}
d_i = \begin{cases} -d_i
    &\mbox{if } p_i > 0 \\
      d_i
    &\mbox{if } p_i < 0 \end{cases} , \forall i=1,\hdots n,
\label{eqn:H_opt2}
\end{equation}
The analytical solution is therefore given as
\begin{equation}
\begin{split}
    H_{const}(s, p) = p^T \hat{f}(s, &\pi(s))  \\
    - &(|p_1| d_1 + |p_2| d_2 + \hdots + |p_n| d_n).
    \end{split}
    \label{eqn:H_opt2}
\end{equation}

We are now ready to present the modified BRT problem formulation as
\begin{equation}
\begin{split}
    D_t \phi(s,t) + \min[0, H_{const}(s, \; D_x \phi(s,t))] &= 0, \\
    \phi(s, 0) &= l(s) \\
\text{where, } H_{const}(s, p) = p^T \hat{f}(s, \pi(s)&) \\
- (|p_1| d_1 + |p_2| d_2 +& \hdots + |p_n| d_n)
\end{split}
\label{eqn:PDE_HJI_v2}
\end{equation}

  \begin{theorem}
  Let $\phi$ be the viscosity solution of the terminal valued PDE in (\ref{eqn:PDE_HJI_v2}), then the zero sublevel set of $\phi$ describes the backward reachable set $B_C(t)$ as
  \begin{equation}
      B_C(t) = \{ s \in \mathbb{R}^n | \phi(s, t) \le 0\}.
      \label{eqn:thm_main}
  \end{equation}
  \label{thm:main}
  \end{theorem}
  %
  \begin{proof}
  For the proof of Theorem \ref{thm:main}, we leverage the result from Theorem \ref{thm1}.
  %
  To do so, we first ensure that our formulation satisfies the necessary requirements for Theorem \ref{thm1} to hold:
  %
 \paragraph{For Requirement 1 } From assumptions \ref{ass1} and \ref{ass2} it has been shown that an ANN modelled dynamics satisfies the properties of boundedness, uniform continuity and Lipschitz continuity for a large choice of activation functions. Therefore, the learnt dynamics $\hat{f}(s, a, d)$ satisfies requirement 1.
 %
 \paragraph{For Requirement 2 } From the discussion in Section \ref{sec:prelim}-C, a bounded and Lipschitz continuous function $l$ can easily be derived for a bounded target set $\mathcal{L}$ by representing it as a union of several well-defined geometric shapes. Therefore, for bounded constrained regions $\mathcal{C}_s$ this requirement is also satisfied.
 \paragraph{For Requirement 3 } Finally, from assumption \ref{ass3} it is established that the input sets $\mathcal{A}, \mathcal{D} \subset \mathbb{R}^m$ are compact sets.

  Since our formulation satisfies all the necessary requirements of Theorem 1, the result of Theorem \ref{thm:main} follows directly from it.
  %
  Therefore, by computing the constrained Hamiltonian $H_{const}$ and solving for the viscosity solution $\phi$ in (\ref{eqn:PDE_HJI_v2}), the backward reachable set at any time $t$ is implicitly defined using (\ref{eqn:thm_main}).
  \end{proof}
  \paragraph{Step 3: \textit{Determine $\mathcal{S}_{safe}$}}

  In the final step, $\mathcal{S}_{safe}$ is determined from the solution of problem \ref{prob2}.
  %
  From the definitions of $\mathcal{S}_{safe}$ and $\mathcal{S}_{unsafe}$ in problem \ref{prob1} and \ref{prob2} respectively the following two statements are evident.
  %
  \begin{equation}
  \begin{split}
    \mathcal{S}_{safe} \cup \mathcal{S}_{unsafe} &= \mathcal{S}_0\\
    \mathcal{S}_{safe} \cap \mathcal{S}_{unsafe} &= \emptyset
    \end{split}
    \label{eqn:set_int}
  \end{equation}
  Therefore, it follows
  \begin{equation}
    \mathcal{S}_{safe} = \overline{\mathcal{S}}_{unsafe} \cap \mathcal{S}_0.
  \end{equation}

After determining $\mathcal{S}_{safe}$ the policy can be classified as either completely safe, completely unsafe or partially safe by checking the conditions in (\ref{eqn:safe_policy}).

\begin{comment}
\subsection{Algorithm Overview}
The proposed algorithm solves problem 1 in two steps.
In step 1 a complementary problem is solved to determine $\mathcal{S}_{unsafe}$ and in step 2 the solution from  problem 2 is used to determine $\mathcal{S}_{safe}$.

\textit{Step 1:}
Instead of directly computing $\mathcal{S}_{safe}$, a complementary problem is framed as follows:
\textbf{Problem 2:}
For given unsafe set $\mathcal{C}_s$ determine the unsafe subset of initial states $\mathcal{S}_{unsafe} \subset \mathcal{S}_0$ such that for any initial state $s_0 \in \mathcal{S}_{unsafe}$ there exists at least one trajectory $\tau$ under deterministic policy $\pi(s)$ which will drive the system state to $\mathcal{C}_s$ within a finite time.

To determine $\mathcal{S}_{unsafe}$ in problem 2, the backward reachable tube $B_R(T)$ from $\mathcal{C}_s$ is computed using (\ref{eqn:PDE_HJI}).
Then $\mathcal{S}_{unsafe} = B_R(T) \cap \mathcal{S}_0$.
The mathematical formulation is given in the following subsection.

\textit{Step 2:}
In this step, $\mathcal{S}_{safe}$ is determined from the solution of problem 2.
From the definitions of $\mathcal{S}_{safe}$ and $\mathcal{S}_{unsafe}$ in problem 1 and 2 respectively the following two statements are evident.
  \begin{equation}
  \begin{split}
    \mathcal{S}_{safe} \cup \mathcal{S}_{unsafe} &= \mathcal{S}_0\\
    \mathcal{S}_{safe} \cap \mathcal{S}_{unsafe} &= \emptyset
    \end{split}
    \label{eqn:set_int}
  \end{equation}
Therefore, it follows
\begin{equation}
    \mathcal{S}_{safe} = \overline{\mathcal{S}}_{unsafe} \cap \mathcal{S}_0.
\end{equation}

\subsection{Mathematical Formulation}
Since the actual system dynamics is unknown and only a learnt model $\hat{f}(s, a)$ is available, we first determine the bounds on the modelling error $d \in \mathcal{D} \subset \mathbb{R}^n$.
The bounds are determined for each state in the $n$-dimensional input vector such that $d_{min} < d < d_{max}$.
From the formulation in (\ref{eqn:model_error}) the two player system dynamics can be represented as
\begin{equation}
    \hat{f}(s, a, d) = \hat{f}(s, a) + d.
    \label{eqn:dyn_ad}
\end{equation}
The HJI PDE formulation remains the same however the hamiltonian is modified as follows:
\begin{equation}
\begin{split}
    H_{const}(s, p) &= \max_{a \in \mathcal{A}} \min_{d \in \mathcal{D}} p^T \hat{f}(s, a, d) \\
    \text{s.t. } &a = \pi(s)
\end{split}
\label{eqn:HJI_appx_con}
\end{equation}
Since the action $a$ is computed by the deterministic policy $\pi(s)$, the outer maximization optimization is carried over a singleton set thereby simplifying (\ref{eqn:HJI_appx_con}) to
\begin{equation}
    H_{const}(s, p) = \min_{d \in \mathcal{D}} p^T [ \hat{f}(s, \pi(s)) + d ]
\end{equation}
where the modelled dynamics from (\ref{eqn:dyn_ad}) are substituted above.
It is observed that the vector product $p^T \hat{f}(s, \pi(s))$ is independent of the design variable $d$ and can be taken out of the optimization.
\begin{equation}
    H_{const}(s, p) = p^T \hat{f}(s, \pi(s)) \; \min_{d \in \mathcal{D}} p^T d
    \label{eqn:H_opt1}
\end{equation}
Further, it is possible to derive an analytical solution of (\ref{eqn:H_opt1}) by expanding the objective function $p^T d$ as
\begin{equation}
    H_{const}(s, p) = p^T \hat{f}(s, \pi(s)) \; \min_{d \in \mathcal{D}} (p_1 d_1 + p_2 d_2 + \hdots + p_n d_n).
    \label{eqn:H_opt2}
\end{equation}
The analytical solution is therefore given as
\begin{equation}
    H_{const}(s, p) = p^T \hat{f}(s, \pi(s))  - (|p_1| d_1 + |p_2| d_2 + \hdots + |p_n| d_n).
    \label{eqn:H_opt3}
\end{equation}

We are now ready to present the modified backward reachable tube problem formulation as
\begin{equation}
\begin{split}
    \frac{\partial J}{\partial t} + \min[0, H_{const}(s, \; \frac{\partial J}{\partial s})] &= 0, \\
    J(s, 0) &= l(s) \\
\text{where,} \; H_{const}(s, p) = p^T \hat{f}(s, \pi(s)) &- (|p_1| d_1 + |p_2| d_2 + \hdots + |p_n| d_n)
\end{split}
\label{eqn:PDE_HJI_v2}
\end{equation}

\end{comment}

\section{Problem Setting}
Let ${\mathbb{S}} \subset \mathbb{R}^n$ denote the set of states and ${\mathbb{A}} \subset \mathbb{R}^m$ denote the set of feasible actions for the RL agent.
Let ${\mathbb{S}}_0 \subset {\mathbb{S}}$ denote the set of bounded initial states and $\xi := \{ ({\bm{s}}_t, {\bm{a}}_t) \}_{t=0}^T$ represent a trajectory generated over a finite time $T$, as a sequence of state and action tuples, where subscript $t$ denotes the instantaneous time.
Additionaly, let ${\bm{s}}(\cdot)$ and ${\bm{a}}(\cdot)$ represent a sequence of states and actions, respectively.
The state constraints imposed on the system are represented as unsafe regions using bounded sets ${\mathbb{C}}_s = \cup_{i=1}^p {\mathbb{C}}_s^{(i)}$, where ${\mathbb{C}}_s^{(i)} \subset {\mathbb{S}}, \; \forall i \in \{1, 2, \hdots p \}$.
The true system dynamics is given by a non-linear function $f: {\mathbb{S}} \times {\mathbb{A}} \rightarrow {\mathbb{R}}^n$ such that, $\dot{{\bm{s}}} = f({\bm{s}}, {\bm{a}})$, and is unknown to the agent.

A model-based RL algorithm is used to find an optimal controller $\pi: {\mathbb{S}} \rightarrow {\mathbb{A}}$, to reach a set of target states ${\mathbb{T}} \subset {\mathbb{S}}$ within some finite time $T$, while avoiding constraints ${\mathbb{C}}_s$.
An NN model, $\hat{f}_{\bm{\theta}}: {\mathbb{S}} \times {\mathbb{A}} \rightarrow {\mathbb{R}}^n$ parameterized by weights ${\bm{\theta}}$, is trained to learn the true, but unknown, system dynamics from the observed state transition data tuples $D = \{ ( {\bm{s}}_t, {\bm{a}}_t, \Delta {\bm{s}}_{t+1})^{(i)} \}_{i=1}^{N}$.
However, due to sampling bias, the learnt model $\hat{f}_{\bm{\theta}}$ may not be accurate.
We assume that it is possible to estimate a bounded set ${\mathbb{D}} \subset {\mathbb{R}}^n$ such that, at any state ${\bm{s}} \in {\mathbb{S}}$, augmenting the learnt dynamics $\hat{f}_{\bm{\theta}}$ with some ${\bm{d}} \in {\mathbb{D}}$ results in a closer approximation of the true system dynamics at that particular state.
Using this notation, we now define the problem of safety verification of a given controller $\pi({\bm{s}})$.
\begin{problem} \label{prob1_safety_verification}
(Safety verification): Given a set of initial states ${\mathbb{S}}_0$, determine if  $\,\, \forall {\bm{s}}_0 \in {\mathbb{S}}_0$, all the trajectories $\xi$ executed under $\pi({\bm{s}})$ and following the system dynamics $f$, satisfy the constraints ${\mathbb{C}}_s$ or not.
\end{problem}
The solution to Problem \ref{prob1_safety_verification} will only provide a binary yes or no answer to whether $\pi({\bm{s}})$ is safe or not with respect to ${\mathbb{S}}_0$.
In the case where the policy is unsafe, a stronger result is the identification of safe initial states ${\mathbb{S}}_{safe} \subset {\mathbb{S}}_0$ from which $\pi({\bm{s}})$ executes trajectories which always satisfy the constraints ${\mathbb{C}}_s$.
This problem is stated below.
\begin{problem} \label{prob2_safe_states}
(Safe initial states): Given $\pi({\bm{s}})$, find ${\mathbb{S}}_{safe}$, such that, any trajectory $\xi$ executed under $\pi({\bm{s}})$ and following the  system dynamics $f$, starting from any ${\bm{s}}_0 \in {\mathbb{S}}_{safe}$, satisfies the constraints ${\mathbb{C}}_s$.
\end{problem}

\section{Safety Verification}
To address Problems \ref{prob1_safety_verification} and \ref{prob2_safe_states}, we use  \textit{reachability analysis}.
Reachability analysis is an exhaustive verification technique which tracks the evolution of the system states over a finite time, known as the \textit{reachable tube}, from a given set of initial states ${\mathbb{S}}_0$ \citep{maler2008computing}.
If the evolution is tracked starting from ${\mathbb{S}}_0$, then it is called the \textit{forward reachable tube} and is denoted as ${\mathbb{F}}_R(T)$.
Analogously, if the evolution is tracked starting from ${\mathbb{C}}_s$ to ${\mathbb{S}}_0$, then it is called the \textit{backward reachable tube} and is denoted as ${\mathbb{B}}_R(T)$.

In the following sections, we will formulate a reachable tube problem for NN-based models and controllers, and then propose a verification framework that (a) can determine whether or not a given policy $\pi$ is safe, and (b) can compute ${\mathbb{S}}_{safe}$ if $\pi$ is unsafe.
To do so, there are two main questions that need to be answered.
First, since the true system dynamics $f$ is unknown, how can we determine a conservative bound on the modeling error, to augment the learnt model $\hat{f}_{\bm{\theta}}$ and better model the true system dynamics when evaluating a controller $\pi$?
Second, how do we formulate the \textit{forward} and \textit{backward reachable tube} problems over the NN modeled system dynamics?

\subsection{Model-Based Reinforcement Learning}
In this section, we focus on the necessary requirements for the modeled dynamics $\hat f_{\bm{\theta}}$ and discuss the estimation of the modeling error set ${\mathbb{D}}$.
A summary of the model-based RL framework is presented in Appendix \ref{sec:modelRL}.
Recall that the learnt model $\hat{f}_{\bm{\theta}}$ is represented using an NN

and predicts the change in states $\Delta \hat{{\bm{s}}}_{t+1} \in {\mathbb{R}}^n$.
To learn $\hat f_\theta$, an observed data set $D= \{ ( {\bm{s}}_t, {\bm{a}}_t, \Delta {\bm{s}}_{t+1})^{(i)} \}_{i=1}^{N}$ is first split into
a training data set $D_t$ and a validation data set $D_v$.
A supervised learning technique is then used to train $\hat f_\theta$ over $D_t$ by minimizing the prediction error
  $ E = \frac{1}{N_t} \sum_{i=1}^{N_t} || \Delta \hat{{\bm{s}}}_{t+1}^{(i)} - \Delta {\bm{s}}_{t+1}^{(i)} ||^2,$
where $\Delta \hat{{\bm{s}}}_{t+1}$ is the change predicted by the learnt model, $\Delta {\bm{s}}_{t+1}$ is the true observed change in state and $N_t = |D_t|$.
With this notation, we now formalize the following necessary assumption for this work.
This assumption is required to ensure boundedness during analysis and is easily satisfied by NNs that use common activation functions like $\texttt{tanh}$ or $\texttt{sigmoid}$ \citep{usama2018towards}.

\begin{assumption} $\hat{f}_{{\bm{\theta}}}$ is Lipschitz continuous and $\forall j \in \{1, ..., n\}, \Delta \hat{{\bm{s}}}_j \in
[-c, c]$, where $|c| < \infty$.
  \label{ass:boundedness}
\end{assumption}

\textbf{Modeling error:}
As mentioned earlier, the accuracy of the learnt model $\hat{f}_{{\bm{\theta}}}$ depends on the quality of data and the NN being used, thereby resulting in some modeling error ${\bm{d}}$ in the prediction of the next state.
Estimating modeling errors is an active area of research and is required for several existing works on safe RL \citep{akametalu2014reachability, gillula2012guaranteed}, and is complementary to our goal.
Since the primary contribution of this work is the development of a reachable tube formulation for model-based controllers that use NNs, we rely on existing techniques \citep{moldovan2015optimism} to estimate a conservative modeling error bound.
We leverage the error estimates $\hat {\bm{d}}_j = \hat \Delta s_{t+1}^{(j)} - \Delta s^{(j)}_{t+1},$ of $\hat f_\theta$, for the transition tuples in the validation set $D_v$ to construct the upper confidence bound ${\bm{d}}^+ = [d_1, \; d_2, \hdots d_n]$ and lower confidence bound ${\bm{d}}^- = [-d_1, -d_2, \hdots -d_n]$  for ${\bm{d}}$, for each state dimension.
Let a high-confidence bounded error set ${\mathbb{D}}$ be defined such that ${\mathbb{D}} = \{{\bm{d}}: \forall i \in \{1,...,n\}, \,\, {\bm{d}}^-_i < {\bm{d}}_i < {{\bm{d}}}^+_i\}$.
We then use ${\mathbb{D}}$ to represent the augmented learnt system dynamics as
\begin{equation}
    \hat{f}^{(r)}_{{\bm{\theta}}} := \hat{f}_{{\bm{\theta}}}({\bm{s}}, {\bm{a}}) + {\bm{d}}, \; \;  {\bm{d}} \in {\mathbb{D}}.

    \label{eqn:aug_r_dynamics}
\end{equation}

\subsection{Reachable Tube Formulation}

For an exhaustive verification technique on a continuous state and action space, it is infeasible to sample trajectories from \textit{every} point in the given initial state and further, to verify whether \textit{all} these trajectories satisfy the safety constraints.
Therefore, reachable sets are usually (approximately) represented as convex polyhedrons and their evolution is tracked by pushing the boundaries of this polyhedron according to the system dynamics.
However, as convex polyhedrons can lead to large approximation errors, we leverage the level set method \citep{mitchell2005time} to compute the boundaries of the reachable tube for NN-based models and controllers at every time instant.
With the level set method, even non-convex boundaries of the reachable tube can be represented, thereby ensuring an accurate computation of the reachable tube \citep{mitchell2005time}.
Since the non-linear, non-convex structure of NNs is not suitable for analytical verification techniques, the reachability analysis gives an efficient, simulation-based approach to analyze the safety of the controller.
Therefore, in this subsection, we formulate the reachable tube for a NN modeled system dynamics $\hat{f}_{{\bm{\theta}}}$, but first we formally define the forward reachable tube and backward reachable tube for a policy $\pi({\bm{s}})$.

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.8\textwidth]{images/Flowchart_MBRL_tiny_black_dots.png}
\end{center}
\caption{Flowchart of the proposed safety verification algorithm. The rectangle on the left represents the flow for computing the forward reachable tube (FRT), which can only state if $\pi({\bm{s}})$ is safe or unsafe for ${\mathbb{S}}_0$. The rectangle on the right presents the flow for computing the backward reachable tube (BRT), which is invoked if $\pi({\bm{s}})$ is unsafe. The BRT analysis can compute the subset of safe initial states ${\mathbb{S}}_{safe}$ for $\pi({\bm{s}})$, if such a set exists.}
\label{fig:flowchart}
\end{figure}

\textbf{Forward reachable tube (FRT):}
It is the set of all states that can be \textit{reached from} an initial set ${\mathbb{S}}_0$, when the trajectories $\xi$ are executed under policy $\pi({\bm{s}})$ and system dynamics $\hat{f}^{(r)}_{{\bm{\theta}}}({\bm{s}}, {\bm{a}}, {\bm{d}})$.
The FRT is computed over a finite length of time $T$ and is formally defined as
\begin{equation}
    \!\!\!{\mathbb{F}}_R(T) := \{ {\bm{s}} : \forall {\bm{d}} \in {\mathbb{D}}, {\bm{s}}(\cdot) \text{ satisfies } \dot{{\bm{s}}} = \hat{f}^{(r)}_{{\bm{\theta}}}({\bm{s}}, {\bm{a}}, {\bm{d}}), \text{ where } {\bm{a}} = \pi({\bm{s}}), \; {\bm{s}}_{t_0} \in {\mathbb{S}}_0, t_f = T \},
    \label{eqn:FRT_def}
\end{equation}
where $t_0$ and $t_f$ denote the initial and final time of the trajectory $\xi$, respectively.

\textbf{Backward reachable tube (BRT):}
It is the set of all states which \textit{can reach} a given bounded target set ${\mathbb{T}} \subset {\mathbb{R}}^n$, when the trajectories $\xi$ are executed under policy $\pi({\bm{s}})$ and system dynamics $\hat{f}^{(r)}_{{\bm{\theta}}}({\bm{s}}, {\bm{a}}, {\bm{d}})$.
The BRT is also computed for a finite length of time, with the trajectories starting at time $t_0 = -T$ and ending at time $t_f = 0$.
It is denoted as
\begin{equation}
\begin{split}
{\mathbb{B}}_R(-T) := \{ {\bm{s}}_0 : \forall {\bm{d}} \in {\mathbb{D}}, {\bm{s}}(\cdot) &\text{ satisfies } \dot{{\bm{s}}} = \hat{f}^{(r)}_{{\bm{\theta}}} ({\bm{s}}, {\bm{a}}, {\bm{d}}), \text{ where } {\bm{a}} = \pi({\bm{s}})\\ &\text{ with } {\bm{s}}_{t_0} = {\bm{s}}_{-T}; \; {\bm{s}}_{t_f} \in {\mathbb{T}}, t_f \in [-T, 0] \}.
\end{split}
\end{equation}
The key difference between the FRT and BRT is that, for the former, the initial set of states are known, whereas for the latter, the final set of states are known.

\paragraph{Outline:}
The flowchart of the safety verification framework proposed in this work is presented in Fig. \ref{fig:flowchart}.
Given a model-based policy $\pi({\bm{s}})$, the set of initial states ${\mathbb{S}}_0$ and the set of constrained states ${\mathbb{C}}_s$, the first step is to estimate the bounded set of modeling error ${\mathbb{D}}$, as discussed in Section 3.1.
Using $\hat{f}^{(r)}_{\bm{\theta}}$, the FRT is constructed from the initial set ${\mathbb{S}}_0$ and it contains all the states reachable by $\pi({\bm{s}})$ over a finite time $T$.
Thus, if the FRT contains any state from the unsafe region ${\mathbb{C}}_s$, $\pi({\bm{s}})$ is deemed \textit{unsafe}.
Therefore, the solution to Problem \ref{prob1_safety_verification}, is determined by analyzing the set of intersection of the FRT with ${\mathbb{C}}_s$ as
\begin{equation}
    \pi = \begin{cases} \text{ safe}
    &\mbox{if } \,\, {\mathbb{F}}_R(T) \cap {\mathbb{C}}_s = \emptyset, \\
    \text{ unsafe} &\mbox{if }\,\, {\mathbb{F}}_R(T) \cap {\mathbb{C}}_s \neq \emptyset.  \end{cases}
\end{equation}
If $\pi({\bm{s}})$ is classified as \textit{safe} for the entire set ${\mathbb{S}}_0$, then no further analysis is required.
However, if $\pi({\bm{s}})$ is classified as \textit{unsafe}, we proceed to compute the subset  ${\mathbb{S}}_\text{safe} \subset {\mathbb{S}}_0$ of initial states for which $\pi({\bm{s}})$ generates safe trajectories.
${\mathbb{S}}_{safe}$ is the solution to Problem \ref{prob2_safe_states} and allows an unsafe policy to be deployed on a real system, with restrictions on the starting states.
To this end, the BRT is computed from the unsafe region ${\mathbb{C}}_s$, to determine the set of trajectories (and states) which terminate in ${\mathbb{C}}_s$.
The intersection of the BRT with ${\mathbb{S}}_0$ determines the set of unsafe initial states ${\mathbb{S}}_{unsafe}$.
To determine ${\mathbb{S}}_{safe}$, we utilize the following properties, (a) ${\mathbb{S}}_{safe} \cup {\mathbb{S}}_{unsafe} = {\mathbb{S}}_0$, and (b) ${\mathbb{S}}_{safe} \cap {\mathbb{S}}_{unsafe} = \emptyset$,
and compute
 \begin{equation}
    {\mathbb{S}}_{safe} = \overline{{\mathbb{S}}}_{unsafe} \cap {\mathbb{S}}_0.
\end{equation}
If ${\mathbb{S}}_{safe} \neq \emptyset$, then we have identified the safe initial states for $\pi({\bm{s}})$, otherwise, it is concluded that there are no initial states in ${\mathbb{S}}_0$ from which $\pi({\bm{s}})$ can generate safe trajectories.
\paragraph{Mathematical Formulation:}
This section presents the mathematical formulation to compute the BRT.
The FRT can be computed with a slight modification to the BRT formulation and this is discussed in the end of this section.

Recall, for the BRT problem, there exists a target set ${\mathbb{T}} \subset {\mathbb{R}}^n$ which the agent has to reach in finite time, i.e., the condition on the final state is given as ${\bm{s}}_{t_f} \in {\mathbb{T}}$.
Conventionally, for the BRT formulation, the final time $t_f = 0$ and the starting time $t_0 = -T$, where $0 < T < \infty$.
When evaluating a policy $\pi({\bm{s}})$, the controller input is computed by the given policy as ${\bm{a}} = \pi({\bm{s}})$.
However, following the system dynamics $\hat{f}_{\bm{\theta}}^{(r)}$ in (\ref{eqn:aug_r_dynamics}), the modeling error ${\bm{d}}$ is now included in the system as an adversarial input, whose value at each state is determined so as to maximize the controller's cost function.
We use the HJ PDE to formulate the effect of the modeling error on the system for computing the BRT, but first we briefly review the formulation of the HJ PDE with an NN modeled system dynamics $\hat{f}_{\bm{\theta}}$ in the following.

For an optimal controller, we first define the cost function which the controller has to minimize.
Let $C({\bm{s}}_t, {\bm{a}}_t)$ denote the running cost of the agent, which is dependent on the state and action taken at time $t \in [-T, 0]$.
Let $g({\bm{s}}_{t_f})$ denote the cost at the final state ${\bm{s}}_{t_f}$.
Then, the goal of the optimal controller is to find a series of optimal actions such that
\begin{equation}
\begin{split}

    \min_{{\bm{a}}_{\tau}(\cdot)} &\left(\int_{-T}^{0} C({\bm{s}}_{\tau}, {\bm{a}}_{\tau}) d \tau + g({\bm{s}}_{t_f}) \right)\\
    \text{subject to  } \dot{{\bm{s}}} &= \hat{f}_{{\bm{\theta}}}({\bm{s}}, {\bm{a}}), \,\,\,

    {\bm{s}}_{t_f} \in {\mathbb{T}},
\end{split}
\label{eqn:J_v1}
\end{equation}
where ${\bm{a}}_{\tau} \in {\mathbb{A}}$ and $\hat{f}_{\bm{\theta}}$ is the NN modeled system dynamics.
The above optimization problem is solved using the \textit{dynamic programming} approach \citep{smith1991dynamic}, which is based on the \textit{Principle of Optimality} \citep{troutman2012variational}.
Let $V({\bm{s}}_t, t)$ denote the value function of a state ${\bm{s}}$ at time $t \in [-T, 0]$, such that
\begin{equation}
    \!\!\! V({\bm{s}}_t, t) =
    \min_{{\bm{a}}_\tau(\cdot)} \left[ \int_{t}^{0} C({\bm{s}}_{\tau}, {\bm{a}}_{\tau}) d \tau + g({\bm{s}}_{t_f}) \right]  = \min_{{\bm{a}}_\tau(\cdot)} \left[ \int_{t}^{t + \delta} \!\!\! C({\bm{s}}_{\tau}, {\bm{a}}_{\tau}) d \tau + V({\bm{s}}_{t+\delta}, t + \delta ) \right],
    \label{eqn:V_def}
\end{equation}
where $\delta > 0$.
$V({\bm{s}}_t, t)$ is a quantitative measure of being at a state ${\bm{s}}$, described in terms of the cost required to reach the goal state from ${\bm{s}}$.

Then, using the Taylor series expansion, $V({\bm{s}}_{t + \delta}, t + \delta)$ is approximated around $V({\bm{s}}_t, t)$ in (\ref{eqn:V_def}) to derive the HJ PDE as
\begin{equation}
\begin{split}
\frac{dV}{dt} + \min_{{\bm{a}}} \left[ \nabla V \cdot \hat{f}_{{\bm{\theta}}}({\bm{s}}, {\bm{a}}) + C({\bm{s}}, {\bm{a}}) \right] &= 0, \\
V({\bm{s}}_{t_f}, t_f) &= g({\bm{s}}_{t_f}),
\end{split}
\label{eqn:HJI_PDE_value}
\end{equation}
where $\nabla V \in {\mathbb{R}}^n$ is the spatial derivative of $V$.
Additionally, the time index has been dropped above and the dynamics constraint in (\ref{eqn:J_v1}) has been included in the PDE.
Equation (\ref{eqn:HJI_PDE_value}) is a terminal value PDE, and by solving (\ref{eqn:HJI_PDE_value}), we can compute the value of a state $V({\bm{s}}_t, t)$ at any time $t$.

We now discuss how the formulation in (\ref{eqn:HJI_PDE_value}) can be modified to obtain the BRT.
It is noted that along with computing the value function, the formulation in (\ref{eqn:HJI_PDE_value}) also computes the optimal action ${\bm{a}}$.
However, in Problems \ref{prob1_safety_verification} and \ref{prob2_safe_states}, the optimal policy $\pi({\bm{s}})$ is already provided.
Therefore, the constraint ${\bm{a}} = \pi({\bm{s}})$ should be included in problem (\ref{eqn:J_v1}), thereby avoiding the need of minimizing over actions ${\bm{a}} \in {\mathbb{A}}$ in (\ref{eqn:HJI_PDE_value}).
Additionally, as discussed in Section 3.1, the NN modeled system dynamics $\hat{f}_{\bm{\theta}}$ may not be a good approximation of the true system dynamics $f$.
Instead, the augmented learnt system dynamics $\hat{f}^{(r)}_{\bm{\theta}}$ in (\ref{eqn:aug_r_dynamics}) is used in place of $\hat{f}_{\bm{\theta}}$ in (\ref{eqn:HJI_PDE_value}), since it better models the true dynamics at a given state.
However, by including $\hat{f}_{\bm{\theta}}^{(r)}$ in (\ref{eqn:HJI_PDE_value}), the modeling error ${\bm{d}}$ is now included in the formulation.
The modeling error ${\bm{d}} \in {\mathbb{D}}$ is treated as an adversarial input which is trying to drive the system away from it's goal state by taking a value which maximizes the cost function at each state.
Thus, to account for this adversarial input, the formulation in (\ref{eqn:HJI_PDE_value}) is now maximized over ${\bm{d}}$.

Lastly, the BRT problem is posed for a set of states and not an individual state.
Hence, an efficient representation of the target set is required to propagate an entire set of trajectories at a time, as opposed to propagating individual trajectories.
\begin{assumption} \label{ass:level_set}
The target set ${\mathbb{T}} \subset {\mathbb{R}}^n$ is closed and can be represented as the zero sublevel set of a bounded and Lipschitz continuous function $l: {\mathbb{R}}^n \
\rightarrow {\mathbb{R}}$, such that, ${\mathbb{T}} = \{ {\bm{s}} : l({\bm{s}}) \le 0\}$.
\end{assumption}
The above assumption defines a function $l$ to check whether a state lies inside or outside the target set.
If ${\mathbb{T}}$ is represented using a regular, well-defined geometric shape (like a sphere, rectangle, cylinder, etc.), then deriving the level set function $l({\bm{s}})$ is straight forward, whereas, an irregularly shaped ${\mathbb{T}}$ can be represented as a union of several well-defined geometric shapes to derive $l({\bm{s}})$.

For the BRT problem, the goal is to determine all the states which can reach ${\mathbb{T}}$ within a finite time.
The path taken by the controller is irrelevant and only the value of the final state is used to determine if any state ${\bm{s}} \in {\mathbb{B}}_R(-T)$.
From Assumption \ref{ass:level_set}, the terminal condition ${\bm{s}}_{t_f} \in {\mathbb{T}}$ can be restated as $l({\bm{s}}_{t_f}) \le 0$.
Thus, to prevent the system from reaching ${\mathbb{T}}$, the adversarial input ${\bm{d}}$ tries to maximize $l({\bm{s}}_{t_f})$, thereby pushing ${\bm{s}}_{t_f}$ as far away from ${\mathbb{T}}$ as possible.
Therefore, the cost function in (\ref{eqn:J_v1}) is modified to $J = l({\bm{s}}_{t_f})$.
Additionally, any state which can reach ${\mathbb{T}}$ \textit{within} a finite time interval $T$ is included in the BRT.
Therefore, if any trajectory reaches ${\mathbb{T}}$ at some $t_f<0$, it shouldn't be allowed to leave the set.
Keeping this in mind, the BRT optimization problem can be posed as
\begin{equation}
\begin{split}
\max_{{\bm{d}}(\cdot)} & \left( \min_{t \in [-T, 0]} l({\bm{s}}_{t_f}) \right)\\
\text{subject to: } \dot{{\bm{s}}} &= \hat{f}^{(r)}_{{\bm{\theta}}} ({\bm{s}}, {\bm{a}}, {\bm{d}}), \,\,\,
{\bm{a}} = \pi({\bm{s}}), \,\,\,\,
l({\bm{s}}_{t_f}) \le 0,
\end{split}
\label{eqn:J_level}
\end{equation}
where the inner minimization over time prevents the trajectory from leaving the target set.
Then, the value function for the above problem is defined as
\begin{equation}
    V_{R}({\bm{s}}_t, t) := \max_{{\bm{d}}(\cdot)} l({\bm{s}}(t_f)).
\label{eqn:v_r}
\end{equation}
Comparing this with (\ref{eqn:V_def}), it is observed that the value of a state ${\bm{s}}$ is no longer dependent on the running cost $C({\bm{s}}, {\bm{a}})$.
This doesn't imply that the generated trajectories are not optimal w.r.t. action ${\bm{a}}$, because the running cost is equivalent to the negative reward function, for which $\pi({\bm{s}})$ is already optimized.
Instead, $V_R$ solely depends on whether the final state of the trajectory lies within the target set or not, i.e., whether or not ${\bm{s}}_{t_f} \in {\mathbb{T}}$.
Thus, the value function $V_R$ for any state ${\bm{s}}$ is equal to $l({\bm{s}}_{t_f})$, where ${\bm{s}}_{t_f}$ is the final state of the trajectory originating at ${\bm{s}}$.
Then, the HJ PDE for the problem in (\ref{eqn:J_level}) is stated as
\begin{equation}
\begin{split}
\frac{dV_R}{dt} + min \{ 0, H^*({\bm{s}}, \nabla V_R({\bm{s}}_t, t), t) \} &= 0, \\
V_R({\bm{s}}_{t_f}, t_f) &= l({\bm{s}}_{t_f}), \\
\text{where } H^* &= \max_{{\bm{d}}} \left( \nabla V_R \cdot \hat{f}^{(r)}_{{\bm{\theta}}}({\bm{s}}, \pi({\bm{s}}), {\bm{d}}) \right), \\
\end{split}
\label{eqn:BRT_v1}
\end{equation}
where $H^*$ represents the optimal Hamiltonian.
Since we are computing the BRT, $min\{ 0, H^*({\bm{s}}, \nabla V_R({\bm{s}}_t, t), t) \}$ in the PDE above ensures that the tube grows only in the backward direction, thereby preventing a trajectory which has reached ${\mathbb{T}}$ from leaving.
In Problem (\ref{eqn:BRT_v1}), the optimal action has been substituted by $\pi({\bm{s}})$ and the augmented learnt dynamics is used instead of $\hat{f}_{\bm{\theta}}$.
The Hamiltonian optimization problem can be further simplified to derive an analytical solution for the modeling error.
By substituting the augmented dynamics from (\ref{eqn:aug_r_dynamics}), the optimization problem can be re-written as
\begin{equation}
    H^* = \nabla V_R \cdot \hat{f}_{{\bm{\theta}}}({\bm{s}}, {\bm{a}}) + \max_{{\bm{d}}} \nabla V_R \cdot {\bm{d}}.
\end{equation}
Expanding $\nabla V_R = [p_1, p_2, \hdots, p_n]^T \in {\mathbb{R}}^n$, the vector product $\nabla V_R \cdot {\bm{d}} = p_1 d_1 + p_2 d_2 + \hdots + p_n d_n$.
Therefore, to maximize $\nabla V_R \cdot {\bm{d}}$, the disturbance control is chosen as
\begin{equation}
d_i = \begin{cases} d_i
    &\mbox{if } p_i > 0 \\
      -d_i
    &\mbox{if } p_i < 0 \end{cases} , \forall i=1,\hdots n.
\label{eqn:H_opt2}
\end{equation}
With this analytical solution, the final PDE representing the BRT is stated as
\begin{equation}
\begin{split}
\frac{dV_R}{dt} + min \{ 0, H^*({\bm{s}}, \nabla V_R({\bm{s}}_{t}, t), t) \} &= 0, \\
V_R({\bm{s}}_{t_f}, t_f) &= l({\bm{s}}_{t_f}), \\
\text{where } H^* = \nabla V_R \cdot \hat{f}_{{\bm{\theta}}}({\bm{s}}, {\bm{a}}) + &(|p_1| d_1 + |p_2| d_2 + \hdots + |p_n| d_n).
\end{split}
\label{eqn:BRT_final}
\end{equation}

The value function $V_R({\bm{s}}_t, t)$ in (\ref{eqn:BRT_final}) represents the evolution of the target level set function backwards in time.
By finding the solution to $V_R$ in the above PDE, the level set function is determined at any time instant $t \in [-T, 0]$, thereby determining the BRT.
From the result of Theorem 2 in \citet{mitchell2005time}, it is proved that the solution of $V_R$ in (\ref{eqn:BRT_final}) at any time $t$ gives the zero sublevel set for the BRT.
Thus,
\begin{equation}
    {\mathbb{B}}_R(-T) = \{ {\bm{s}} : V_R({\bm{s}}_t, t) \le 0, \; t \in [-T, 0] \}.
\end{equation}
The solution to $V_R$ can be computed numerically by
using existing solvers for the level set method.
A brief note on the implementation of the algorithm is included in subsection A.3 in the Appendix.

There are a few things to note about the formulation in (\ref{eqn:BRT_final}).
First, Equation (\ref{eqn:BRT_final}) assumes that ${\mathbb{T}}$ is a desired goal state.
However, the formulation can be modified if ${\mathbb{T}}$ is an unsafe set, in which case, the adversarial modeling error tries to minimize the Hamiltonian.
Similarly, the input ${\bm{d}}$ can represent any other disturbance in the system, either adversarial or cooperative.
Second, to compute the FRT, the formulation in (\ref{eqn:BRT_final}) is modified from a final value PDE to an initial value PDE.

redpajama_set_name: RedPajamaArXiv