\subsection{Proof of Lemma \ref{lemma:recoveredsafeset}}\label{appendix:lemma_recoveredsafeset}
\begin{proof}
Assume for a contradiction that we compute $\Tilde{\mathcal{S}}$ as in Equation \eqref{eq:reachrecoveredsafeset} and that
\vspace{-0.25em}
\begin{equation}
\begin{aligned}
\small{
\Tilde{\mathcal{S}} \not\subseteq \text{BRT}^C}
\end{aligned}
\vspace{-0.25em}
\end{equation}
which implies
\vspace{-0.25em}
\begin{equation}
\begin{aligned}
\small{
\exists x^{'}: x^{'} \in \Tilde{\mathcal{S}}, x^{'} \not\in \text{BRT}^C}
\end{aligned}
\vspace{-0.25em}
\end{equation}
$x^{'} \in \Tilde{\mathcal{S}}$ implies that 
\vspace{-0.25em}
\begin{equation}
\begin{aligned}
\small{
\Tilde{V}(x^{'}, 0) > \delta}
\end{aligned}
\vspace{-0.25em}
\end{equation}
by Equation \eqref{eq:reachrecoveredsafeset}. 

On the other hand, $x^{'} \not\in \text{BRT}^C$ implies that $V(x^{'}, 0) \le 0$. 
We know that $\forall (x,t), J_{\Tilde{\policy}}(x,t) \le V(x, t)$ from Equation \eqref{eq:valuefunc}. 
Thus, $J_{\Tilde{\policy}}(x^{'},0) \le V(x^{'}, 0) \le 0$. 
Then 
\vspace{-0.25em}
\begin{equation}
\begin{aligned}
\small{
\Tilde{V}(x^{'}, 0) \le \delta}
\end{aligned}
\vspace{-0.25em}
\end{equation}
by Equation \eqref{eq:reachSafetyMetric}.

We have found the contradiction that $\Tilde{V}(x^{'}, 0) > \delta$ and $\Tilde{V}(x^{'}, 0) \le \delta$.
Thus, our assumption must be false, and it must be true that if we compute $\Tilde{\mathcal{S}}$ as in Equation \eqref{eq:reachrecoveredsafeset}, then $\Tilde{\mathcal{S}} \subseteq \text{BRT}^C$.

To prove the equivalent of Lemma \ref{lemma:recoveredsafeset} but in the case when $\mathcal{L}$ represents a set of desirable states, we use $\text{BRT}$ instead of $\text{BRT}^C$ and flip the value and cost inequalities.
\end{proof}

\subsection{Lemma \ref{lemma:scenarioOptimization}}\label{sec:lemma_scenarioOptimization}
\begin{lemma}
\label{lemma:scenarioOptimization}
Select a violation parameter $\epsilon \in (0, 1)$ and a confidence parameter $\beta \in (0, 1)$. Pick $N$ such that
\vspace{-0.25em}
\begin{equation}
\begin{aligned}
\small{
N \ge \frac{2}{\epsilon}\left(\ln{\frac{1}{\beta}}+1\right)}
\end{aligned}
\vspace{-0.25em}
\end{equation}
Then, with probability at least $1-\beta$, the solution $\hat{\delta}$ to Algorithm \ref{alg:scenarioOptimization} executed with $N$ satisfies the following condition.
\vspace{-0.25em}
\begin{equation} \label{eq:metricbound}
\begin{aligned}
\small{
\underset{\{x: x\in X, \Tilde{V}(x,0) > \hat{\delta}\}}{\mathbb{P}}\left(\left(\Tilde{V}(x,0): J_{\Tilde{\policy}}(x,0) \le 0\right) > \hat{\delta}\right) \le \epsilon}
\end{aligned}
\vspace{-0.25em}
\end{equation}
\end{lemma}
\begin{proof}
Lemma \ref{lemma:scenarioOptimization} follows from Theorem 1 in Scenario Optimization \cite{campi2009scenario}. To use the theorem, we need to prove the following conditions:
\begin{enumerate}
    \item computing $\delta$ can be converted into a standard Scenario Optimization problem
    \item Algorithm \ref{alg:scenarioOptimization} obtains $\hat{\delta}$ by sampling i.i.d from the same space that $\delta$ is optimized over, $X$, but with a probability distribution that is uniform over $\{x: x\in X, \Tilde{V}(x,0) > \hat{\delta}\}$
\end{enumerate}
$\delta$ can be formalized as the following optimization problem:
\begin{align*}
    &\min{g} \\
    &\text{s.t. } \forall x\in X, \left(\Tilde{V}(x,0): J_{\Tilde{\policy}}(x,0) \le 0\right) \le g
\end{align*}
This is a semi-infinite optimization problem where the constraints are linear, and thus convex, in the optimization variable $g$ for any given $x$. Thus, Lemma \ref{lemma:scenarioOptimization} follows from Theorem 1 in Scenario Optimization \cite{campi2009scenario} by replacing $c$ by $1$, $\gamma$ by $g$, $\Delta$ by the state space of interest $X$, and $f$ by $\left(\Tilde{V}(x,0): J_{\Tilde{\policy}}(x,0) \le 0\right) - g$. To use Theorem 1, we also require that i.i.d samples are chosen according to the uniform distribution over $\{x: x\in X, \Tilde{V}(x,t) > \hat{\delta}\}$. This can be proven by observing that, in Algorithm \ref{alg:scenarioOptimization}, the last iteration of the while loop samples $x$ randomly and independently from $\{x: x\in X, \Tilde{V}(x,t) > \hat{\delta}\}$ where $\hat{\delta}$ is exactly the final metric returned. Importantly, the algorithm does not update $\hat{\delta}$ in the final loop (since the loop subsequently breaks).

To prove the equivalent of Lemma \ref{lemma:scenarioOptimization} but in the case when $\mathcal{L}$ represents a set of desirable states, we flip the value and cost inequalities. We also take a maximum instead of a minimum and flip the inequality with $g$ when we formalize $\delta$ as an optimization problem.
\end{proof}

\subsection{Proof of Theorem \ref{theorem:scenarioOptimization}}\label{appendix:theorem_scenopt}
\begin{proof}
We rewrite the LHS of Inequality \eqref{eq:metricbound} of Lemma \ref{lemma:scenarioOptimization} in Appendix\ref{sec:lemma_scenarioOptimization} by substituting in $\hat{\mathcal{S}}$ by Equation \eqref{eq:reachapproximatelyrecoveredsafeset}:
\vspace{-0.25em}
\begin{align*}
\underset{\{x: x\in X, \Tilde{V}(x,0) > \hat{\delta}\}}{\mathbb{P}}&\left(\left(\Tilde{V}(x,0): J_{\Tilde{\policy}}(x,0) \le 0\right) > \hat{\delta}\right) \\
= \underset{\{x: x\in X, \Tilde{V}(x,0) > \hat{\delta}\}}{\mathbb{P}}&\left(\left(\Tilde{V}(x,0) > \hat{\delta}\right) \cap \left(J_{\Tilde{\policy}}(x,0) \le 0\right)\right) \\
= \underset{\{x: x\in X, \Tilde{V}(x,0) > \hat{\delta}\}}{\mathbb{P}}&\left(J_{\Tilde{\policy}}(x,0) \le 0\right) \\
= \underset{x \in \hat{\mathcal{S}}}{\mathbb{P}}&\left(J_{\Tilde{\policy}}(x,0) \le 0\right)
\end{align*}
\vspace{-0.25em}
where the second-to-last line follows because the probability is over the set of states that are already given to have a learned value greater than $\hat{\delta}$.

We know that $\forall (x,t), J_{\Tilde{\policy}}(x,t) \le V(x, t)$ from Equation \eqref{eq:valuefunc}. Thus, it follows that:
\vspace{-0.25em}
\begin{equation}
\begin{aligned}
\small{
\underset{x \in \hat{\mathcal{S}}}{\mathbb{P}}\left(V(x,0) \le 0\right) \le \underset{x \in \hat{\mathcal{S}}}{\mathbb{P}}\left(J_{\Tilde{\policy}}(x,0) \le 0\right) \le \epsilon}
\end{aligned}
\vspace{-0.25em}
\end{equation}

To prove the equivalent of Theorem \ref{theorem:scenarioOptimization} but in the case when $\mathcal{L}$ represents a set of desirable states, we flip the value and cost inequalities.
\end{proof}

\section{Background: Hamilton-Jacobi (HJ) Reachability and DeepReach} \label{sec:hj_reachability}
In this work, we build upon Hamilton-Jacobi reachability analysis to compute an approximation of the safe set. Here, we provide a quick overview of Hamilton-Jacobi reachability analysis and a specific toolbox to compute high-dimensional reachable sets, DeepReach.

\subsection{Hamilton-Jacobi (HJ) Reachability.} 
In HJ reachabilty, the computation of $\text{BRT}$ is formulated as an optimal control problem. We will explain it in the context of $\mathcal{L}$ being a set of undesirable states. In the end, we will comment on the case where $\mathcal{L}$ represents a set of desirable states. 

To construct this problem, we first define a target function $l(x)$ such that the sub-zero level of $l(x)$ yields $\mathcal{L}$:
\vspace{-0.25em}
\begin{equation}
\small{
    \mathcal{L} = \{x: l(x) \le 0\}}
    \vspace{-0.25em}
\end{equation}
$l(x)$ is commonly a signed distance function to $\mathcal{L}$. For example, we can choose $l(x) = \sqrt{p_x^2+p_y^2}-R$ for our Dubins3D Avoid running example.

Next, we define the cost function of a state corresponding to some policy $u(\cdot)$ to be the minimum of $l(x)$ over its trajectory:
\vspace{-0.75em}
\begin{equation}
    \label{eq:costfunctional}
    \small{
    J_{u(\cdot)}(x,t) = \min_{\tau \in [t, T]} l(\traj_{\state,\tvar}^{\ctrl}(\tau)).}
    \vspace{-0.25em}
\end{equation}
\normalsize
Since the system wants to avoid $\mathcal{L}$, our goal is to maximize $J_{u(\cdot)}(x,t)$.  
Thus, the value function corresponding to this optimal control problem is:
\vspace{-0.25em}
\begin{equation}
\small{
    \label{eq:valuefunc}
    V(x,t) = \sup_{u(\cdot)} \{J_{u(\cdot)}(x,t)\}.}
    \vspace{-0.25em}
\end{equation}
\normalsize

By defining our optimal control problem in this way, we can recover $\text{BRT}$ using the value function. In particular, the value function being sub-zero implies that the target function is sub-zero somewhere along the optimal trajectory, or in other words, has reached $\mathcal{L}$. Thus, $\text{BRT}$ is given as the sub-zero level set of the value function at the initial time:

\small
\vspace{-0.8em}
\begin{equation}
    \label{eq:BRT_from_valfunc}
    \text{BRT} = \{x: x\in X, V(x,0) \le 0 \}
\end{equation}
\normalsize
The value function in Equation \eqref{eq:valuefunc} can be computed using dynamic programming, resulting in the following final value Hamilton-Jacobi-Bellman Variational Inequality (HJB-VI):

\small
\vspace{-0.5em}
\begin{equation}
\begin{aligned}
    \label{eq:HJBVI}
    \min\Big\{D_t V(x,t)+ H(x,t), l(x)-V(x,t)\Big\} = 0,
    \end{aligned}
    \vspace{-0.25em}
\end{equation}
\normalsize
with the terminal value function $V(x,T) = l(x)$. 
$D_t$ and $\nabla$ represent the time and spatial gradients of the value function. 
$H$ is the Hamiltonian that encodes the role of dynamics and the optimal control.

\small
\vspace{-0.5em}
\begin{equation}
    \label{eq:ham}
    \begin{aligned}
    H(x,t) = \max_u & \langle \nabla V(x,t), \dyn(x,u)\rangle.
        \end{aligned}
    \vspace{-0.25em}
\end{equation}
\normalsize

The value function in Equation \eqref{eq:valuefunc} induces the optimal safety controller:

\small
\vspace{-0.25em}
\begin{equation}
    \label{eq:opt_ctrl}
    \begin{aligned}
    u^*(x,t) = \underset{u}{\arg\max} & \langle \nabla V(x,t), \dyn(x,u)\rangle.
        \end{aligned}
    \vspace{-0.25em}
\end{equation}
\normalsize
Intuitively, the safety controller aligns the system dynamics in the direction of the value function's gradient, thus steering the system towards higher-value states.

We have just explained the case where $\mathcal{L}$ represents a set of undesirable states. When the system instead wants to reach $\mathcal{L}$, a infimum is used instead of supremum in Equation $\eqref{eq:valuefunc}$.
The control wants to reach $\mathcal{L}$, hence there is a minimum instead of a maximum in Equation \eqref{eq:ham} and Equation \eqref{eq:opt_ctrl}.

Traditionally, the value function is computed by solving the HJB-VI over a discretized grid in the state space. Unfortunately, doing so involves computation whose memory and time complexity scales exponentially with respect to the system dimensionality, making these methods practically intractable for high-dimensional systems, such as those beyond 5D. Fortunately, a recent deep learning approach, \text{DeepReach}{}, has been proposed to enable Hamilton-Jacobi reachability for high-dimensional systems. 

\subsection{\text{DeepReach}{} Approximate Solutions.} 
Instead of solving the HJB-VI over a grid, DeepReach represents the value function as a sinusoidal deep neural network (DNN) to learn a parameterized approximation of the value function \cite{bansal2021deepreach}. Thus, memory and complexity requirements for training scale with the value function complexity rather than grid resolution.
To train the DNN, DeepReach uses self-supervision on the HJB-VI itself. 
Ultimately, it takes as input a state $x$ and time $t$, and it outputs a learned value $\Tilde{V}(x,t)$.
$\Tilde{V}(x,t)$ also induces a corresponding policy $\Tilde{\pi}(x,t)$.
We refer interested readers to \cite{bansal2021deepreach} for further details.

The learned $\Tilde{V}(x,t)$ from DeepReach can be used to obtain a BRT as in Equation \eqref{eq:BRT_from_valfunc}, but it will only be as accurate as $\Tilde{V}(x,t)$. Unfortunately, like any learning method, $\Tilde{V}(x,t)$ can deviate substantially from the true $V(x,t)$. 

Our goal is to recover the biggest safe set $\hat{\mathcal{S}}$ using $\Tilde{V}(x,t)$ that we can (probabilistically) guarantee to be fully contained within the true safe set.
Although we work with \text{DeepReach}{} solutions in particular for this problem setup, our proposed approach in Section \ref{sec:approach} can verify any general $\Tilde{V}(x,t)$ and $\Tilde{\pi}(x,t)$, regardless of whether \text{DeepReach}{}, a level-set method, or some other tool is used to obtain them.
 
\section{Discussion and Future Work}
In this work, we present an approach to compute an error bound for DeepReach solutions to recover a provably safe approximation of the true reachable tube.
This allows us to utilize the power of learning-based reachability methods for providing safety assurances for high-dimensional dynamical systems.
We apply our method to obtain reachable tubes for high-dimensional rocket-landing and multi-vehicle collision-avoidance problems which traditional methods struggle with. 

In the future, we will explore a more refined approach to error correction (as opposed to a uniform correction), as discussed briefly in the end of Section \ref{sec:approach}. 
Other directions include additionally considering worst-case disturbances in the system dynamics and using the verification approaches proposed here for a targeted refinement of DeepReach solutions.

\section{Introduction}
Imagine the real-world motivating example of an aerial vehicle that wants to reach a destination while avoiding any collisions. There might even be adversarial disturbances, such as wind or another aircraft nearby, that have the potential to hinder the vehicle from achieving its objective. Since such autonomous systems are becoming increasingly common in our everyday lives, it is vital that we can compute controllers with provable safety guarantees. There already exist many approaches in the literature for computing these. For example, ... However, they are limited by ... The popular Hamilton-Jacobi reachability approach, on the other hand, addresses many of these limitations.

\textbf{Hamilton-Jacobi Reachability Methods.} Hamilton-Jacobi (HJ) reachability analysis is a popular formal verification method for autonomous systems that computes both a safe controller and its corresponding safe state configurations. Its advantages over other approaches includes its ability to handle general nonlinear system dynamics, bounded adversarial system disturbances, and state and input constraints. 

In reachability analysis, we are interested in the \textit{Backward Reachable Tube (BRT)} of a dynamical system, which is defined in the following Section \ref{sec:problem_setup} in Equation \eqref{eq:reachBRT}, since it directly relates to system safety.
Hamilton-Jacobi reachability analysis computes a system's BRT by solving a Hamilton-Jacobi (HJ) partial differential equation (PDE) \eqref{eq:HJIVI}, as explained in Section \ref{sec:problem_setup}. Level set methods can compute BRTs accurately up to numerical precision on a grid representing a discretization of the state space \cite{mitchell2005time, lygeros2004reachability, Margellos11, fisac2015reach, bansal2017hamilton}.
Unfortunately, this process results in an exponential scaling of memory and computational complexity with respect to system dimensionality, making it challenging to scale beyond 5D systems \cite{bansal2017hamilton}.

Several methods have been proposed to overcome this challenge by trading off between the class of dynamics they can handle, the approximation quality of the BRT, and the required computation. For example, the methods presented in \cite{Greenstreet1998, Frehse2011, Kurzhanski00, Kurzhanski02, Maidens13, girard2005reachability, althoff2010computing} analyze high-dimensional affine systems using sets of pre-specified shapes, sum-of-squares optimization \cite{Majumdar13, majumdar2017funnel} and parallelotopes \cite{Dreossi16} rely on specific dynamics for computing BRTs for polynomial dynamics \cite{henrion2014convex}; leverage monotonicity of dynamics \cite{coogan2015efficient}; use decoupled or reduced order dynamics to decompose the BRT computation \cite{chen2018decomposition, holmes2020reachable, liu2020leveraging, Chen2016d, somil2020sequential}; leverage the Hopf-Lax formula to compute BRTs for convex dynamics \cite{chow2017algorithm}; or deal with affine dynamics \cite{bak2019numerical, Nilsson2016}.
Control Lyapunov (and Barrier) functions are also used to design safety controller \cite{ames2016control, ames2019control}; however, finding feasible Lyapunov functions can be challenging, especially in the presence of control or disturbance bounds.

Sampling and learning-based methods have also been explored for computing BRTs \cite{liebenwein2018sampling}, including scenario optimization \cite{sutter2017data, devonport2020estimating, ghosh2019new} and adversarial sampling \cite{lew2020sampling}, Support Vector Machines \cite{allen2014machine} and DNNs \cite{rubies2019classification, fisac2019bridging, darbon2020overcoming, niarchos2006neural, djeridane2006neural} leveraging data to compute approximate BRTs for nonlinear systems; however, they often struggle in the presence of adversarial disturbances and state constraints.
Other methods require explicit supervision for learning which is hard to generate for general high-dimensional systems.

Unfortunately, all these methods either need strong assumptions on dynamics, cannot handle presence of exogenous inputs or hard constraints, or lead to a very conservative approximation of the reachable set.
Additionally, very few of these approaches can handle adversarial inputs; yet, considering adversarial inputs (e.g., other agents or the model uncertainty) and state constraints is critical in analysis 
and synthesis of safe controllers. Thus, in this work, we focus on a recently developed deep learning-based approach to approximately solve high-dimensional Hamilton-Jacobi reachability problems: \text{DeepReach}{}, which overcomes all of these limitations.

\textbf{DeepReach.} What sets \text{DeepReach}{} apart is its ability to compute BRTs as well as the \textit{corresponding safety controller} for \textit{general} nonlinear dynamical systems in the presence of \textit{disturbances} and \textit{state and input constraints}. 
In contrast to level-set methods, DeepReach learns continuous, rather than discretized, solutions to high-dimensional systems which would otherwise be intractable. 

Instead of solving this PDE over a grid, \text{DeepReach}{} draws inspiration from the recent progress in solving PDEs using deep learning, and represents the value function as a sinusoidal deep neural network (DNN) to learn a parameterized approximation of the value function via self-supervision on the PDE. 
\text{DeepReach}{} can also be interpreted as an approximate dynamic programming method \cite{bertsekas1995dynamic} where a DNN is used to \textit{learn} the basis features for the value function. 
Thus, the computation and memory requirements for obtaining the value function do not scale with the grid resolution, but rather the complexity of the value function.
This enables \text{DeepReach}{} to scale to general, high-dimensional dynamical systems.

Even if we can now obtain an approximate BRT for high-dimensional systems for the first time with DeepReach, they are limited in usefulness by their approximate nature. Recall that one of the primary use cases of Hamilton-Jacobi reachability is to provide formal, not approximate, safety guarantees. Thus, our work bridges this gap by proposing a verification method to transform a DeepReach solution into a formal guarantee.

\textbf{Existing Verification Methods.}

\textbf{\text{Scenario Optimization Verification}{} Method.} 
\begin{enumerate}
    \item list specific contributions of this work
\end{enumerate}


\section{Approach}\label{sec:approach}
Our approach is to compute a minimal error correction to $\Tilde{V}(x,t)$ such that the corrected value function can be used to extract a (probabilistically) safe set.
We will explain our approach in the context of $\mathcal{L}$ being a set of undesirable states. In the end of each subsection, we will comment on the case where $\mathcal{L}$ represents a set of desirable states. In Section \ref{subsec:errormetric}, we propose an error metric for correcting the value function. In Section \ref{subsec:scenopt}, we propose a practical method to compute this error metric through scenario optimization.

\subsection{Error Metric}\label{subsec:errormetric}
We propose a new error metric $\delta_{\Tilde{\vfunc},\Tilde{\policy}}$ for the value function correction. It is defined as the maximum learned value of an empirically unsafe initial state under the induced policy $\Tilde{\pi}$:
\vspace{-0.25em}
\begin{equation}
\label{eq:reachSafetyMetric}
\begin{aligned}
\small{
\delta_{\Tilde{\vfunc},\Tilde{\policy}} \coloneqq \max_{x\in X}\{\Tilde{V}(x,0): J_{\Tilde{\policy}}(x,0) \le 0\}},
\end{aligned}
\end{equation}
where $J_{\Tilde{\policy}}(x,0)$ is the cost function associated with the trajectory obtained by using the policy $\Tilde{\pi}(x,t)$ from an initial state $x$ and initial time $t=0$ until $T$ (see Equation \eqref{eq:costfunctional}).
Here on, we use $\delta$ as a shorthand for $\delta_{\Tilde{\vfunc},\Tilde{\policy}}$ for brevity purposes.

Intuitively, $\delta$ finds the tightest level of the learned value function that separates the states that are safe under the induced policy $\Tilde{\pi}(x,t)$ from the ones that are not.
Thus, any initial state within the super-$\delta$ level set of $\Tilde{V}(x,0)$ is guaranteed to be safe under the (possibly sub-optimal) policy $\Tilde{\pi}(x,t)$. 
Lemma \ref{lemma:recoveredsafeset} formalizes this claim.
\begin{lemma}
\label{lemma:recoveredsafeset}
If we compute $\Tilde{\mathcal{S}}$ as:
\vspace{-0.25em}
\begin{equation}
\label{eq:reachrecoveredsafeset}
\begin{aligned}
\small{
\Tilde{\mathcal{S}} = \{x\in X: \Tilde{V}(x,0) > \delta\}}
\end{aligned}
\vspace{-0.25em}
\end{equation}
then $\Tilde{\mathcal{S}} \subseteq \text{BRT}^C$.
\end{lemma}

The proof for Lemma \ref{lemma:recoveredsafeset} can be found in Appendix\ref{appendix:lemma_recoveredsafeset}.
Intuitively, Lemma \ref{lemma:recoveredsafeset} states that if we could exactly compute $\delta$, then we can recover a $\Tilde{\mathcal{S}}$ that is guaranteed to be a subset of the true safe set. 
Furthermore, it is easy to see that, by definition, $\delta$ is the smallest (uniform) value adjustment required on $\Tilde{V}(x,t)$ to guarantee its safety. 
Thus, $\Tilde{\mathcal{S}}$ is the largest safe set we can recover from the learned value function by such a uniform error correction procedure.

Unfortunately, computing $\delta$ is a challenging optimization problem, since both the value function and the cost function in the metric definition \eqref{eq:reachSafetyMetric} are typically non-convex functions of $x$. 
In the next section, we propose a scenario optimization-based approach to compute a high-confidence bound on $\delta$ and thereby generate a high-confidence safe set.

\begin{remark}
We assume there exists some unsafe states in the system so that $\delta_{\Tilde{\vfunc},\Tilde{\policy}}$ is well-defined. $\delta_{\Tilde{\vfunc},\Tilde{\policy}}$ otherwise can trivially be defined as $-\infty$.
\end{remark}

So far, we have discussed an error correction metric for the case where $\mathcal{L}$ represents a set of undesirable states. 
When $\mathcal{L}$ represents a set of desirable states, we use a minimum instead of a maximum and flip the cost inequality in the metric definition \eqref{eq:reachSafetyMetric}.
To recover the safe set, we extract the sub-$\delta$ level set of the value function in Equation \eqref{eq:reachrecoveredsafeset}.

\subsection{\text{Scenario Optimization Verification}{} Method}\label{subsec:scenopt}
We will compute a high-confidence bound on $\delta$ by utilizing a random sampling procedure referred to as \textit{scenario optimization} in the systems and control design literature \cite{campi2009scenario}. The \text{Scenario Optimization Verification}{} Method to compute an approximation $\hat{\delta}$ is summarized in Algorithm \ref{alg:scenarioOptimization}.

\setlength{\textfloatsep}{0pt}
\begin{algorithm}
\caption{\text{Scenario Optimization Verification}{}}
\label{alg:scenarioOptimization}
\begin{algorithmic}[1]
\REQUIRE $X$, $N$, $M$, $\Tilde{V}(x,0)$, $J_{\Tilde{\policy}}(x,0)$
\STATE $\delta_0 \gets -\infty$
\FOR{$i= 0, 1, \ldots, M-1$} 
    \STATE $\mathcal{D}_i \gets$ Sample $N$ states IID from $\{x:x\in X, \Tilde{V}(x,0) > \delta_i\}$ \label{line:sample}
    \IF{$\exists x \in \mathcal{D}_i : J_{\Tilde{\policy}}(x,0) \le 0$}
        \STATE $\delta_{i} \gets \max_{x \in \mathcal{D}_i}\{\Tilde{V}(x,0): J_{\Tilde{\policy}}(x,0) \le 0\}$ \label{line:assign}
    \ELSE
        \STATE \textbf{break}
    \ENDIF
\ENDFOR
\RETURN $\hat{\delta} := \delta_i$
\end{algorithmic}
\end{algorithm}

At a high-level, Algorithm \ref{alg:scenarioOptimization} computes a converging sequence of $\delta_i$ that approximates $\delta$ via random sampling until we no longer find any safety violations or reach a maximum number of iterations, $M$.
Specifically, at each iteration $i$, we randomly sample $N$ initial states within the super-$\delta_i$ level set of $\Tilde{V}$ (Line \ref{line:sample}) and compute the costs of associated trajectories $J_{\Tilde{\policy}}(x,0)$ obtained by using $\Tilde{\pi}(x,t)$ as the controller (Equation \eqref{eq:costfunctional}). We next compute the maximum learned value among the states that violate safety and use that as the new estimate of $\delta$ (Line \ref{line:assign})), essentially discarding a level region of $\Tilde{V}(x, 0)$ that were empirically found to be unsafe under $\Tilde{\pi}$.
Thus, with each iteration, we obtain a tighter approximation of $\delta$.
We repeat these steps until no more safety violations are encountered or a maximum number of iterations are achieved. 

While the method described in Algorithm \ref{alg:scenarioOptimization} is simple to understand and execute, scenario optimization provides a formal guarantee for bounding $\delta$ with $\hat{\delta}$.
Thus, we can use the approximated $\hat{\delta}$ instead of $\delta$ to compute an approximated safe set $\hat{\mathcal{S}}$ similar to the one defined in Lemma \ref{lemma:recoveredsafeset}:
\vspace{-0.25em}
\begin{equation}
\label{eq:reachapproximatelyrecoveredsafeset}
\begin{aligned}
\small{
\hat{\mathcal{S}} = \{x: x\in X, \Tilde{V}(x,0) > \hat{\delta}\}}
\end{aligned}
\vspace{-0.25em}
\end{equation}
Crucially, we can make a formal probabilistic safety guarantee for $\hat{\mathcal{S}}$. This is summarized in Theorem \ref{theorem:scenarioOptimization} and proven in Appendix\ref{appendix:theorem_scenopt}.
\begin{theorem}[Scenario Optimization Verification Theorem]
\label{theorem:scenarioOptimization}
Select a violation parameter $\epsilon \in (0, 1)$ and a confidence parameter $\beta \in (0, 1)$. Pick $N$ such that
\vspace{-0.25em}
\begin{equation} \label{eq:prescribedN}
\begin{aligned}
\small{
N \ge \frac{2}{\epsilon}\left(\ln{\frac{1}{\beta}}+1\right)}
\end{aligned}
\vspace{-0.25em}
\end{equation}
Then, with probability at least $1-\beta$, the recovered safe set $\hat{\mathcal{S}}$ from Equation \eqref{eq:reachapproximatelyrecoveredsafeset} satisfies the following condition:
\begin{equation}\label{eq:guarantee}
\begin{aligned}
\small{
\underset{x \in \hat{\mathcal{S}}}{\mathbb{P}}(V(x,0) \le 0) \le \epsilon}
\end{aligned}
\vspace{0.50em}
\end{equation}
\end{theorem}

Note that here we state our results for the case where Algorithm \ref{alg:scenarioOptimization} terminates prior to reaching the maximum number of iterations (i.e., the sequence $\delta_i$ converges).
Similar results can be proved for the other case as well (and are not here, for the lack of space). 
In practice, we find that Algorithm \ref{alg:scenarioOptimization} converges within 3-4 iterations.

Disregarding the confidence parameter $\beta$ for a moment, the above theorem states that the volume of the unsafe states within $\hat{\mathcal{S}}$ is smaller than or equal to the prescribed $\epsilon$ value, as long as we sample enough states (to satisfy Inequality \eqref{eq:prescribedN}) during the computation of $\hat{\delta}$. 
As $\epsilon$ approaches $0$, the number of safety violations in the recovered safe set also approaches $0$. In
turn, the simulation effort grows unbounded since $N$ is inversely
proportional to $\epsilon$.

To interpret the confidence parameter $\beta$, note that $\hat{\mathcal{S}}$ is a random variable that depends on a randomly sampled set of initial states. It may be the case that we just happen to draw a poorly representative sample, in which case the $\epsilon$ bound does not hold. $\beta$ controls the probability of this adverse event happening, which regards the correctness of the final guarantee in Equation \eqref{eq:guarantee}. Fortunately, $N$ only grows logarithmically with $\frac{1}{\beta}$, so $\beta$ can be chosen to be an extremely small value such as $10^{-12}$. $1-\beta$ should then be so close to $1$ that it does not have any practical importance.

When $\mathcal{L}$ represents a set of desirable states, we initialize $\delta_0$ to be infinity instead of negative infinity, flip the value and cost inequalities, and take a minimum instead of maximum in Algorithm \ref{alg:scenarioOptimization}. Also, we flip the value inequalities in Equation \eqref{eq:reachapproximatelyrecoveredsafeset} and Equation \eqref{eq:guarantee}. 
\vspace{0.2em}

\noindent \textit{\textbf{Running example: Dubins3D Avoid.}} We now demonstrate the \text{Scenario Optimization Verification}{} method on the Dubins3D Avoid running example introduced in Section \ref{sec:problem_setup}. We choose system parameters $v=0.6 m/s, u_{\min}=-1.1 rad/s, u_{\max}=1.1 rad/s, R=0.25m$, and use DeepReach to learn the system's value function. 
The used neural network parameters and architecture are the same as in the DeepReach paper \cite{bansal2021deepreach}, i.e., a DNN with 3 hidden layers and 512 neurons in each layer.

For this first example, we purposefully sabotage the training of DeepReach by permitting it to train for only 8,000 epochs, far less than $100$K epochs used in the DeepReach paper. 
We do this to illustrate our method's utility even with low-quality learned solutions. 
We apply our method to a well-trained DeepReach solution in Section \ref{sec:dubins3davoid}.

The overall training took less than 30 minutes on an NVIDIA GeForce RTX 3090Ti.
We execute our verification approach choosing $\beta = 10^{-16}, \epsilon = 10^{-3}$, resulting in an $N = 75683$ as per Theorem 1.
Thus, we will be $1-10^{-16}$ confident that at least $1-10^{-3}$ of our recovered set will be provably safe.
Algorithm \ref{alg:scenarioOptimization} converges in 4 iterations in this case, resulting in a $\hat{\delta} = 0.3728$.
Slices of the trained BRT (yellow), the recovered BRT (blue), and the ground truth BRT boundary (black) are shown in Figure \ref{fig:baddubinsavoid}. 
The ground truth is computed by a state-of-the-art PDE solver, \textit{Level Set Toolbox (LST)} \cite{LSToolbox}, which computes the value function over a discrete grid of size $101\times 101\times 101$. 

\begin{figure}[ht]
    \vspace{-0.65em}
    \centering
    \includegraphics[width=1 \columnwidth]{initial_sub/figures/Figure_1.png}
    \vspace{-1.5 em}
    \caption{Dubins3D Avoid: Slices of the trained, recovered, and ground truth border BRTs for three values of $\theta$ from a sabotaged DeepReach solution. The recovered BRT completely encompasses the ground truth BRT and provides a provably safe approximation of the safe set.}
    \vspace{-0.65em}
    \label{fig:baddubinsavoid}
\end{figure}

As evident from the figure, the recovered BRT slices completely encompass the true unsafe states (the blue region encompasses the states within the black boundary). 
Hence, the complement of the recovered BRT is safe, as we expect based on Theorem \ref{theorem:scenarioOptimization}. 
In the rightmost slice, observe that the trained BRT is smaller than the ground truth BRT and is thus unsafe as it is. 
Our verification approach correctly expands the trained BRT past the ground truth boundary, ensuring safety. 
Validation by sampling 1M states within the recovered safe set reveals a violation rate of $3\times 10^{-6} << \epsilon$.

Note that the trained solution's safety violations are concentrated near $\theta=\pi$. 
It is likely these errors which are responsible for the large expansion of the recovered BRT even for slices that are almost accurate (left and middle slices in Figure \ref{fig:baddubinsavoid}).
This is because our method expands the BRT uniformly.
However, too much expansion is undesirable, since we want to recover as much of the safe set as possible. 

One way to overcome this challenge is to apply refined error correction to $\Tilde{V}(x,t)$ by independently verifying separate regions of the states-space (i.e., separately consider each region as our $X$).
This allows us to recover sets with the same safety guarantee but expanded differently depending on their own region's error. 
However, this method's efficacy depends heavily on the choice of sub-regions, which are hard to prespecify, especially for high-dimensional systems.
Instead, we can resort to a data-driven approach.

Specifically, we trained a simple multilayer perceptron (MLP) on 1M randomly sampled initial states and their empirical safety violation costs (the training took an hour on a standard GPU). 
Even though the predictor may be inaccurate (after all, predicting the safe set is the original challenging problem), we hypothesize that it will have learned something useful enough about the general distribution of errors throughout the state space. 
We divided the output of the MLP in 10 bins, each corresponding to a different range of empirical cost.
Thus, we expect the last bin to correspond to the states with maximal safety violations (in this case, that will correspond to states with $\theta=\pi$).
We next ran Algorithm \ref{alg:scenarioOptimization} independently for each bin to compute an error correction.
The corresponding results are shown in Figure \ref{fig:binneddubinsavoid}.
As evident, we get a much tighter approximation of the safe set near the slices that have small errors, allowing us to recover a much bigger safe set overall.
In this case, the recovered safe set volume increased by 35\%. 

\begin{figure}[ht]
    \vspace{-0.65em}
    \centering
    \includegraphics[width=1 \columnwidth]{initial_sub/figures/Figure_3.png}
    \vspace{-1.5em}
    \caption{Dubins3D Avoid: Slices of the trained, recovered, and ground truth border BRTs for three values of $\theta$ from a sabotaged DeepReach solution when MLP binning is used to refine recovery. The recovered BRT completely encompasses the ground truth BRT and provides a provably safe approximation of the safe set that is larger than when binning is not used.}
    \vspace{-0.65em}
    \label{fig:binneddubinsavoid}
\end{figure}

As discussed earlier, the performance improvement of this approach depends heavily on the selection of subregions.
We defer this to future work, and for the rest of this paper, we only focus on uniform value function correction.
\section{Problem Setup} \label{sec:problem_setup}
Consider a dynamical system with state $x \in X \subseteq \mathbb{R}^n$, control $u \in \mathcal{U}$, and dynamics $\dot{x} = \dyn(x, u)$ governing how $x$ evolves over time until a final time $T$. 
Let $\traj_{\state,\tvar}^{\ctrl}(\tau)$ denote the state achieved at time $\tau \in [t, T]$ by starting at initial state $x$ and time $t$ and applying control $u(\cdot)$ over $[t,\tau]$. 
Let $\mathcal{L}$ represent a target set of states that the agent wants to either reach (e.g. goal states) or entirely avoid (e.g. obstacles).
\vspace{0.2em}

\noindent \textit{\textbf{Running example: Dubins3D Avoid.}} As a running example, consider a simple low-dimensional system in the literature known as Dubins3D. It involves a car with position $(p_x, p_y)$, heading $\theta$, velocity $v$, and steering control $u_1 \in [u_{\min}, u_{\max}]$. The car's state $x$ evolves according to
\begin{align*}
    \dot{p_x} = v\cos{\theta}, \quad \dot{p_y} = v\sin{\theta}, \quad
    \dot{\theta} = u_1
\end{align*}
It wants to avoid a circle with radius $R$ centered at the origin. Thus, we define the target set $\mathcal{L}$ of this system to be:
\begin{align*}
    \mathcal{L} = \{x: \sqrt{p_x^2 + p_y^2} \le R\}
\end{align*}
We use Dubins3D as a benchmark because it is tractable for traditional reachability methods to compute a solution for, and we can thus compare our results with a ground truth. 
In Section \ref{sec:examples}, we present high-dimensional problems which traditional methods struggle with.

In this setting, we are interested in computing the system's initial-time Backward Reachable Tube within our state space of interest $X$, which we denote as $\text{BRT}$. 
We define $\text{BRT}$ as the set of all initial states in $X$ from which the agent will eventually reach $\mathcal{L}$ within the time horizon $[0,T]$, despite best control efforts:
\vspace{-0.25em}
\begin{equation}
\label{eq:avoidBRT}
\begin{aligned}
\small{
\text{BRT} = \{x: x\in X, \forall u(\cdot), \exists \tau \in [0, T], \traj_{\state,0}^{\ctrl}(\tau) \in \mathcal{L}\}}
\end{aligned}
\vspace{-0.25em}
\end{equation}
When $\mathcal{L}$ represents unsafe states for the system, staying outside of $\text{BRT}$ is desirable. 
When $\mathcal{L}$ instead represents the states that the agent wants to reach (e.g., a goal set), $\text{BRT}$ is defined as the set of all initial states in $X$ from which the agent, acting optimally, can eventually reach $\mathcal{L}$ within $[0,T]$. Thus, staying within $\text{BRT}$ is desirable.

Our goal in this work is to compute a provable approximation of the safe set. Specifically, we want to compute an approximation $\hat{\mathcal{S}}$ such that $\hat{\mathcal{S}}\subseteq \text{BRT}^C$ ($\hat{\mathcal{S}}\subseteq \text{BRT}$ when $\mathcal{L}$ represents goal states). 
We are especially interested in settings where the system is high-dimensional, with which current state-of-the-art methods struggle.

\section{Case Studies} \label{sec:examples}
We will evaluate our approach on various reachability problems.
First, we show results for the avoid version and the reach version of the low-dimensional Dubins3D system for which we have a ground truth BRT available for comparison. 
Next, we recover safe sets for a 9D multivehicle collision avoidance problem and for a 6D rocket landing problem, with which traditional methods struggle.
Like our running example, we use $\beta = 10^{-16}, \epsilon = 10^{-3}, N = 75683$.

\subsection{Dubins3D Avoid} \label{sec:dubins3davoid}
In Figure \ref{fig:dubinsavoid}, we show the results of applying our method to a well-trained DeepReach solution for the same Dubins3D Avoid system as in the running example from Section \ref{sec:approach}.
The solution is trained to consider the periodicity of $\theta$ and for 100K epochs (as opposed to 8K for the running example solution).
The proposed algorithm results in a very small error correction of -0.0016. Interestingly, this shrinks the trained BRT to recover a larger safe set than initially learned.
Validation by sampling 1M states within the recovered safe set reveals a violation rate of $3.2\times 10^{-5} << \epsilon$. 

The fact that a much smaller error correction is found for the trained solution shown in Figure \ref{fig:dubinsavoid} than the one in Figure \ref{fig:baddubinsavoid} indicates it is of much higher quality.
This illustrates how the error correction metric $\delta$ can also be used to evaluate the quality of different approximate value function solutions.
For example, we can evaluate the relative performance of different DNN hyperparameters, which is especially helpful when the ground truth BRT is not available for comparison.  

\begin{figure}[ht]
    \vspace{-0.65em}
    \centering
    \includegraphics[width=1 \columnwidth]{initial_sub/figures/Figure_2.png}
    \vspace{-1.5em}
    \caption{Dubins3D Avoid: Slices of the trained, recovered, and ground truth border BRTs for three values of $\theta$ from a well-trained DeepReach solution. The trained BRT completely encompasses the recovered BRT, so we plot the recovered BRT border on top for visualization purposes.}
    \vspace{-0.65em}
    \label{fig:dubinsavoid}
\end{figure}

\subsection{Dubins3D Reach} \label{sec:dubins3dreach}
In Figure \ref{fig:dubinsreach}, we show the results of applying our method to a DeepReach solution trained with the same scheme as for the Dubins3D Avoid solution in Section \ref{sec:dubins3davoid}, but for the reach version. 
That is, $\mathcal{L}$ is now a set of desirable states, so a safe approximate BRT should be fully contained by the true BRT. 
The trained BRT is just barely crossing the ground truth boundary, and the proposed method cuts it down to just behind the boundary, ensuring system safety.
\begin{figure}[ht]
    \vspace{-0.65em}
    \centering
    \includegraphics[width=1 \columnwidth]{initial_sub/figures/Figure_4.png}
    \vspace{-1.5em}
    \caption{Dubins3D Reach: Slices of the trained, recovered, and ground truth border BRTs for three values of $\theta$ from a well-trained DeepReach solution. The recovered BRT is correctly a subset of the ground truth BRT.}
    \vspace{-0.65em}
    \label{fig:dubinsreach}
\end{figure}

\subsection{Multivehicle Collision Avoidance} \label{sec:multicollision}
We now consider a 9D collision avoidance system involving 3 independent Dubins3D cars. 
The $i$th car has position $(p_{xi}, p_{yi})$, heading $\theta_i$, velocity $v$, and steering control $u_i \in [u_{\min}, u_{\max}]$. 
The dynamics of vehicle $i$ are given as:
\begin{align*}
    \dot{p}_{xi} = v\cos{\theta_i}, \quad \dot{p}_{yi} = v\sin{\theta_i}, \quad
    \dot{\theta_i} = u_i
\end{align*}
The (undesirable) target set is given by the states where any of the vehicle pairs is in collision:
\begin{equation*}
    \mathcal{L} = \{x: \min\{d(Q_1, Q_2), d(Q_1, Q_3), d(Q_2, Q_3)\} \le R\}
\end{equation*}
where $d(Q_i, Q_j)$ is the distance between cars $i, j$.
We choose $v=0.6, u_{\min}=-1.1, u_{\max}=1.1, R=0.25$ for our case study.
The results of applying our method to a DeepReach solution are shown in Figure \ref{fig:multi}.

\begin{figure}[ht]
    \vspace{-0.65em}
    \centering
    \includegraphics[width=1 \columnwidth]{initial_sub/figures/Figure_5.png}
    \vspace{-1.5em}
    \caption{Multivehicle Collision Avoidance: Slices of the trained and recovered BRTs for three values of $\theta_1$ from a DeepReach solution.}
    \vspace{-0.65em}
    \label{fig:multi}
\end{figure}

This high-dimensional example is typically difficult to compute with traditional methods, yet here we demonstrate successful recovery of a safe set with formal safety guarantees.
While the recovery preserves much of the safe set that is learned, a significant amount is pruned off. 
The percent reduction in safe set size (calculated from 1M samples) is much larger in this example (26\%) than for the Dubins3D solutions in Section \ref{sec:dubins3davoid} (0\%) and Section \ref{sec:dubins3dreach} (10\%), suggesting that the learned value function is more inaccurate for higher dimensional systems. 

To validate safety of the recovered BRT, in Figure \ref{fig:multihist}, we plot the minimum pairwise spare distance between vehicles along trajectories 
spawning from states sampled within the trained and recovered safe sets. 
Note that the trained safe set contains states which result in a collision, as shown by the mass to the left of the dotted line representing the collision distance. After verification, the recovered safe set shows no such mass. 
Validation by sampling 1M samples in the recovered safe set reveals a violation rate of $5\times 10^{-6} << \epsilon$.

\begin{figure}[ht]
    \vspace{0.5em}
    \centering
    \includegraphics[width=1 \columnwidth]{initial_sub/figures/Figure_6.png}
    \vspace{-1.5em}
    \caption{Multivehicle Collision Avoidance: Histogram of the minimum pairwise spare distance between vehicles along trajectories which spawn from states sampled within the trained and recovered safe sets. The collision distance is indicated by the dotted line.}
    \vspace{0.25em}
    \label{fig:multihist}
\end{figure}

\subsection{Rocket Landing} \label{sec:rocketlanding}
Consider a 6D rocket landing system with position $(p_{x}, p_{y})$, heading $\theta$, velocity $(v_{x}, v_{y})$, angular velocity $\omega$, and torque controls $\tau_1, \tau_2\in [-250, 250]$. 

\noindent The dynamics are:
\begin{gather*}
    \dot{p_x} = v_x,~\dot{p_y} = v_y, ~\dot{\theta} = \omega,~\dot{\omega} = 0.3\tau_1,\\ \dot{v_x} = \tau_1 \cos{\theta} - \tau_2 \sin{\theta},~
    \dot{v_y} = \tau_1 \sin{\theta} + \tau_2 \cos{\theta} - g, 
\end{gather*}
where $g = 9.81$ is acceleration due to gravity.
The target set is the set of states where the rocket reaches a rectangular landing zone of side length 20 centered at the origin:
\begin{align*}
    \mathcal{L} = \{x: |p_x| < 20.0, p_y < 20.0\}
\end{align*}

In Figure \ref{fig:rocket}, we show the results of applying our method to a DeepReach solution. 
We also plot the trajectories emanating from 20 randomly sampled initial states within the recovered safe set, demonstrating that they do indeed reach the target set (the green region) by following the induced policy. Validation by sampling 1M states within the recovered safe set reveals a violation rate of $5\times 10^{-6} << \epsilon$.

\begin{figure}[ht]
    \vspace{-0.65em}
    \centering
    \includegraphics[width=1 \columnwidth]{initial_sub/figures/Figure_7.png}
    \vspace{-1.5em}
    \caption{Rocket Landing: Slices of the trained and recovered BRTs for three values of $v_y$ from a DeepReach solution. The system trajectories safely reach the target set (green) from the initial states within the recovered BRT.}
    \vspace{-0.65em}
    \label{fig:rocket}
\end{figure}

\section{Introduction}
As autonomous systems become increasingly common in our everyday lives, it is vital that we can design provably safe controllers for these systems.
Hamilton-Jacobi (HJ) reachability analysis provides a powerful framework to design such controllers for general nonlinear dynamical systems \cite{lygeros2004reachability, mitchell2005time}.
In reachability analysis, the safe states for the system are characterized through the \textit{Backward Reachable Tube (BRT)} of the system.
This is the set of states that spawn trajectories that will eventually reach some given target set despite the best control effort.
Thus, if the target set represents the set of undesirable states, the BRT represents unsafe states for the system and should be avoided. 
Along with the BRT, reachability analysis also provides a safety controller to keep the system outside the BRT. 

Traditionally, the BRT computation in HJ reachability is formulated as an optimal control problem.
The BRT can then be obtained as a sub-zero level solution of the corresponding value function.
Obtaining the value function requires solving a partial differential equation (PDE) over a state-space grid, resulting in an exponentially scaling computation complexity with the number of states \cite{bansal2017hamilton}.
To overcome this challenge, a variety of solutions have been proposed that trade off between the class of dynamics
they can handle, the approximation quality of the BRT,
and the required computation.
These include specialized methods for linear and affine dynamics \cite{Greenstreet1998, Frehse2011, Kurzhanski00, Kurzhanski02, Maidens13, girard2005reachability, althoff2010computing, bak2019numerical, Nilsson2016}, polynominal dynamics \cite{Majumdar13, majumdar2017funnel, Dreossi16, henrion2014convex}, monotonic dynamics \cite{coogan2015efficient}, and convex dynamics \cite{chow2017algorithm}
(see \cite{bansal2017hamilton, bansal2021deepreach} for a survey).
Owing to the success of deep learning, there has also been a surge of interest in approximating high-dimensional BRTs \cite{rubies2019classification, fisac2019bridging, djeridane2006neural, niarchos2006neural, darbon2020overcoming} and optimal controllers \cite{onken2022neural} through deep neural networks (DNN).
Building upon this line of work, the authors in \cite{bansal2021deepreach} have proposed DeepReach -- a toolbox that leverages recent advances in neural implicit representations and neural PDE solvers to compute a value function and a safety controller for high-dimensional systems.
Compared to the aforementioned methods, DeepReach can handle general nonlinear dynamics, the presence of exogenous disturbances, as well as state and input constraints during the BRT computation.
However, even if we can now obtain an approximate BRT for high-dimensional systems, it is limited in usefulness by its approximate nature. 
In particular, the learned BRT can be overly optimistic; thus, no formal safety guarantees can be provided on the obtained BRT and controller. 

In this work, our goal is to compute a BRT and a safety controller for high-dimensional systems with provable safety guarantees.
Specifically, we build upon DeepReach and propose a verification method to provide safety assurances on the DeepReach solution.
The key insight of our work is to rely on the consistency between the learned value function and the implicit safety controller induced by this value function to compute a \textit{uniform correction bound} for the value function.
Essentially, if a particular state is outside the BRT (i.e., ``safe'' as per the learned value function), then starting from this state, the corresponding controller should keep the system trajectory outside the BRT at all times (i.e., the state should also be ``safe'' under the prescribed controller).
Otherwise, we cannot ensure safety of this state and the value function should be corrected such that this state is inside the BRT.
Once all such states have been added to the learned BRT, we have obtained a provably safe approximation of the true BRT.

However, in general, it is challenging to tractably compute the required value function correction across the entire state space, since the learned value function can be highly nonlinear and hard to optimize.
We propose a scenario optimization-based method for computing this correction. 
Scenario optimization is a sampling-based method to solve semi-infinite optimization problems, and has been widely used for system and control design \cite{campi2011sampling, calafiore2006scenario, campi2009scenario}.
The proposed method has general applicability and is not restricted to a specific class of
system dynamics or value functions. 
Given that the value function correction is obtained via randomization and, hence,
is a random quantity, we provide probabilistic guarantees on the safety of the recovered BRT.
However, this confidence is a design parameter and can be chosen as close to 1 as desired (within a simulation budget).

To summarize, the key contributions of this paper are:
\begin{itemize}
    \item an error correction mechanism for DeepReach that allows us to compute a provable BRT and safety controller for general dynamical systems;
    \item a practical method to compute this error correction that is not restricted to a specific class of systems; and
    \item a demonstration of the proposed approach for various dynamical systems, inspired by rocket landing and multi-vehicle collision avoidance problems.
\end{itemize}


\section*{Appendix}
\input{initial_sub/sections/appendix}

\end{document}
