text:
\begin{document}

\ifdefined\ARXIVVERSION

  \title{Collision Probabilities for Continuous-Time Systems\\Without Sampling\\
  \Large{[with Appendices]}}
\else

  \title{Collision Probabilities for Continuous-Time Systems\\Without Sampling}
\fi

\author{
  \IEEEauthorblockN{Kristoffer M.\ Frey\IEEEauthorrefmark{1}\IEEEauthorrefmark{2}, Ted J.\ Steiner\IEEEauthorrefmark{2}, and Jonathan P.\ How\IEEEauthorrefmark{1}}
  \IEEEauthorblockA{
    \IEEEauthorrefmark{1}Department of Aeronautics and Astronautics, MIT\\
    \IEEEauthorrefmark{2}The Charles Stark Draper Laboratory, Inc.\\
    Cambridge, MA 02139\\
    Email: kfrey@mit.edu
  }
}

\maketitle

\begin{abstract}
Demand for high-performance, robust, and safe autonomous systems has grown substantially in recent years.
Fulfillment of these objectives requires accurate and efficient risk estimation that can be embedded in core decision-making tasks such as motion planning.

On one hand, Monte-Carlo (MC) and other sampling-based techniques can provide accurate solutions for a wide variety of motion models but are cumbersome to apply in the context of continuous optimization.
On the other hand, ``direct'' approximations aim to compute (or upper-bound) the failure probability as a smooth function of the decision variables, and thus are widely applicable.
However, existing approaches fundamentally assume discrete-time dynamics and can perform unpredictably when applied to continuous-time systems operating in the real world, often manifesting as severe conservatism.
State-of-the-art attempts to address this within a conventional discrete-time framework require additional Gaussianity approximations that ultimately produce inconsistency of their own.
In this paper we take a fundamentally different approach, deriving a risk approximation framework directly in continuous time and producing a lightweight estimate that actually improves as the discretization is refined.
Our approximation is shown to significantly outperform state-of-the-art techniques in replicating the MC estimate while maintaining the functional and computational benefits of a direct method.
This enables robust, risk-aware, continuous motion-planning for a broad class of nonlinear, partially-observable systems.
\end{abstract}

\section{Introduction}

Robotic motion planning is a decision-making problem that must balance optimality and safety.
In the real world, these decisions are complicated by the presence of uncertainty due to imperfect sensing, partial observability, and stochastic dynamics.
This uncertainty is often difficult or impossible to explicitly bound, and safety cannot be guaranteed against all realizations of noise and disturbance.

This motivates the replacement of deterministic safety constraints with \emph{risk constraints} \cite{li2002probabilistically} that seek to compute or bound the probability of failure.
Unfortunately, exact evaluation of this probabilistic risk is challenging and computationally intractable for generic nonlinear systems.
While Monte-Carlo (MC) estimation techniques \cite{calafiore2006scenario,blackmore2010probabilistic,janson2018monte} provide a general and powerful workaround, they are still computationally-demanding and difficult to embed within a continuous motion planner.
A number of ``directly-computable'' risk approximations have been proposed \cite{li2002probabilistically,blackmore2009convex,patil2012estimating,ono2015chance}, but all fundamentally assume the system evolves in discrete time.
Many systems we care to control in practice evolve continuously, and while application of discrete-time methods is possible in these settings, the ensuing risk estimates will be highly sensitive to the chosen time discretization.
As recognized by \citet{ariu2017chance}, \citet{janson2018monte}, and others, they may be either too lax (allowing ``corner-cutting''), too conservative (leading to severe sub-optimality or artificial infeasibility), or both simultaneously.

\begin{figure}[t]
  \centering
  \includegraphics[width=0.9\linewidth]{us_v_them}
  \caption{Consider a planar system maneuvering from the green initial distribution to the red goal position in the presence of a rectangular obstacle.
           Shaded ellipses represent the state distribution at a set of discrete timesteps, $\pT = \{t_0, t_1, \ldots, t_k\}$, and reflect uncertainty under closed-loop execution.
           A common technique for computing the total collision probability leverages Boole's inequality to simply sum the violation probabilities at each step (magenta).
           However, the resulting estimate will be sensitive to the choice of $\pT$ -- too coarse, and it will underestimate; too fine, and it will ``double-count'' probability mass corresponding to trajectories that remain in collision across multiple timesteps.
           This leads to over-conservatism, artificial infeasibility, and ultimately brittle planning.
           In this paper we introduce a \emph{continuous}-time approximation (cyan) that instead aggregates probability mass \emph{entering} collision over each \emph{interval}.
           As $\pT$ is refined, our result actually improves, eventually converging to the true failure probability.

          }
  \label{fig:intro_comparison}
\end{figure}

This paper addresses this problem at its source, taking a rigorous look at the evolution of failure probability directly in continuous time.
Our approximation is general, applying to partially-observable\footnote{
          Note that under partial-observability, the process $\x(t)$ can only be considered Markov if it is ``augmented'' with the internal state of the estimator/controller.
          For the purposes of this paper we leave this implicit -- because constraints are assumed to involve only on the state of the physical plant, augmentation will contribute no analytic or computational cost.
  } stochastic systems $\x(t)$ in $\manX = \R^{n_x}$ with nonlinear It\^o dynamics
\begin{align}
  \diff \x(t) &= \vc{f}\big(t, \x(t), \u(t)\big) \diff t + \mx{G}\big(t, \x(t), \u(t) \big) \diff \w(t)  \label{eq:dynamics}
\end{align}
where $\u(t)$ and $\w(t)$ represent a vector-valued control process and Brownian motion, respectively.
Under partial observability, $\x(t)$ is not directly available to the controller, and feedback must be accomplished via a parallel observation process represented by the filtration $\Yt$ (the information available at time $t$).
The special case of fully-observable systems is easily captured here as well.
In order to simplify presentation, our main results impose some additional restrictions between the dynamics (\ref{eq:dynamics}) and the constraints that define the \emph{safe set} $\Xsafe \subset \manX$.
Even so, they accommodate a wide variety of nonlinear systems and constraints, including the ubiquitous case of second-order systems with position constraints.
Moreover, relaxation of some of these requirements will be straightforward, though beyond the scope of this paper.

The types of planning problems we address take the form
\begin{align}
  \min_{\u(\cdot)} \qquad &\E \int_0^T l\big(t, \u(t), \x(t) \big) \diff t  \label{eq:nom_cost} \\
  \text{subject to: }  & P\bigg( \bigvee_{t \in [0, T]} \x(t) \not \in \Xsafe \bigg) \leq \Delta  \label{eq:ct_constraint}
\end{align}
where the $\bigvee$ symbol is a logical OR, implying existence of a satisfying event among a (possibly uncountable) collection, and $\Delta \in (0,1)$ represents a given risk tolerance over the finite horizon $[0,T]$.
The constraint (\ref{eq:ct_constraint}) upper-bounds the probability of failure at \emph{any} time in the planning horizon.
This \emph{risk} is challenging to evaluate and optimize against because it couples states across the planning horizon as a whole.
As pointed out by \citet{ono2015chance} and others, joint constraints of this nature can be approached via Lagrangian relaxation; that is, converting the risk-\emph{constrained} problem to a risk-\emph{minimizing} one with objective
\begin{equation}\label{eq:lagrangian_opt}
  \E \int_0^T \!\!\!l\big(t, \u(t), \x(t)\big) \diff t + \lambda \Big[ P\big( \!\!\!\bigvee_{t \in [0, T]} \!\!\!\!\x(t) \not \in \Xsafe \big) - \Delta \Big]
\end{equation}
for some $\lambda \geq 0$.
However, the augmented objective in (\ref{eq:lagrangian_opt}) does not possess the time-additive \emph{Bellman} structure of (\ref{eq:nom_cost}), precluding the application of many optimal control techniques such as dynamic programming \cite{todorov2005generalized,van2012motion,ono2015chance}.
Like \cite{blackmore2009convex,ono2015chance} and other direct approximations, the approximation presented in this paper restores this time-additive structure while achieving superior accuracy for continuous-time systems.

This paper takes a fresh look at the time-evolution of the failure probability in (\ref{eq:ct_constraint}), specifically in continuous-time.
Related work is outlined in Section \ref{s:related_work}.
Section \ref{s:time_additive} leverages the language of first-exit times to produce a time-additive framework for continuous-time risk estimation.
The classic theory is rich and well-explored, but to the best of the authors' knowledge no techniques yet exist to enable computation in the context of generic nonlinear systems and nonlinear constraints.
To address this, we propose a piecewise-continuous approximation in Sections \ref{s:piecewise} and \ref{s:computation} that provably converges as the time discretization is refined -- this allows us to losslessly extend classic results for constant-coefficient systems.
As a second contribution, we identify a lightweight method to account for ``safety-thus-far'' that avoids attempting to approximate an explicit posterior and ensures conservatism without being excessively so.
Finally, we identify a large class of systems for which the requisite numerical quadratures can be computed exactly at significantly-reduced dimension, which is critical to ensure computational feasibility.
The resulting risk approximation is empirically demonstrated in Section \ref{s:experiments} to well-approximate MC estimates while retaining the computational simplicity and general applicability of direct methods, paving the way for risk-aware, continuous motion planning onboard real-world systems.

\section{Related Work}\label{s:related_work}
  Risk-aware motion planning is far from a new problem and has received much attention over the years.
  For example, when probabilistic constraints $P\big( \x(t) \not \in \Xsafe \big) \leq \Delta$ are enforced \emph{independently} for each $t$, suitable extensions of classic algorithms such as RRT \cite{luders2010chance,bry2011rapidly,van2011lqg} and differential dynamic programming (DDP) \cite{van2012motion} have been proposed.
  Alternatively, for problems with discrete time and action spaces, risk-constrained search methods such as RAO$^\star$ \cite{santana2016rao, huang2018hybrid} have shown promise.

  This paper considers problems continuous in time and action, where safety is most naturally expressed by the joint constraint (\ref{eq:ct_constraint}).
  The field of robust control provides some solutions against \emph{worst-case} (\ie, bounded) disturbances, for example \citet{majumdar2013robust} and \citet{lopez2019dynamic}.
  However, these approaches are generally restricted to fully-observable systems.
  Furthermore, their performance may be severely conservative in the average case, motivating a quantification of \emph{probabilistic} risk to allow explicit trade-off between safety and expedience during planning.

  \subsection{Sampling-Based Methods}
    MC techniques provide a general and powerful method for estimating failure probabilities, at the cost of having to run a potentially large number of simulation rollouts.
    As closed-loop partially-observable systems require simulation of plant, estimator, and controller, this can represent non-negligible amounts of computation.
    Aside from the issue of sample-complexity, which has been partially addressed by \citet{calafiore2006scenario} and \citet{janson2018monte}, sample-based estimates are discontinuous and therefore fundamentally cumbersome to incorporate into a continuous optimization framework.
    For example, \cite{calafiore2006scenario} proposes enforcement of a deterministic constraint for each sample, \cite{janson2018monte} resorts to iterative obstacle inflation in RRT, and \citet{blackmore2010probabilistic} apply Mixed-Integer techniques.
    This limitation of MC motivates the search for efficient, ``direct'' risk approximations amenable to online, \emph{continuous} optimization over motion plans.

  \subsection{Direct Risk Estimates (in Discrete Time)}
    A number of direct techniques have been developed in a discrete-time setting, where (\ref{eq:ct_constraint}) simplifies to
    \begin{equation}\label{eq:dt_constraint}
      P\big( \vee_{i = 0}^k \x(t_i) \not \in \Xsafe \big) \leq \Delta.
    \end{equation}

    Early work by \citet{li2002probabilistically} recognized that, under Linear-Time-Varying (LTV) dynamics, dispersions will be distributed as an $n_x (k+1)$-dimensioned normal distribution.
    The joint probability in (\ref{eq:dt_constraint}) can then evaluated as a high-dimensional integral via quadrature methods or sampling, both of which have complexity exponential in dimension.
    To avoid this unfavorable scaling with $k$, \citet{blackmore2009convex} use Boole's inequality (a.k.a.\ the union bound) to decompose the probability in (\ref{eq:dt_constraint}) over time as
    \begin{equation}\label{eq:dt_booles}
      P\big( \vee_{i = 0}^k \x(t_i) \not \in \Xsafe \big) \leq \sum_{i = 0}^k P( \x(t_i) \not \in \Xsafe ).
    \end{equation}
    This decoupling over timesteps is highly convenient.
    As pointed out by \citet{ono2015chance}, the right-hand side of (\ref{eq:dt_booles}) is time-additive, allowing the use of dynamic programming to minimize (\ref{eq:lagrangian_opt}).
    Alternatively, \citet{ono2008iterative} introduce the \emph{risk-allocation} formulation
    \begin{equation}\label{eq:dt_risk_allocation}
      P( \x(t_i) \not \in \Xsafe ) \leq \Delta_i \,\, \forall i \leq k \quad \text{and} \quad \sum_{i=0}^{k} \Delta_i = \Delta
    \end{equation}
    which explicitly allocates a risk ``budget'' between timesteps.
    This framework has inspired a series of other works including \cite{vitus2011feedback,ma2012fast,dai2019pcheckov}, all fundamentally dependent on decomposition (\ref{eq:dt_booles}).

    Despite its popularity, the use of Boole's inequality to decompose joint constraints over time can be severely conservative, as illustrated in Fig.~\ref{fig:intro_comparison}.
    By ignoring correlations between states at adjacent timesteps, the sum in (\ref{eq:dt_booles}) will ``double-count'' violation probabilities, particularly as the time discretization is refined.
    This problem is addressed by \citet{patil2012estimating}, who recognize that the state distributions should be \emph{conditioned} on the safety of prior timesteps.
    Because capturing this conditioning exactly is challenging, they approximate the projected and truncated distribution corresponding to each constraint as a single-dimensional Gaussian, allowing closed-form update of the state distribution parameters.
    However, Gaussianity here is inexact and the resulting estimate may not remain statistically consistent or result in a conservative (upper-bounding) risk estimate.

    When addressing continuous-time systems, a standard numeric approach to risk-aware optimization approximates the original system as discrete-time under a given time partition $\pT = \{t_0 = 0, t_1, \ldots, t_k = T\}$ and then substitutes constraint (\ref{eq:dt_constraint}) for (\ref{eq:ct_constraint}).
    This allows well-studied discrete-time techniques to be applied, but, as mentioned earlier, the results will be highly sensitive to the choice of discretization.
    Depending on the complexity and movement speed of the robot, these effects can lead to unsafe or overly-cautious behavior, or both.
    This motivates the fundamental contribution of this paper -- a direct risk approximation derived fundamentally in continuous time.

  \subsection{Stochastic Processes in Continuous-Time}
    In an effort to explicitly address the continuous-time constraint (\ref{eq:ct_constraint}) and avoid ``corner-cutting,'' \citet{ariu2017chance} apply Boole's inequality over \emph{intervals} rather than instantaneous states as
    \begin{equation}\label{eq:interval_booles}
      P\Big( \!\!\! \bigvee_{s \in [0,T]} \!\!\! \x(s) \not \in \Xsafe \Big) \leq \sum_{i=0}^{k-1} P\Big( \!\!\! \bigvee_{s \in [t_i,t_{i+1})} \!\!\!\!\! \x(s) \not \in \Xsafe \Big).
    \end{equation}
    Note that continuity allows the endpoint $\x(T)$ to be dropped.
    As with (\ref{eq:dt_booles}), this ``interval Boole's'' neglects correlations between events.
    In particular, their method assumes the state process itself follows a Brownian motion, ultimately resulting in a \emph{doubly}-conservative estimate compared to its discrete-time counterpart.
    In contrast, the approximation proposed in this paper also operates on intervals, but it addresses the nonlinear dynamics directly and avoids the use of Boole's inequality in this way, producing a much tighter risk bound.

    This paper leans heavily on the classical notion of first-passage times, a well-studied topic in the field of continuous stochastic processes.

    Indeed, the piecewise approximation presented in Section \ref{s:piecewise} of this paper is similar in spirit to that of \citet{jin2017first}, although we consider the case where the process (rather than purely the boundary) is nonlinear.

\section{How Risk Evolves in Time}\label{s:time_additive}
  Consider a \cadlag\footnote{a.k.a.\ right-continuous with left limits.}, Markov process $\z(t) \in \R^m$ with associated probability space $(\Omega, \mc{F}, P)$.
  Here, $\Omega$ represents the \emph{outcome} space, $\mc{F}$ is a \emph{sigma-algebra} over $\Omega$ (a collection of events $E \subseteq \Omega$) such that $\z(t)$ is measurable for all $t$, and $P$ is a probability measure over $\mc{F}$.
  Note that for now we do \emph{not} assume $\z(t)$ has continuous sample paths.
  We begin by exploring the ``nearly'' time-additive evolution of the \emph{exit cumulant}
  \begin{equation}
    \Fz[\z](t) \triangleq P\big(\bigvee_{s \in [0,t]} \z(s) \not \in D \big)
  \end{equation}
  for some closed set $D \subset \R^m$.

  For clarity in the following discussion, we adopt the language of \emph{passage times}, also known as \emph{hitting} or \emph{exit times} \cite{karatzsas1988,yong1999}.
  Define a parameterized family of exit times with respect to process $\z(t)$ as
  \begin{equation}
    \tc[\z]_t \triangleq \inf \big\{s \in [t, T] \,\big|\, \z(s) \not \in D\big\}
  \end{equation}
  where the infinum of the empty set is assigned to $\infty$.
  That is, $\tc[\z]_t$ refers to the first time (after $t$) that $\z(s)$ ``exits'' $D$.
  For both $\Fz[\z]$ and $\tc[\z]_t$ we will drop the explicit $\z$ specification when the corresponding process is unambiguous.

  The following useful properties can be identified.
  \begin{lemma}\label{lemma:tc_props}
    For any $0 \leq t_1 < t_2 \leq T$, the following hold:
    \begin{itemize}
      \item $\tc_{t_1} \geq t_2 \iff \z(s) \in D \quad \forall s \in [t_1, t_2)$.

      \item $\tc_0 \in [t_1, t_2) \iff \tc_0 \geq t_1, \tc_{t_1} < t_2$.

      \item For each outcome $\omega \in \Omega$, either
        \begin{enumerate}
          \item $\tc_0 = \infty$ (``no exit''), or
          \item $\z(\tc_0) \in \partial D \cup D^c$ (``smooth exit'' or a ``jump out'').
        \end{enumerate}
        where $\partial D$ and $D^c$ are the boundary and complement of $D$, respectively.
    \end{itemize}
  \end{lemma}

  It is straightforward to verify that the total exit probability can be written as $\Fz(T) = P(\tc_0 \leq T)$, and thus $\tc_0$ offers a means by which to analyze the time-evolution of $\Fz(t)$.
  In some cases $\Fz(t)$ is known to be time-differentiable \cite{jin2017first}, and thus computing this derivative (called the \emph{first-passage density}) would seem to be a natural goal.
  However, this computation is challenging for generic nonlinear processes, and instead we will settle for an interval-based ``integration'' scheme that reflects how $\Fz(T)$ can be approximated in practice.
  In later sections, through both analysis and experiment we demonstrate that this approximation indeed converges to the true value as the time discretization is refined.

  Proceeding, assume a given partition $\pT = \{t_0 = 0, t_1, t_2, \ldots, t_k = T\}$ of the fixed horizon $[0,T]$.
  A crucial advantage of the language of the \emph{first}-exit time is that it provides a natural disjointness between events, and in particular
  \begin{align}
    \Fz(T)
           &= \sum_{i=0}^{k-1} P\big(\tc_0 \in [t_i, t_{i+1})\big) + P(\tc_0 = T)  \label{eq:ival_cond1}.
  \end{align}
  Note that, in contrast to (\ref{eq:interval_booles}), the relation (\ref{eq:ival_cond1}) holds with equality.
  Proceeding from here and adopting the shorthand $\z_i \triangleq \z(t_i)$, $\tc_i \triangleq \tc_{t_i}$, and so on, the above can be written
  \begin{align}
    &\sum_{i=0}^{k-1} \Big( P\big(\tc_0 \in [t_i, t_{i+1}), \z_i \in D\big) + P(\tc_0 = t_i, \z_i \not \in D) \Big) \nonumber \\
      &\qquad \qquad  + P(\tc_0 = T, \z(T) \in D^c)  \label{eq:pt:marginalize_endpoint}  \\
    &= \sum_{i=0}^{k-1} P\big(\tc_0 \in [t_i, t_{i+1}), \z_i \in D \big) + \sum_{i=0}^k P(\tc_0 = t_i, \z_i \in D^c)  \label{eq:Fz_sum0}
  \end{align}
  where in (\ref{eq:pt:marginalize_endpoint}) we split the probability over the event that $\z_i \in D$.
  Note that the right-hand summation in (\ref{eq:Fz_sum0}) involves the probabilities that discontinuous sample paths ``jump out'' of $D$ at each of the partition points $t_i$.
  Examining the first set of terms, Lemma \ref{lemma:tc_props} allows
  \begin{align}
    &P\big(\tc_0 \in [t_i, t_{i+1}), \z_i \in D \big)    \label{eq:Phi_deriv} \\

    &= P\big(\tc_i < t_{i+1} \cbar \tc_0 \geq t_i, \z_i \in D \big) P(\tc_0 \geq t_i, \z_i \in D). \nonumber
  \end{align}
  The conditioning on $\tc_0 \geq t_i$ and $\z_i \in D$ in (\ref{eq:Phi_deriv}) implies that the process has not exited ``yet,'' imposing a specific posterior over $\z_i$ that we will call the \emph{anthropic} distribution\footnote{
    In cosmology, the \emph{anthropic principle} remarks that life can only observe universes that themselves allow for the existence of life.
  }
  \begin{equation}\label{eq:pibar}
    \pibar_t(\diff \z) \triangleq P(\z(t) \in \diff \z \cbar \tc_0 \geq t, \z(t) \in D).
  \end{equation}

  Before proceeding further, define the function
  \begin{equation}\label{eq:Phiz_def}
    \Phiz[\z](t_i, t_{i+1}; \mu_i) \triangleq P(\tc_i < t_{i+1} \cbar \z_i \sim \mu_i)
  \end{equation}
  for any (not necessarily normalized) measure $\mu_i$ over $\R^m$.
  (\ref{eq:Phiz_def}) captures the probability of an exit (not necessarily the first) in the interval $[t_i, t_{i+1})$, given that $\z_i$ is distributed according to $\mu_i$.
  As before, the $\z$ subscript will be left implicit where possible.
  Because the process $\z(t)$ is assumed Markov and applying (\ref{eq:Phi_deriv}), sum (\ref{eq:Fz_sum0}) can be re-written in terms of $\Phiz$ and $\pibar_t$ as
  \begin{align}
    \Fz(T) &= \sum_{i=0}^{k-1} \Phiz(t_i, t_{i+1}; \pibar_i) P( \tc_0 \geq t_i, \z_i \in D) \notag \\
           &~~+ \sum_{i=0}^k P(\tc_0 = t_i, \z_i \in D^c). \label{eq:Fz_sum2}
  \end{align}

  A key challenge in evaluation of (\ref{eq:Fz_sum2}) is computation of $\pibar_i$.
  For one thing, its support is clearly limited to $D$, which is sufficient to ensure non-Gaussianity.
  This motivates us to avoid attempting to approximate or bound $\pibar_t$ directly and instead decompose it via Baye's rule, producing
  \begin{align}
    \pibar_t( \diff \z )
                  &= \frac{ P\big( \tc_0 \geq t, \z(t) \in D \cbar \z(t) = \z \big) P\big(\z(t) \in \diff \z \big) }{P( \tc_0 \geq t, \z(t) \in D )}  \nonumber \\
                  &\triangleq \frac{ \Psibar_t(\z) \bbar_t(\diff \z) }{ P(\tc_0 \geq t, \z(t) \in D). }  \label{eq:pi_decomp}
  \end{align}
  where $\bbar_t$ is the \emph{a priori} distribution of $\z(t)$, and we refer to $\Psibar_t : \R^m \mapsto [0, 1]$ as the \emph{anthropic likelihood}.
  It is straightforward to show
  \begin{equation}\label{eq:Psi_identity}
    \Phiz(t_i, t_{i+1}; \pibar_i) P( \tc_0 \geq t_i, \z_i \in D ) = \Phiz(t_i, t_{i+1}; \Psibar_i \bbar_i),
  \end{equation}
  and therefore (\ref{eq:Fz_sum2}) can be written
  \begin{equation}\label{eq:Fz_sum_apriori}
    \Fz(T) = \sum_{i=0}^{k-1} \Phiz(t_i, t_{i+1}; \Psibar_i \bbar_i) + \sum_{i=0}^k P(\tc_0 = t_i, \z_i \in D^c).
  \end{equation}
  As discussed in Section \ref{ss:anthro_approx}, and in contrast to the case of $\pibar_t$, identifying conservative approximations for $\Psibar_t$ will be both straightforward and effective.

  \subsection{Analogous Development Under Filtration $\mc{F}_t$}
    Though not the primary focus of this paper, we note that the preceding development can be applied analogously in the context on a filtration $\mc{F}_t$ representing information that becomes available \emph{during} execution as opposed to \emph{a priori}.
    For example, in the context of $\x(t)$ as defined in (\ref{eq:dynamics}), we might consider the observation filtration $\Ft = \Yt$.

    Though the information generated by $\Ft$ itself evolves randomly, modification of the above discussion is straightforward and we simply provide some analogous definitions here for clarity.
    As in (\ref{eq:pibar}), we can define the \emph{anthropic belief}
    $  \pi_t(\diff \z) \triangleq P\big( \z(t) \in \diff \z \cbar \tc_0 \geq t, \z(t) \in D, \Ft \big)$,
    and like (\ref{eq:pibar}) it has the structure
    \begin{equation}\label{eq:anthro_belief}
      \pi_t(\diff \z) = \frac{\Psi_t(\z) b_t(\diff \z)}{ P\big(\tc_0 \geq t, \z(t) \in D \cbar \Ft \big) }
    \end{equation}
    where the familiar estimation belief $b_t$ and anthropic likelihood $\Psi_t$ are adapted to $\Ft$.
    An analogous identity to (\ref{eq:Psi_identity}) can be established, yielding
    \begin{equation}\label{eq:Fz_sum}
      \Fz(T) = \E \Big[ \sum_{i = 0}^{k-1} \Phiz(t_i, t_{i+1}; \Psi_i b_i) \Big] + \sum_{i=0}^k P(\tc_0 = t_i, \z_i \in D^c).
    \end{equation}
    where the expectation is taken over $\Ft$.

\section{A Piecewise-Continuous Approximation}\label{s:piecewise}
  The preceding section introduced a framework for computing the failure probability based on a helper function $\Phiz$ defined in (\ref{eq:Phiz_def}).
  Rather than attempting to evaluate $\Phiz[\x]$ directly for the nonlinear process $\x(t)$, this section introduces an approximating process \emph{in the constraint space}, for which computation is made tractable.
  This simplified process is piecewise-continuous according to the given time discretization $\pT$, which controls the ``accuracy'' of the approximation.
  Though $\pT$ in practice will likely be dictated by the computational resources available, we prove that as it is refined the corresponding failure probability estimate converges to $\Fz[\x](T)$.
  That is to say, the approximation is asymptotically lossless.

  Assume the feasible set $\Xsafe \subset \manX$ is defined by the sublevel sets $g_j(\x) \leq 0$ for a set of $m$ twice-differentiable, real-valued functions $\{g_j\}$.
  Let $\vc{g}(\x)$ refer to the stacked vector in $\R^m$, and let $\vc{a}_j(\x)$ and $\mx{H}_j(\x)$ refer to the gradient vector and Hessian of each $g_j$, respectively.
  Note that $\x \in \Xsafe$ is equivalent to the statement that $g_j(\x) \leq 0$ for all $j$.
  This naturally motivates a consideration of $\y(t) \triangleq \vc{g}\big( \x(t) \big)$ as a process with respect to the non-positive orthant $\setOm \subset \R^m$.

  Before we proceed, however, computation and analysis will require some regularity conditions on the control process $\u(t)$, which evolves according to $\Yt$ and some (potentially non-deterministic) policy.
  Rather than complicate the discussion by attempting to account for all such possibilities, we make the following simplifying assumption.
  \begin{assumption}[$\y(t)$ Locally Independent of Control]\label{ass:loc_ind_u}
    We require that either
    \begin{enumerate}
      \item we have state-feedback $\u(t) = \kappa\big(t, \x(t)\big)$, with $\kappa$ deterministic and Lipschitz, or
      \item $\vc{a}_j^\T(\x) \vc{f}(t, \x, \u)$ and $\vc{a}_j^\T(\x) \mx{G}(t, \x, \u)$ are independent of $\u$ for all constraints $j$.
    \end{enumerate}
  \end{assumption}
  Assumption \ref{ass:loc_ind_u} ensures that the dynamics of process $\y(t)$ depend on the control only via the state $\x(t)$, and notably accommodates second-order systems under workspace constraints (\ie, obstacle avoidance).
  Partial relaxation of this assumption will in many cases require only minor modifications of the analysis and computations presented in this paper.

  Proceeding, we consider a piecewise approximation of this process constructed from the time discretization $\pT$.
  Using the shortcuts $\vc{f}_t \triangleq \vc{f}\big(t, \x(t), \u(t)\big)$, $\mx{G}_t \triangleq \mx{G}\big(t, \x(t), \u(t)\big)$, and so on, define the process $\yhat(s)$ such that
  \begin{align}
    \yhat(s) &\triangleq \vc{g}\big( \x(t) \big) + \int_t^s \vc{h}_t \diff \sprime + \int_t^s\!\!\! \mx{\Sigma}_t \diff \vc{w}(\sprime)  \\
    h_{t,j} &\triangleq \vc{a}_j\big(\x(t)\big)^\T \vc{f}_t + \frac{1}{2} \tr \Big( \mx{G}_t^\T \mx{H}_j\big(\x(t)\big) \mx{G}_t \Big)   \label{eq:h_def} \\
    \mx{\Sigma}_t &\triangleq \begin{bmatrix} \vc{a}_1^\T\big(\x(t)\big) \\ \vdots \\ \vc{a}_m^\T\big(\x(t)\big) \end{bmatrix} \mx{G}_t
                  \triangleq \begin{bmatrix} \vc{\sigma}_1\big(t,\x(t)\big) \\ \vdots \\ \vc{\sigma}_m\big(t,\x(t)\big) \big) \end{bmatrix}
  \end{align}
  where $t = \max \{ t_i \in \pT \cbar t_i \leq s \}$ corresponds to the start point of the interval $[t_i, t_{i+1})$ containing $s$.
  Here $\vc{h}_s$ and $\mx{\Sigma}_s$ can be identified as the drift and diffusion coefficients of $\y(s)$ as prescribed by It\^o's formula \cite[Thm.~5.5]{yong1999}.
  Thus, the piecewise-continuous $\yhat(s)$ has constant coefficients over each interval and coincides with $\y(s)$ for every $s = t \in \pT$.

  Approximating $\Fz[\yhat]$ will be much easier than working directly with $\Fz[\x]$, and in the following we show that, under mild conditions, the substitution of $\yhat(s)$ for $\y(s)$ can be considered asymptotically lossless.

  \begin{figure}[t]

    \centering
    \begin{tikzpicture}[scale=1, every node/.style={transform shape}]
      \node[anchor=south west,inner sep=0] (image) at (0,0)
          {\includegraphics[width=0.97\linewidth,height=2.5cm]{yhat}};

      \node at (0.95, 2.4)   {\footnotesize $y_j(t)$};
      \node at (0.95, 1.97) {\footnotesize $\yjhat(t)$};

      \node at (-0.1, 1.5) {\footnotesize $0$};

      \node at (0.05, 1.15) {\footnotesize $0$};
      \node at (8.58, 1.15) {\footnotesize $T$};

    \end{tikzpicture}
    \caption{
        An illustration of the continuous process $y_j(t) = g_j(\x(t))$ and our approximation $\yjhat(t)$.
        The two processes coincide at each discrete timestep $t_i \in \pT$, and as $\pT$ is refined, Prop.~\ref{prop:yhatn_puc} establishes that the two converge pathwise-uniformly.
        Zero-crossings of $y_j(t)$ imply constraint violations, and therefore exit-times of $\yhat(t)$ from the non-positive orthant approximate those of $\x(t)$ from $\Xsafe$.
    }
    \label{fig:yhat}
  \end{figure}

  \subsection{Convergence Analysis}
    Consider evaluating $\Fz[\yhat](T)$ for the process $\yhat(t)$ -- the result will, of course, depend on $\yhat(t)$ and therefore the underlying time discretization $\pT$.
    We show that as $\pT$ is refined, $\yhat(t)$ converges to $\y(t)$ in a specific sense, and that $\Fz[\yhat](T)$ converges to $\Fz[\y](T) = \Fz[\x](T)$.
    Proofs are left to \refAppProofs.

    \begin{definition}[Pathwise-Uniform Convergence]\label{defn:puc}
      For a sequence of processes $\zn(t)$ and an associated real-valued sequence $\eps_n \rightarrow 0$, we say $\zn(t)$ converges $\eps_n$-\emph{pathwise-uniformly} to process $\z(t)$ if there exists for each $\omega \in \Omega$ an $0 < M_\omega < \infty$ almost surely such that
      \begin{equation*}
        P \Big( \sup_{s \in [0,T]} \norm{\zn(s) - \z(s)} > \eps_n M_\omega \Big) \rightarrow 0.
    \end{equation*}
    \end{definition}

    Note that pathwise-uniform convergence (PUC) does not imply that $\zn(\cdot; \omega)$ converges to $\z(\cdot; \omega)$ for any fixed outcome $\omega \in \Omega$.
    Additionally, it is specific to the given sequence $\eps_n$.
    Nonetheless, we will see that this condition is sufficient to ensure a non-trivial convergence-of-probabilities result.

    \begin{definition}[Supremum Process]\label{defn:sup_proc}
      For a process $z(t)$, the corresponding \emph{supremum process} is
      $
        z^\star(t) \triangleq \sup_{s \in [0,t]} z(s).
      $
    \end{definition}
    \begin{theorem}\label{thm:puc_implies_P_conv}
      If a given sequence of $\R^m$-valued processes $\zn(t)$ converges $\eps_n$-pathwise-uniformly to $\z(t)$, and we have that $P( z_j^\star(T) = 0 ) = 0$ for all $j$, then
      \begin{equation*}
        P\Big(\bigvee_{s \in [0,T]} \zn(s) \not \in \setOm \Big) \rightarrow P \Big( \bigvee_{s \in [0,T]} \z(s) \not \in \setOm \Big).
      \end{equation*}

    \end{theorem}

    For technical reasons, we will be required to make the following additional assumption relating the diffusion matrix $\mx{G}$ and the constraint gradients $\{\vc{a}_j\}$.
    \begin{assumption}[Paired Lipschitz:]\label{ass:paired_Lipschitz}
      With respect to a given function $g : \manX \mapsto \R$ and diffusion matrix $\mx{G}(t, \x, \u)$, either
      \begin{enumerate}
        \item $\mx{G}$ is constant and the function gradient, $\vc{a}$, is Lipschitz,
        \item $\vc{a}$ is constant (\ie, $g$ is linear) and $\mx{G}$ is Lipschitz, or
        \item both $\vc{a}$ and $\mx{G}$ are bounded and Lipschitz.
      \end{enumerate}
    \end{assumption}

    \begin{remark}
      The requirements that Assumption \ref{ass:paired_Lipschitz} imposes on $\mx{G}$ in no way imply that disturbances are bounded.
      Furthermore, when constraints are linear $\mx{G}$ is only required to be Lipschitz, so cases of multiplicative noise can be handled here as well.
    \end{remark}

    This leads to our main result.
    \begin{proposition}\label{prop:yhatn_puc}
      Let Assumptions \ref{ass:loc_ind_u} and \ref{ass:paired_Lipschitz} hold for the $m$ constraint functions $\{g_j\}$.
      Additionally, assume:
      \begin{enumerate}
        \item $\vc{f}_t$ is pathwise-bounded \as, and $\E \int_0^T \norm{\vc{f}_t}^2 \diff t < \infty$.
        \item $\E \norm{ \mx{G}_t }^2$ bounded over $[0,T]$, and $\E \int_0^T \norm{\mx{G}_t}^2 \diff t < \infty$.
      \end{enumerate}
      Then, for a sequence of partitions $\pT_n = \{0 = t_0, t_1, \ldots, t_{k_n} = T\}$ over the compact interval $[0, T]$ with mesh $\delta_n \rightarrow 0$ and $k_n \delta_n \leq c T$ for fixed $c > 0$,
      there exists a sequence $\eps_n \rightarrow 0$ such that the sequence of processes $\yhatn(s)$ converges $\eps_n$-pathwise-uniformly to $\y(s)$.
    \end{proposition}

    From Thm.~\ref{thm:puc_implies_P_conv} and the definition of $\y(t)$ the following corollary immediately follows.
    \begin{corollary}\label{cor:Fz_convergence}
      If the assumptions of Prop.~\ref{prop:yhatn_puc} hold and $P\big( y^\star(T) = 0 \big) = 0$ for all $j$
      then $\Fz[\yhatn](T) \rightarrow \Fz[\x](T)$.
    \end{corollary}

    \begin{remark}

      The condition $P\big( y_j^\star(T) = 0 \big) = 0$ is somewhat difficult to verify; while conditions over which the variable $y_j(t)$ admits a density have been somewhat examined, the case of the supremum process is less explored.
      For the purposes of the paper we consider this, as well as the regularity conditions required by Prop.~\ref{prop:yhatn_puc}, to be fairly mild.
    \end{remark}

    Corollary \ref{cor:Fz_convergence} makes clear that the substitution of $\yhat(s)$ for $\y(s)$ is (asymptotically) loss-less for the purposes of computing the risk probability $\Fz[x](T)$.
    It is worth noting that this result does not anywhere assume that the dynamics or observation model are Gaussian or LTV, although practical computation of $b_t$ may still require this assumption.

\section{Computation}\label{s:computation}
  The results of Section \ref{s:piecewise} indicate that we can perform computation with the simplified process $\yhat(t)$ and the result will approximate the true $\Fz[\x](T)$ in an asymptotic sense.
  While $\yhat(t)$ is only piecewise-continuous and thus ``jump out'' probabilities in (\ref{eq:Fz_sum}) may be non-zero, our empirical results indicate that they can be safely ignored.
  Identifying a rigorous proof is left as a topic for future work.
  Additionally, for the sake of presentation we make the simplifying assumption that $P(\x_0 \not \in \Xsafe) = 0$.
  In light of these considerations, we focus on computing
  \begin{equation}\label{eq:Fyhat}
    \Fhat(T) \triangleq \sum_{i = 0}^{k-1} \Phiz[\yhat](t_i, t_{i+1}; \Psibar_i \bbar_i).
  \end{equation}

  In order to compute (\ref{eq:Fyhat}), we require both a means of computing the quantity $\Phiz[\yhat](t_i, t_{i+1}; \mu_i)$ for arbitrary measure $\mu_i$ and a means of approximating $\Psibar_i \bbar_i$.
  The former is the subject of \ref{ss:phi_compute} and the latter of \ref{ss:anthro_approx}.

  \subsection{Bounding $\Phiz[\yhat]$}\label{ss:phi_compute}
    The constant-coefficient approximation $\yhat(t)$ allows us to leverage classic results in the study of first-exit times.
    First, recalling (\ref{eq:Phiz_def}) we apply the union bound to decouple individual constraints, producing a sum over element-wise diffusions $\yjhat$
    \begin{align}
      &\Phiz[\yhat](t_i, t_{i+1}; \mu_i)

                                        \leq \sum_{j=1}^m \Phiz[\yjhat](t_i, t_{i+1}; \mu_i).    \label{eq:Phi_yhat}
    \end{align}

    \xx{Distinguish from dt booles}

    Letting $\dXsafe$ refer to the boundary of $\Xsafe$, a classic result (see \citet[Ch.~3, Eq.~5.12]{karatzsas1988}) provides the following lemma.
    \begin{lemma}\label{lemma:psi}
      Under Assumption \ref{ass:loc_ind_u}, for any (not-necessarily normalized) measure $\mu_i$ over $\Xsafe$, and such that $\mu_i(\dXsafe) = 0$,
      \begin{align}
        &\Phiz{\yjhat}(t_i, t_{i+1}; \mu_i) = P\Big( \sup_{s \in [t_i, t_{i+1})} \!\!\! \yjhat(s) > 0  \cBar \yjhat(t_i) \sim g_j( \mu_i ) \Big)  \nonumber \\
        &\, = \int_{\Xsafe}\!\!\! \psi\Big(g_j(\x), h_j(t_i,\x), \norm{\vc{\sigma}_j(t_i,\x)}, \Delta_i\Big) \diff \mu_i(\x)  \label{eq:Phi_zj}
      \end{align}
      where $\Delta_i = t_{i+1} - t_i$ and we define $\psi(z, h, \sigma, \Dt) \triangleq$
      \begin{equation}\label{eq:psi}
      \quad 1 - \phi\Big(\frac{-h\Dt - z}{\sigma \sqrt{\Dt} } \Big) + \exp \Big\{-\frac{2 h z}{\sigma^2} \Big\} \phi\Big(\frac{-h\Dt + z}{\sigma \sqrt{\Dt} } \Big)
      \end{equation}
      and $\phi$ refers to the CDF of the standard normal distribution.
    \end{lemma}

    Note that the condition $\mu_i(\dXsafe) = 0$ is satisfied for all measures which admit a density (a mild assumption).

    Like the instantaneous probability $P\big( x(t_i) \not \in \Xsafe \big)$, computing (\ref{eq:Phi_zj}) in general requires evaluating a multi-dimensional integral over the state space $\manX = \R^{n_x}$.
    Numerical quadrature techniques have exponential complexity in dimension (the curse of dimensionality).
    Fortunately, depending on the structure of the dynamics and constraints, significant dimensionality reductions can often be found.
    In particular, consider Gaussian-distributed second-order systems with workspace constraints (such as obstacle avoidance).
    If system noise is injected only in the velocity and acceleration dynamics (corresponding to $\vc{\sigma}_j = 0$ for all $j$), then the integral in (\ref{eq:Phi_zj}) reduces to workspace dimension $n_p$ (commonly 2 or 3).
    This ensures that $\Phiz[\yhat]$ can be evaluated efficiently for a large class of dynamics.
    Details are provided in \refAppRedQuad.

  \subsection{Approximating Anthropic Belief}\label{ss:anthro_approx}
    Given we can compute (\ref{eq:Phi_yhat}), we still require some means of estimating the anthropic state distribution via $\pibar_t$ or $\Psibar_t \bbar_t$.
    Though inexact, Gaussian approximations for $\bbar_i$ are applicable and widely-used whenever the state dynamics are sufficiently smooth and, as can be seen in our experimental results, accurately model the dispersion of even nonlinear MC trajectories.
    However, modeling the anthropic information represented by $\Psibar_i$ or $\pibar_i$ is not so straightforward.

    \subsubsection{Approximating $\pibar_i$ as a Gaussian}\label{sss:Gaussian_pi}
      In existing work, \citet{patil2012estimating} propose a Gaussian approximation for the anthropic distribution $\pibar_i$ (or analogously, anthropic belief $\pi_i$) at each timestep $t_i$.
      This is attractive because it implies $\pibar_i$ can be propagated in parallel with $\bbar_i$.
      Likewise, we can consider approximating $\pibar_i$ by conditioning on safety (via a recursive Gaussian approximation) at each prior timestep
      \begin{equation}\label{eq:pibar_gauss}
        \hat{\pibar}_i(\diff \x) \approx P\big(\x_i \in \diff \x \cbar \x_l \in \Xsafe \, \forall l \leq i \big).
      \end{equation}

      We refer the interested reader to their paper for implementation details.

      Again considering (\ref{eq:Fyhat}), we can apply identity (\ref{eq:Psi_identity}) and directly produce an upper-bound in terms of $\pibar_i$
      \begin{equation}\label{eq:ival_gauss}
        \Fhat(T) \leq \sum_{i=0}^{k-1} \Phiz[\yhat](t_i, t_{i+1}; \pibar_i).
      \end{equation}
      Together with our Gaussian approximation $\hat{\pibar}_i$ from (\ref{eq:pibar_gauss}), we refer to the resulting estimate as \texttt{ival\_gauss}.

      While the independence and Gaussianity assumptions involved here are convenient, they are also heuristic.
      This can lead to inconsistent $\pibar_t$ estimates and ultimately unpredictable risk estimation, as will be seen in Section \ref{s:experiments}.

    \subsubsection{A Simple Alternative}\label{sss:ival_safe}
      Rather than attempting to estimate $\pibar_t$, we propose a simple, yet surprisingly powerful alternative based on the anthropic likelihood $\Psibar_t$.
      A clear upper-bound follows directly from the definition as
      \begin{equation}\label{eq:Psi_dumb}
        \Psibar_t(\x) = P\big(\tc[\x]_0 \geq t, \x(t) \in \Xsafe \cbar \x(t) = \x \big) \leq \Ind_{\Xsafe}(\x).
      \end{equation}
      Applying (\ref{eq:Psi_dumb}) to (\ref{eq:Fyhat}) yields a conservative estimate
      \begin{equation}\label{eq:ival_safe}
        \Fhat(T) \leq \sum_{i = 0}^{k-1} \Phiz[\yhat](t_i, t_{i+1}; \Ind_{\Xsafe} \bbar_i)
      \end{equation}
      which accumulates ``new'' exits over each interval -- note the contrast with the naive ``interval Boole's'' approximation given by (\ref{eq:interval_booles}).
      Because it restricts exit flow to probability mass which is ``safe'' at the start of each interval, we refer to approximation (\ref{eq:ival_safe}) as \texttt{ival\_safe}.

      In principle, (\ref{eq:ival_safe}) shares one of the same weaknesses of the discrete-time Boole's approximation in (\ref{eq:dt_booles}) -- it is not bounded above, and indeed may diverge as $\pT$ is refined.
      However, it avoids ``locally'' double-counting by counting only probability mass \emph{leaving} $\Xsafe$ over each interval.
      This interpretation is illustrated in Fig.~\ref{fig:intro_comparison}, and our empirical results suggest \texttt{ival\_safe} does very well in practice.
      Furthermore, compared to \cite{patil2012estimating} or (\ref{eq:ival_gauss}) it requires no ``extra'' belief propagation at all, making this approximation particularly lightweight.

  \begin{figure}[t]
    \centering
    \begin{tikzpicture}[scale=1.0, every node/.style={transform shape}]
      \node[anchor=south west,inner sep=0] (image) at (0,0)
          {\includegraphics[width=0.45\linewidth,trim={0.5cm, 0.4cm, 0.4cm, 0.37},clip]{cc_demo1}};
      \node at (2.15, 6) {\footnotesize $\pT_1 \sim 10$ Hz};
      \node[rotate=90] at (-0.2, 1.7) {\scriptsize Estimated $\Fz[\x](t)$};
    \end{tikzpicture}
    \hfill
    \begin{tikzpicture}[scale=1.0, every node/.style={transform shape}]
      \node[anchor=south west,inner sep=0] (image) at (0,0)
          {\includegraphics[width=0.45\linewidth,trim={0.5cm, 0.4cm, 0.4cm, 0.37},clip]{cc_demo3}};
      \node at (2.15, 6) {\footnotesize $\pT_2 \sim 60$ Hz};
    \end{tikzpicture}
    \vspace{-0.25cm}
    \caption{
      We simulate a noisy Dubin's car system navigating in a narrow passageway, computing collision probabilities under coarse (left) and fine (right) time discretizations, $\pT_1$ and $\pT_2$.
      MC trajectories (red) simulate a fixed 60Hz feedback control rate, and a corresponding \emph{a priori} state distribution $\bbar_i$ (blue) is propagated via an LTV-Gaussian assumption.
      The anthropic estimates $\{\hat{\pibar}_i\}_j$ (yellow) are computed according to each $\pT_j$.
      The underlying plots confirm the divergence of \texttt{dt\_booles} under fine discretizations, and illustrate the susceptibility of \texttt{dt\_gauss} and \texttt{ival\_gauss} to inconsistency in the Gaussian $\hat{\pibar}_i$ estimates.
      In contrast, our method \texttt{ival\_safe} avoids estimating $\pibar$ entirely, producing a lightweight and consistent risk estimate.
    }
    \label{fig:cc_demo}
  \end{figure}

\section{Results}\label{s:experiments}
  \begin{figure*}[t]
    \centering

    \begin{tikzpicture}[scale=1.0, every node/.style={transform shape}]
      \node[anchor=south west,inner sep=0] (image) at (0,0)
          {\includegraphics[width=0.225\linewidth,trim={0.5cm, 0.5cm, 0.5cm, 0.45cm},clip]{open51_62}};
      \node[rotate=90] at (-0.1, 1.7) {\scriptsize Estimated $\Fz[\x](t)$};
    \end{tikzpicture}
    \begin{tikzpicture}[scale=1.0, every node/.style={transform shape}]
      \node[anchor=south west,inner sep=0] (image) at (0,0)
          {\includegraphics[width=0.225\linewidth,trim={0.5cm, 0.5cm, 0.5cm, 0.45cm},clip]{open51_78}};
    \end{tikzpicture}

    \hfill
    \vrule
    \hfill

    \begin{tikzpicture}[scale=1.0, every node/.style={transform shape}]
      \node[anchor=south west,inner sep=0] (image) at (0,0)
          {\includegraphics[width=0.225\linewidth,trim={0.5cm, 0.5cm, 0.5cm, 0.45cm},clip]{closed17_04}};
      \node[rotate=90] at (-0.2, 1.7) {\scriptsize Estimated $\Fz[\x](t)$};
    \end{tikzpicture}
    \begin{tikzpicture}[scale=1.0, every node/.style={transform shape}]
      \node[anchor=south west,inner sep=0] (image) at (0,0)
          {\includegraphics[width=0.225\linewidth,trim={0.5cm, 0.5cm, 0.5cm, 0.45cm},clip]{closed17_01}};
    \end{tikzpicture}
    \vspace{-0.2cm}
    \caption{
      Sample Dubin's car trajectories in randomly-generated environments over a $T = 2.5$ [s] horizon under 60 Hz LQG feedback control towards the (green) goal.
      The two scenarios on the left demonstrate \emph{nominally}-safe plans, which can experience significant risk upon stochastic execution.
      On the right are risk-constrained plans (generated using our \texttt{ival\_safe} approximation with $\Delta = 10\%$).
      MC trajectories are shown in red, and dispersions are well-approximated by the \emph{a priori} distributions $\{\bbar_i\}$ shown in blue.
      The Gaussian anthropic belief estimates $\{\hat{\pibar}_i\}$ are shown in yellow.
      The lower plots show $\Fz[\x](t)$ estimates as compared to MC (estimated standard error is indicated by the shaded region).
      Note that \texttt{dt\_booles} is prone to extreme over-estimation, while the two $\{\hat{\pibar}_i\}$-dependent methods \texttt{dt\_gauss} and \texttt{ival\_gauss} suffer from the inconsistency of the underlying Gaussianity assumption.
      In contrast, our proposed method \texttt{ival\_safe} achieves the best approximation performance and avoids explicit $\pibar_i$ approximation.
    }
    \label{fig:eval_samples}
  \end{figure*}

  \begin{table*}[h]
    \centering
    \caption{Compute times for Dubin's car MATLAB implementation on consumer laptop. QD refers to quadrature dimension.}
    \vspace{-0.3cm}
    \label{table:compute_times}
    \begin{tabular}{c  @{\quad\qquad\qquad}  c @{\qquad} c  @{\quad\qquad\qquad}  c @{\qquad} c}
      \toprule
      \vspace{0.07cm}
      Monte Carlo     &                          &                              &  $\sum_{i=0}^k P\big(x(t_i) \not \in \Xsafe\big)$  & $\sum_{i=0}^{k-1} \Phiz(t_i, t_{i+1}; \Ind_\Xsafe \bbar_i)$  \\
      ($N = 1000$)    &  LTV-Gaussian $\bbar_i$  &   Gaussian $\hat{\pibar}_i$  &  QD: $(n_p-1) = 1$                            & QD: $n_p = 2$ \\
      \midrule
         39.1 [s]     &                0.034 [s] &                  0.054 [s]   &                                      0.034 [s]  &                   0.316 [s]   \\
      \bottomrule
    \end{tabular}
    \vspace{-0.3cm}
  \end{table*}

  To evaluate the accuracy of our approximation, we simulate a second-order Dubin's car system in the plane.
  System details are provided in \refAppSys.
  Noise in dynamics and observations render the system partially-observable, and an LQG-style feedback control policy stabilizes around a nominal trajectory.
  We evaluate four principle approximations.
  Discrete-time methods \texttt{dt\_booles} and \texttt{dt\_gauss} refer to the naive Boole's approximation (\ref{eq:dt_booles}) and the conditioned variant proposed by \cite{patil2012estimating}, respectively.
  Our methods \texttt{ival\_gauss} and \texttt{ival\_safe} are described in the preceding section.

  First, we successfully use our \texttt{ival\_safe} estimate to optimize a risk-constrained ($\Delta = 2 \%$) plan in a tight environment, shown in Fig.~\ref{fig:cc_demo}.
  We then computed risk estimates for this plan under multiple time discretizations.
  As predicted, \texttt{dt\_booles} diverges severely as $\pT$ is refined while \texttt{ival\_safe} converges correctly.
  Interestingly, \texttt{dt\_gauss} and \texttt{ival\_gauss} converge to incorrect values, likely due to inconsistency in the underlying Gaussian $\hat{\pibar}_t$ estimate.

  Next, we perform a larger statistical evaluation over random, programmatically-generated environments.
  We consider two different regimes: \emph{nominally}-safe trajectories, which ensure only that the nominal trajectory is collision-free, and \emph{risk}-constrained trajectories optimized using \texttt{ival\_safe} with $\Delta = 10\%$.
  This segregation allows us to identify potentially distinct statistical performance in scenarios representing the full range of risk values $[0,1]$ and over the smaller range of ``planning-relevant'' values in $[0, \Delta]$.
  These environments are non-convex, and the locally-optimal trajectories were generated via MATLAB's \texttt{fmincon} routine \cite{MatlabOTB}.
  To ensure that the batch of experiments is not dominated by ``uninteresting'' cases, we specifically reject scenarios where the \emph{unconstrained} optimal has negligible failure probability (\ie, when the obstacles are not relevant).
  Computation times for our MATLAB implementation of various methods are shown in Table \ref{table:compute_times} -- as predicted, our method is significantly faster than MC while representing some additional complexity compared to discrete-time methods.
  A more detailed discussion is presented in \refAppRedQuad.

  Some representative scenarios are shown in Fig.~\ref{fig:eval_samples}, and statistics are reported in Table \ref{table:batch_eval} over 100 nominally-safe and 50 risk-constrained scenarios.
  For reference, the average MC-evaluated risk in each of the two test batches is listed in the top of the column.
  The \textbf{Bias} column lists the mean (signed) difference between the estimate and a $1000$-sample MC estimate, and $P(\text{\textbf{Conservative}})$ reports the percentage of cases where the estimate was greater than (or within 5\%) of the MC ``truth.''
  The proposed method \texttt{ival\_safe} outperforms or matches all others in bias, root-mean-squared error (\textbf{RMSE}), and median relative error (\textbf{MRE}), while remaining conservative (\ie, safe) in a significant majority of trials.
  Both \texttt{dt\_gauss} and \texttt{ival\_gauss} outperform \texttt{dt\_booles} but suffer from inconsistency of the Gaussian $\hat{\pibar}_i$ estimate.

  \begin{table}[t]
    \centering
    \caption{$\Fz[\x](T)$ estimation statistics for Dubin's car system.}
    \label{table:batch_eval}
    \vspace{-0.3cm}
    \begin{tabular}{l@{\qquad}r@{\quad}r@{\quad}r@{\quad}r}
      \toprule
      \textbf{Batch / Method}         & \textbf{Bias} & \textbf{RMSE} & \textbf{MRE} & $P($\textbf{Conservative}$)$ \\
      \midrule

      Nominally-safe          & MC : 0.2649  & \\
      \texttt{dt\_booles}     &       +2.9780   &  4.908       &    948 \%        & \textbf{100} \%   \\
      \texttt{dt\_gauss}      &       +0.3156   &  0.923       &     70 \%        &  89 \%   \\
      \texttt{ival\_gauss}    &  \textbf{+0.0073}   &  0.257       &     49 \%        &  46 \%   \\

      \texttt{ival\_safe}     &       +0.0915   &  \textbf{0.173} & \textbf{32 \%}   &  92 \%   \\
      \midrule

      Risk-constrained        & MC : 0.0321 & \\
      \texttt{dt\_booles}     &       +0.4849  &  0.957       &   1162 \%      &  \textbf{100} \%  \\
      \texttt{dt\_gauss}      &       +0.0510  &  0.067       &    224 \%      &  96 \%  \\
      \texttt{ival\_gauss}    &       -0.0138  &  0.025       &     \textbf{39} \%      &  30 \%  \\

      \texttt{ival\_safe}     &  \textbf{+0.0057}   &  \textbf{0.019} &   \textbf{41} \%      &  86 \%  \\
      \bottomrule
    \end{tabular}
  \end{table}

\section{Conclusion}\label{sec:conclusion}

This paper addresses the challenging problem of efficiently and accurately estimating failure probabilities to enable risk-aware, continuous motion planning.
By developing a rigorous framework directly in continuous-time and leaning heavily on the classic study of first-exit times, it becomes straightforward to identify a lightweight approximation (\texttt{ival\_safe}) that dramatically outperforms existing methods.

Furthermore, our approximation restores a convenient Bellman structure required for optimal control, enabling practical application for a wide variety of nonlinear systems.

Ultimately, the framework and concepts introduced in this paper (particularly Section \ref{s:time_additive}) motivate a number of future investigations.
Of particular interest are robust methods of estimating anthropic belief via $\pi_t$ or $\Psi_t$, as ``survival-thus-far'' may represent a useful and yet-uncaptured source of information for aspects of autonomous decision-making not limited to risk estimation.
Also, fusing MC and Gaussian-based risk-approximations in a hybrid approach may provide the best of both worlds: high accuracy in the face of nonlinearity and compatibility with continuous optimization.

\section*{Acknowledgments}
This work was supported by the Education Office at The Charles Stark Draper Laboratory, Inc.\ and by ARL DCIST under Cooperative Agreement Number W911NF-17-2-0181.

\bibliographystyle{plainnat}
\bibliography{IEEEabrv,ref}

\ifdefined\ARXIVVERSION
\clearpage

\newcommand*\zhatt[1][t]{{}^{#1}\!\hat{z}}

\appendix
\subsection{Proofs of Stated Results}\label{app:proofs}
    \textbf{Lemma \ref{lemma:tc_props}}
    \begin{IEEEproof}
      The claims follow almost directly from the definitions; the proof here simply makes this explicit.

      Let $S_{t_1} \triangleq \{s \in [t_1, T] \cbar \z(s) \not \in D\}$ be the set of ``violation times'' after (and possibly including) $t_1$.
      Note that $\tc_{t_1} \geq t_2$ is equivalent to the statement $t_1 \leq t_2 \leq \tc_{t_1} \leq S_1$ (where the last inequality refers to a set lower bound), and the first claim is immediate.
      The contrapositive is then also established, namely
      \begin{align}
        \tc_{t_1} < t_2 &\iff \exists s \in [t_1, t_2) \text{ s.t. } \z(s) \not \in D  \\
                        &\iff S_{t_1} \cap [t_1, t_2) \neq \emptyset.
      \end{align}

      Moving on to the second claim, note that $S_{t_1} = S_0 \cap [t_1, T]$, where $S_0$ is defined analogously to $S_{t_1}$ above.
      Therefore,
      \begin{equation}
        S_0 \cap [t_1, t_2) = S_0 \cap [t_1, T] \cap [t_1, t_2) = S_{t_1} \cap [t_1, t_2).
      \end{equation}
      Thus,
      \begin{align}
        \tc_0 \geq t_1, \tc_{t_1} < t_2 &\iff t_1 \leq S_0, S_{t_1} \cap [t_1, t_2) \neq \emptyset  \\
                                        &\iff t_1 \leq S_0, S_0 \cap [t_1, t_2) \neq \emptyset  \\
                                        &\iff \tc_0 \in [t_1, t_2)
      \end{align}
      and we are done.

      The final set of claims follow directly from $z(t)$ \cadlag and $D$ closed.
    \end{IEEEproof}
    \vspace{0.5cm}

    The following auxiliary lemma does most of the heavy lifting.
    \begin{lemma}\label{lemma:xs_smooth}
      Consider process $\x(t)$ under dynamics (\ref{eq:dynamics}) and $C^2$ function $g : \manX \mapsto \R$.
      Assume there exists a finite $\bar{c}$ such that
      \begin{enumerate}
        \item $\vc{f}_t$ is pathwise-bounded \as\ and $\E \int_0^T \norm{\vc{f}_t}^2 \diff t \leq \bar{c}^2$.
        \item $\E \norm{ \mx{G}_t }^2 \leq \bar{c}^2$ for all $t$ and $\E \int_0^T \tr ( \mx{G}_t \mx{G}_t^\T ) \diff t < \infty$.
        \item Assumptions \ref{ass:loc_ind_u} and \ref{ass:paired_Lipschitz} hold with respect to $\vc{f}, \mx{G}$ and $g$.
      \end{enumerate}
      For any $t$ let $\vc{a}_t$ and $\mx{H}_t$ represent the gradient and Hessian of $g$ at $\x(t)$, and
      \begin{align*}
        \zhatt(s)    &\triangleq g\big( \x(t) \big) + \int_t^s h_t \diff \sprime + \int_t^s \vc{\sigma}_t \diff \w(\sprime)  \\
        h_s           &\triangleq \vc{a}_s^\T \vc{f}_s + \frac{1}{2} \tr \big( \mx{G}_s^\T \mx{H}_s \mx{G}_s \big)  \\
        \vc{\sigma}_s &\triangleq \vc{a}_s^\T \mx{G}_s
      \end{align*}
      be a constant-coefficient approximation to $z(s) \triangleq g\big(\x(s)\big)$ for $s \geq t$.
      Then there exists a random $M_\omega < \infty$ \as{} such that
      \begin{equation}\label{eq:xs_smooth_claim}
        P \Big( \sup_{\sprime \in [t, s]} | \zhatt(\sprime) - z(\sprime) | > (s-t)^p M_\omega \Big) \leq (s-t)^{2 - 2p}
      \end{equation}
      for any $p < 1$ and $s \geq t$.
    \end{lemma}
    \begin{IEEEproof}
      In the following, we will at several points make use of the \textbf{sum-of-squares (SOS) inequality}:
      For any $\vc{a}, \vc{b} \in \R^m$
      \begin{equation}\label{eq:sos_inequality}
        \norm{ \vc{a} + \vc{b} }^2 \leq 2 \norm{ \vc{a} }^2 + \norm{ \vc{b} }^2.
      \end{equation}

      Additionally we rely on the \textbf{It\^o isometry} for Brownian motion $\w(t)$, which follows directly from \citet[Prop.~5.2]{yong1999} (also \citet[Prop.~2.10]{karatzsas1988}) and the independent increments property.
      \begin{equation}\label{eq:ito_isometry}
        \E \norm{ \int_{t_1}^{t_2} \mx{A}_t \diff \w(t) }^2 = \E \int_{t_1}^{t_2} \tr (\mx{A}_t \mx{A}_t)^\T \diff t
      \end{equation}
      for any progressively-measurable $\mx{A}_t$.

      As a first step, we will show
      \begin{equation}\label{eq:claim1}
        \E \norm{ \x(s) - \x(t) }^2 \leq (s-t) \bar{c}^2 \quad \forall s \leq t \in [0,T].
      \end{equation}

      Starting from the left-hand-side,
      \begin{align}
      & \E \norm{\x(s) - \x(t)}^2 = \E \norm{ \int_{t}^{s} \vc{f}_\sprime \diff \sprime + \int_{t}^{s} \mx{G}_\sprime \diff \vc{w}(\sprime) }^2  \\
      & \quad \leq 2\E \norm{\int_{t}^{s} \vc{f}_\sprime \diff \sprime}^2 + 2 \E \norm{ \int_{t}^{s} \mx{G}_\sprime \diff \vc{w}(\sprime) }^2  \label{eq:pt:xterms_vanish}  \\
      & \quad = 2\E \int_t^s \int_t^s \vc{f}_{s_1}^\T \vc{f}_{s_2} \diff s_1 \diff s_2 + 2 \E \int_t^s \tr (\mx{G}_\sprime \mx{G}_\sprime^\T) \diff \sprime  \label{eq:pt:eval2}
      \end{align}
      where in (\ref{eq:pt:xterms_vanish}) we use the SOS inequality, and in (\ref{eq:pt:eval2}) we appeal to the It\^o isometry (\ref{eq:ito_isometry}).
      Evaluating the first term, application of Cauchy-Schwarz to the integrand produces
      \begin{align}
        &2\E \int_t^s \int_t^s \vc{f}_{s_1}^\T \vc{f}_{s_2} \diff s_1 \diff s_2   \nonumber \\
        &\quad \leq 2\E \int_t^s \int_t^s \max \big(\norm{\vc{f}_{s_1}}^2, \norm{\vc{f}_{s_2}}^2 \big) \diff s_1 \diff s_2  \\
        &\quad \leq 2\E \int_t^s \int_t^s \big( \norm{\vc{f}_{s_1}}^2 + \norm{\vc{f}_{s_2}}^2 \big) \diff s_1 \diff s_2  \\
        &\quad \leq 4 (s-t) \bar{c}^2
      \end{align}
      by assumption of square-integrability.
      After leveraging an assumption to bound the integrand in its second term, (\ref{eq:pt:eval2}) can now be bounded
      \begin{equation}
        \E \norm{\x(s) - \x(t)}^2 \leq 4 (s-t) \bar{c}^2 + 2 (s-t) \bar{c}^2
      \end{equation}
      and after absorbing constants into $\bar{c}$ we are done.  \checkmark

      Next we show that the drift term $h_t$ is pathwise-bounded.
      Recall that $g$ being $C^2$ ensures that $\vc{a}_t$ and $\mx{H}_t$ are continuous and therefore pathwise-bounded over $[0,T]$ \as.
      Furthermore, note that Assumption \ref{ass:paired_Lipschitz} implies that either $\mx{G}_t$ is bounded or $\mx{H}_t = 0$, so in either case $\tr (\mx{G}_t^\T \mx{H}_t \mx{G}_t )$ must be pathwise-bounded.
      From our assumption that $\vc{f}_t$ is pathwise-bounded, it clearly follows that $h_t$ is also pathwise-bounded, say by $M_\omega < \infty$ almost surely.

      Now we proceed to show the main claim (\ref{eq:xs_smooth_claim}).
      It\^o's rule \cite[Thm.~5.5]{yong1999} justifies the given definitions of $h_s$ and $\vc{\sigma}_s$, in the sense that
      \begin{equation}
        z(s) = z(t) + \int_t^s h_\sprime \diff \sprime + \int_t^s \vc{\sigma}_\sprime \diff \vc{w}(\sprime).
      \end{equation}
      The approximation gap $|z(s) - \zhatt(s)|$
      \begin{align}
        &\quad = \Bigg\lvert \int_t^s (h_\sprime - h_t) \diff \sprime + \int_t^s (\vc{\sigma}_\sprime - \vc{\sigma}_t) \diff \w(\sprime) \Bigg\rvert  \\
        &\quad \leq \int_t^s \underbrace{ \lvert h_\sprime - h_t \rvert}_{\leq 2 M_\omega} \diff \sprime + \Bigg\lvert \underbrace{ \int_t^s (\vc{\sigma}_\sprime - \vc{\sigma}_t) \diff \w(\sprime) }_{\triangleq \xi_{t,s}} \Bigg\rvert \label{eq:pt:c_omega} \\
        &\quad \leq 2 (s-t) M_\omega + \lvert \xi_{t,s} \rvert  \label{eq:split1}
      \end{align}
      where in (\ref{eq:pt:c_omega}) we've applied our pathwise-bound on $h_s$.
      We've also introduced the shortcut $\xi_{t,s}$ to represent the It\^o integral.

      Note that under Assumptions \ref{ass:loc_ind_u} and \ref{ass:paired_Lipschitz},
      \begin{align}
        \norm{\vc{\sigma}_s - \vc{\sigma}_t} &= \norm{\mx{G}_s^\T \vc{a}_s - \mx{G}_t^\T \vc{a}_t }  \\
                                              &= \norm{ \mx{G}_s^\T (\vc{a}_s - \vc{a}_t) + (\mx{G}_s - \mx{G}_t)^\T \vc{a}_t}  \\
                                              &\leq 2 L \norm{x(s) - x(t)}^2.  \label{eq:sigma_lipschitz}
      \end{align}
      It follows that the second moment of $\xi_{t,s}$ can be bounded via the It\^o isometry (\ref{eq:ito_isometry}) as
      \begin{align}
        \E |\xi_{t,s}|^2 &= \E \int_t^s \norm{ \vc{\sigma}_s - \vc{\sigma}_t }^2 \diff \sprime  \\
                        &\leq \E \int_t^s 4 L^2 \norm{x(\sprime) - x(t)}^2 \diff \sprime  \label{eq:use_sigmatilde_bound} \\
                        &\leq 4 L^2 \int_t^s (\sprime - t) \bar{c}^2 \diff \sprime  \label{eq:xi_mom2_bound} \\
                        &\leq 2 (s - t)^2 \bar{c}^2 L^2
      \end{align}
      where in (\ref{eq:xi_mom2_bound}) we've appealed to the recently-proven (\ref{eq:claim1}).

      We then appeal to Doob's maximal inequality \cite[Thm.~4.5]{yong1999} for the martingale $\xi_{t,s}$, which gives for any $p$
      \begin{align}
        P\Big( \sup_{\sprime \in [t, s]} |\xi_{t,s}| > (s-t)^p \bar{c} L \Big) &\leq \frac{ \E |\xi_{t,s}|^2 }{ \big((s-t)^p \bar{c} L \big)^2}  \\
                                                                                &\leq 2 (s-t)^{2 - 2p}.  \label{eq:p_bound}
      \end{align}
      Let $E \subset \Omega$ refer to the event that the condition in the left-hand-side of (\ref{eq:p_bound}) fails.
      Returning to (\ref{eq:split1}), it is clear that for any $\omega \in E$ and $p < 1$,
      \begin{align}
        \sup_{\sprime \in [t,s]} |z(\sprime) - \zhatt(\sprime)| &\leq 2 (s-t) M_\omega + (s-t)^p \bar{c} L  \\
                                                                &\leq (s-t)^p M_\omega
      \end{align}
      where we've assumed WLOG that $M_\omega \geq \bar{c} L$.
      Recalling from (\ref{eq:p_bound}) that $P(E^c) \leq 2 (s-t)^{2 - 2p}$, we are done.
    \end{IEEEproof}
    \vspace{0.5cm}

    \textbf{Theorem \ref{thm:puc_implies_P_conv}}
    \begin{IEEEproof}
      Define the approximating and target events, respectively, as
        \begin{align*}
          &A_n \triangleq \big\{\omega \in \Omega \,|\, \max_j \sup_s \zjn(s) \leq 0 \big\} \text{ and } \\
          &E \triangleq \big\{\omega \in \Omega \,|\, \max_j \sup_s z_j(s) \leq 0 \big\}.
        \end{align*}
        Then $\eps_n$-PUC of $\zn(s)$ and $\vc{z}(s)$ implies that there exist subsets $S_n \subseteq \Omega$ such that for any $\omega \in S_n$
        \begin{equation}\label{eq:puc_condition}
          \norm{ \zn(s) - \vc{z}(s) } \leq \eps_n M_\omega
        \end{equation}
        and that $P(S_n) \rightarrow 1$.
        Define the joint events
        \begin{equation*}
          B_n \triangleq A_n \cap S_n \text{ and } E_n \triangleq E \cap S_n.
        \end{equation*}
        By the additivity of measure over disjoint sets, it is easy to verify that
        \begin{equation*}
          P(A_n) = P(B_n) + P(A_n \cap S_n^c)
        \end{equation*}
        and because $P(A_n \cap S_n^c) \leq P(S_n^c) \rightarrow 0$, it must follow that $\lvert P(A_n) - P(B_n)\rvert \rightarrow 0$.
        An analogous argument establishes that $P(E_n) \rightarrow P(E)$.
        Thus, if we can simply show that $\lvert P(B_n) - P(E_n) \rvert \rightarrow 0$, we will be done.

        Again, using the additivity of measure, we have
        \begin{align*}
          &P(B_n) = P(B_n \cap E_n) + P(\underbrace{B_n \cap E_n^c}_{\triangleq C_n}) \text{ and }  \\
          &P(E_n) = P(E_n \cap B_n) + P(\underbrace{E_n \cap B_n^c}_{\triangleq D_n})
        \end{align*}
        and it will be sufficient to show that the probability of ``disagreement'' $P(C_n) + P(D_n) \rightarrow 0$.

        From here, we can consider each of the $m$ elements of the vector processes $\z(s)$ and $\zn(s)$ independently.
        This is because $P(C_n)$ can be written
        \begin{align}
          &P\big(S_n \bigwedge \max_j \sup_s \zjn(s) \leq 0 \bigwedge \max_j \sup_s z_j(s) > 0 \big)  \\
          &\leq \sum_{l = 1}^m  P(S_n \bigwedge \max_j \sup_s \zjn(s) \leq 0 \bigwedge \sup_s z_l(s) > 0 \big)  \label{eq:pt:union_bound}  \\
          &\leq \sum_{l = 1}^m  \underbrace{ P(S_n \bigwedge \sup_s \zjn[l](s) \leq 0 \bigwedge \sup_s z_l(s) > 0 \big) }_{P(C_{n,l})}  \label{eq:pt:relax2}
        \end{align}
        where in (\ref{eq:pt:union_bound}) we've used the union bound and in (\ref{eq:pt:relax2}) we've relaxed the probabilistic statement.
        So if we can show that if each term $P(C_{n,j})$ decays to zero (which involves only the $j$-th channel) decays to 0, then $P(C_n)$ must also decay to 0.
        Furthermore, the PUC criterion for the vector process implies a similar criterion for each of the $j$ elemental processes.
        An analogous decomposition can be performed for $D_n$, so for the remainder we consider each of the $j$ channels independently.

        Define the random variable $e_n \triangleq \sup_s \norm{ \zn(s) - \z(s) }$, allowing $e_{n} = \infty$ if the supremum does not exist.
        Then $z_j(s) \in [\zjn(s) - e_n, \zjn(s) + e_n]$ for all $j,s$.

        First, consider $C_{n,j}$, which by construction is a subset of $S_n$, implying (\ref{eq:puc_condition}) and therefore $e_n(\omega) \leq \eps_n M_\omega$ holds in $C_{n,j}$.
        \begin{align*}
          P(C_{n,j})  &= P\big(S_n \bigwedge \sup_s \zjn(s) \leq 0 \bigwedge \sup_s z_j(s) > 0 \big)  \\
                      &\leq P\big( S_n \bigwedge \sup_s z_j(s) \in (0, e_n] \big)  \\
                      &\leq P\big( \sup_s z_j(s) \in (0, \eps_n M_\omega] \big)
        \end{align*}
        Because the interval approaches the empty set monotonically as $n \rightarrow \infty$, it is clear that $P(C_{n,j}) \rightarrow 0$ for each $j$.

        Similarly for $D_{n,j}$, we have
        \begin{align*}
          P(D_{n,j}) &= P\big(S_n \bigwedge \sup_s z_j(s) \leq 0 \bigwedge \sup \zjn(s) > 0 \big)  \\
                     &\leq P\big( S_n \bigwedge \sup_s z_j(s) \in (-e_n, 0] \big)  \\
                     &\leq P\big( \sup_s z_j(s) \in (-\eps_n M_\omega, 0] \big)
        \end{align*}
        In this case, the latter probability converges monotonically to $P\big(\sup z_j(s) = 0\big)$, which by assumption equals 0.
        Therefore, $P(D_n)$ also converges to 0.
        Altogether, we have that
        \begin{align*}
          \lvert P(A_n)  - P(E) \rvert &\leq \lvert P(A_n) - P(B_n) \rvert  \\
                                       &~~+ \lvert P(B_n) - P(E_n) \rvert  \\
                                       &~~+ \lvert P(E_n) - P(E) \rvert \rightarrow 0
        \end{align*}
        completing the proof.
    \end{IEEEproof}
    \vspace{0.5cm}

    \textbf{Proposition \ref{prop:yhatn_puc}}
    \begin{IEEEproof}
      From Lemma \ref{lemma:xs_smooth} we already have a similar statement for each element $\yjhatn(s)$ and $y_j(s) = g_j\big(\x(s)\big)$.
      A fundamental relation between the $l_1$ and $l_2$ norms in $\R^m$ gives
      \begin{equation*}
        |\yjhatn(s) - y_j(s)| \leq \alpha \,\, \forall j \, \implies \, \norm{\yhatn(s) - \y(s)} \leq m \alpha
      \end{equation*}
      for any $\alpha > 0$.
      For any timestep $t_i \in \pT$, Lemma \ref{lemma:xs_smooth} states that for any $j$
      \begin{equation}\label{eq:segment_p_bound}
        P\Big( \underbrace{ \sup_{s \in [t_i, t_{i+1})} |\yjhatn(s) - y_j(s)| > \delta_n^p M_\omega}_{\triangleq S_{n,i,j}^c} \Big) \leq \delta_n^{2 - 2p}.
      \end{equation}
      Letting $S_{n,i,j}$ be the event that error is ``small'' for the $j$-th channel at the $i$-th timestep, define $S_n$ as the event that all such errors are small
      \begin{equation*}
        S_n \triangleq \bigcap_{i = 0}^{k_n-1} \bigcap_{j = 1}^{m} S_{n,i,j}.
      \end{equation*}
      Then by the union bound the probability of $S_n$ failing is upper-bounded
      \begin{equation}
        P(S_n^c) \leq \sum_{i = 0}^{k_n-1} \sum_{j = 1}^m P(S_{n,i,j}^c) \leq m k_n \delta_n^{2 - 2p} \leq m c T \delta_n^{1 - 2p}.
      \end{equation}
      Taking $p = \frac{1}{4}$, $P(S_n^c)$ clearly goes to 0 as $\delta_n \rightarrow 0$.
      Furthermore, by construction $\omega \in S_n$ implies that
      \begin{equation}
        \norm{ \yhatn(s) - \y(s) } \leq m \delta_n^\frac{1}{4} M_\omega \quad \forall s \in [0,T].
      \end{equation}
      and defining $\eps_n \triangleq m \delta_n^\frac{1}{4}$ the proof is complete.
    \end{IEEEproof}

\subsection{Reduced Quadratures: Second-Order Gaussian Systems}\label{app:red_quad}
\newcommand*\p{\vc{p}}
\newcommand*\vjp{v_{j,\vc{p}}}
\newcommand*\Rnp{{\R^{n_p}}}

  We would like to evaluate the following $n_x$-dimensioned integral for each constraint $j$
  \begin{align}
    &\Phiz[\yjhat](t_i, t_{i+1}; \mu_i) =   \label{eq:quad_full}  \\
    &\quad \int_{\manX} \Ind(\x \in \Xsafe) \psi\big(g_j(\x), h_j(\x), \norm{\vc{\sigma}_j(\x)}, \Dt_i\big) \diff \mu_i(\x)  \nonumber
  \end{align}
  which requires some form of numeric quadrature evaluation.

  Note that discrete-time direct methods based on (\ref{eq:dt_booles}) also appeal to evaluation of a similar integral
  \begin{equation}\label{eq:P_violation}
    P\big( \x(t_i) \in \cXsafe \big) = \int_\manX \Ind(\x \not \in \Xsafe) \diff \mu_i(\x).
  \end{equation}
  Upon application, the dimensionality of both (\ref{eq:quad_full}) and (\ref{eq:P_violation}) can often be significantly reduced, particularly when $\mu_i$ refers to a multivariate Gaussian density.

  Many physical systems evolve as nonlinear ``integrators'', that is, the state space can be partitioned as $\x = (\vc{p}, \vc{p}, \star)$ where $\star$ represents ``other'' state, such that
  \begin{equation}\label{eq:p_dynamics}
    \diff \vc{p}(t) = \vc{v}(t) \diff t + \mx{G}_p\big(\vc{p}(t)\big) \diff \vc{w}(t),
  \end{equation}
  where $\vc{p}$ can be interpreted as a ``position'' vector in a workspace $\R^{n_p}$, and $\mx{G}_p$ refers to the corresponding rows of $\mx{G}$.
  It is also often the case that safety constraints depend exclusively on this position state; that is $\vc{g}(\x) = \vc{g}(\vc{p})$.
  Note that (\ref{eq:p_dynamics}) is furthermore independent of control.

  In light of (\ref{eq:p_dynamics}) and the position-only dependence of $g_j$, we can write (\ref{eq:P_violation}) and (\ref{eq:quad_full}) as
  \begin{equation}\label{eq:P_violation_p}
    P\big( \x(t_i) \in \cXsafe \big) = \int_{\R^{n_p}} \Ind\big( \vc{g}(\p) \not \in \setOm \big) \diff \mu_i(\vc{p})
  \end{equation}
  and
  \begin{align}
    &\Phiz[\yjhat](t_i, t_{i+1}; \mu_i) =  \label{eq:quad_pv} \\
    &\quad \int_\Rnp\!\!\!\! \Ind\big( \vc{g}(\p) \in \setOm \big) \int_\Rnp \!\!\!\! \psi\big( \p, h_j(\p, \vc{v}) \big) \diff \mu_i(\vc{v} \cbar \p) \diff \mu_i(\p)  \nonumber
  \end{align}
  where with a slight abuse of notation $\mu_i\big(\dot{\vc{p}} \cbar \vc{p} \big)$ and $\mu_i(\vc{p})$ represent the corresponding conditioned or marginal distributions.
  Recalling from (\ref{eq:h_def}), $h_j(\p, \vc{v})$ is affine in $\vc{v}$
  \begin{equation}
    h_j(\p, \vc{v}) = \frac{1}{2} \tr\big( \mx{G}_p(\p)^\T \mx{H}(\p)_{j,p} \mx{G}_p(\p) \big) + \vc{a}_j^\T(\p) \vc{v}.
  \end{equation}
  When $\mu_i$ represents a multivariate Gaussian, we can consider the scalar drift $\vjp \triangleq \vc{a}_j^\T(\p) \vc{v} \in \R$ and then (\ref{eq:quad_pv}) can be written
  \begin{align}
    \int_\Rnp\!\!\!\! \Ind\big( \vc{g}(\p) \in \setOm \big) \int_\R \!\! \psi(\p, \vjp) \diff \mu_i(\vjp \cbar \p) \diff \mu_i(\p)  \label{eq:quad_ps}
  \end{align}
  which can now be evaluated in dimension $n_p+1$.

  \subsubsection{The Locally-Deterministic Case}\label{app_ss:local_determinism}
    Consider the case $\mx{G}_p = 0$; this corresponds to the case that no noise is ``injected'' directly into the $\p$ channel of the continuous dynamics.
    This is in fact commonly the case for physical systems, as it is often assumed that disturbances act in the space of forces or torques (\ie, the $\vc{v}$ dynamics).
    In this case, $\vc{\sigma_j} = \vc{a}_j^\T \mx{G}_p = 0$ uniformly, and it is straightforward to verify $h_j(\p,\vc{v}) = \vjp$ and that
    \begin{equation}
      \psi(z, v, 0, \Delta t) = \Ind\big( z + v \Delta t > 0 \big).
    \end{equation}
    In this case (\ref{eq:quad_ps}) simplifies to
    \begin{equation}\label{eq:quad_p}
      \int_\Rnp\!\!\!\! \Ind\big( \vc{g}(\p) \in \setOm \big) \phi\Big(\frac{\vc{g}(\p)}{\Delta t}; -\bar{v}_j(\p), \sigma_j(\p) \Big) \diff \mu_i(\p)
    \end{equation}
    where $\bar{v}_j(\p)$, $\sigma_j(\p)$, and $\phi$ refer to the mean, standard deviation, and CDF of the normal distribution $\mu_i(\vjp \cbar \p)$.
    This implies that $\Phiz[\yjhat](t_i, t_{i+1}; \mu_i)$ can be evaluated over only $n_p$ dimensions.

\subsection{Simulated Dubin's System}\label{app:sys}
  Our experimental validation relies on a simulated second-order Dubin's car with state $\x = (\p, \vc{v}, \theta, \omega) \in \R^6$, and control $\u = (c, \alpha)$.
  $\theta$ and $\omega$ capture the heading angle and angular rate, respectively, while $c$ refers to forward-pointing thrust and $\alpha$ the angular acceleration.
  The time-invariant, nonlinear dynamics are given by
  \begin{equation}
    \vc{f}(\x, \u) = \begin{bmatrix} \dot{\p} \\ \dot{\vc{v}} \\ \dot{\theta} \\ \dot{\omega} \end{bmatrix}
                   = \begin{bmatrix} \vc{v} \\ c \mx{R}(\theta) \vc{e}_1 \\ \omega \\ \alpha \end{bmatrix}
  \end{equation}
  with constant noise matrix
  \begin{equation}
    \mx{G} = 0.05 \begin{bmatrix} \mx{0} & 0 & 0 \\
                                  \mx{I}_2 & 0  & 0 \\
                                  \mx{0}  & 0.1 & 0 \\
                                  \mx{0}  &  0  & 1 \end{bmatrix}.
  \end{equation}

  Note that, for this system, we can compute $P(\x(t_i) \not \in \Xsafe)$ (required for discrete-time methods) with a quadrature of dimension $n_p - 1 = 1$.
  Fortunately, because this system satisfies the ``locally-deterministic'' criteria described in Appendix~\ref{app_ss:local_determinism} evaluating $\Phiz[\yhat]$ requires $n_p = 2$ dimensions, an increase of only one.

  To simulate output-feedback control, we introduce a simple observation process at each timestep $k$
  \begin{equation}
    \vc{y}_k  = \x(t_k) + \vc{\nu}
  \end{equation}
  where $\vc{\nu} \sim \mc{N}(0, 0.0001 \mx{I}_6)$.
  Then a corresponding LQG-style controller is constructed to track a nominal trajectory that also defines the linearization point.
\fi

\balance

\end{document}