Learning Density Distribution of Reachable States for Autonomous Systems
Authors: Yue Meng
Dawei Sun
Zeng Qiu
Md Tawhid Bin Waez
Chuchu Fan
Venue: CoRL2021 Poster (2021)
Keywords: Reachability Density Distribution
Learning Density Distribution
Liouville Theorem

Abstract: State density distribution, in contrast to worst-case reachability, can be leveraged for safety-related problems to better quantify the likelihood of the risk for potentially hazardous situations. In this work, we propose a data-driven method to compute the density distribution of reachable states for nonlinear and even black-box systems. Our semi-supervised approach learns system dynamics and the state density jointly from trajectory data, guided by the fact that the state density evolution follows the Liouville partial differential equation. With the help of neural network reachability tools, our approach can estimate the set of all possible future states as well as their density. Moreover, we could perform online safety verification with probability ranges for unsafe behaviors to occur. We use an extensive set of experiments to show that our learned solution can produce a much more accurate estimate on density distribution, and can quantify risks less conservatively and flexibly comparing with worst-case analysis.

Learning Density Distribution of Reachable States for Autonomous Systems
Yue Meng
Dawei Sun UTCC United States
Zeng Qiu Ford
MIT
United States
United States cqiu1@forde com
mengyuečmit.c. edu
daweis20illinoii edu
Md Tawhid Bin Waez Ford
Chuchu Fan MIT
United States
United States chuchu@mit edu
mwaezêford.c com
Abstract: State density distribution, in contrast to worst-case reachability, leveraged for safety-related problems to better quantify the likelihood the risk for potentially hazardous situations. In this work, propose data-driven method to compute the density distribution of reachable states for nonlinear and even black-box systems. Our semi-super approach learns system dynamics and the state density jointly from trajectory data, guided by the fact that the state density evolution follows the Liouville partial differential equation. With the help of neural network reachability tools, our approach can estimate the set of all pos- sible future states well their density. Moreover, we could perform online safety verification with probability ranges for unsafe behaviors to occur. We use extensive set of experiments show that our learned solution can produce a much more accurate estimate on density distribution, and can quantify risks less conservatively and flexibly comparing with worst-case analysis.
Keywords: Reachability Density, Machine Learning, Liouville Theorem
Introduction
Reachability analysis has been a central topic and the key for the verification of safety-critical autonomous systems. The majority of existing reachability approaches either compute worst- case reachable sets [1, 3, 4, 5, 6), or use Monte Carlo simulation to estimate the reach- able states with probabilistic guarantees [7, 8, 9, 2]. There several obvious disadvantages of such methods. 1). Worst-case methods, especially those for nonlinear systems, often over- approximate the reachable sets and produce very conservative results (due to convex set representa- tions [10, 11, 12, 13), wrap effects [14, 15], low-order approximation [3, 16), etc.), not to say they usually very computationally expensive (also known the curse of dimensionality). 2). Exist- ing methods do not care about reachable state concentration and produce the same reachable state estimations for different distributions of initial states if those distributions share the same support 1 That is, current approaches do not compute which states are more likely to be reached. For example, in Fig. 1, the state density distribution of Van der Pol oscillator evolves over time and concentrates on certain states (the highlighted part in the figure), even if the initial states are uniformly distributed. one uses an existing worst-case reachability algorithm, most likely the results in Fig. (b)(c) will show almost the entire black region inside the highlighted part is reachable, as those methods often use convex sets represent the reachable sets.
oo
(a) Os. (b) 0.35s. (c) 0.80s.
Figure 1: Reachability density of Van der Pol.
1 'The probabilistic guarantees the sampling based methods rely form of the initial distri- bution. However, the final reachable sets estimate the same for different distributions with the same support.
5th Conference on Robot Learning (CoRI. 2021), London, UK.

---

In this paper, we aim to tackle the above problems and propose a learning- based method that can compute the density distribution of the reachable states from given initial distributions. The state density is much more powerful than worst-case reachability and can better quantify risks. Our proposed method is based on the Liouville theorem [17, 18, 19], which is from classical Hamiltonian mechanics and asserts that the state density distribution function (which is the measurement of state concentration) is constant along the trajectories of the system. Given an autonomous system á = f(x) that is locally Lipschitz continuous, the evolution of the state density p(r, t) (i.e. the density of state x at time t) is characterized by a Liouville partial differential equation (PDE). We learn the state density as neural network (NN) while respecting the laws of physics described by the nonlinear Liouville PDE, in a way that is similar to the physics-informee NN [20].
Furthermore, we make two major improvements to make the learned density NN suitable for ver- ification. 1). Instead of learning / for fixed initial distribution p( 0), we learn the density concentration function which specifies the multiplicative change of density from any p( 0). 2). We use single NN to jointly learn both the reachable state ((r,,) t) and its corresponding density. Moreover, we use Reachable Polyhedral Marching (RPM) [21]-an exact ReLU NN reachability tool-to parse our learned NN as linear mappings from input polyhedrons to output polyhedrons. Using such parsed polyhedrons, we can perform online forward and backward reachability analysis and get the range of density bounds [pmiu: Pmax] for each output polyhedron. Together, our method can perform online safety verification by computing the probability of safety (instead of single yes or no answer) under various initial conditions, even with unknown system dynamics f (i.e. a black- box system where only have access to a simulator of ztt)). In this way, it also has the potential to collect environmental data in run-time and update its distribution for online safety verification.
We conduct experiments on 10 different benchmarks covering systems from low dimension aca- demic examples to high dimension black-box simulators equipped with either hand-crafted or NN controllers. Surprisingly, without using any ground-truth density data in the learning process, our approach can achieve up to 99.89% of error reduction in KL divergence with respect to the ground- truth value, when compared to sampling- Processes. Morcover, our learned density concentration function can also be used for reachability distribution analysis. We show that in several systems, more than 90% of the states can actually just reside in a small region (less than 10% of the volume of the convex hull for those states) in the state space, which also points out the conservativeness of the worst-case reachability analysis methods in terms of quantifying risks. We also show that different initial distributions can lead to a drastic change in the safety probability, which can help in cases when unsafe is inevitable but can be designed to happen with a very low probability.
methods like kernel density estimation and Gaussian
Our major contributions are: (1) we are the first to provide explicit probability density function for reachable states of dynamical systems characterized as ordinary differential equations; this den- sity function can be used for online safety analysis, (2) we propose the first data-driven method to learn the state density evolution and give accurate state density estimation for arbitrary (bounded support) initial density conditions without retraining fine-tuning, and (3) we use variety of ex- amples to show the necessity to perform reachability distribution analysis instead of pure worst-case reachability, to flexibly and less conservatively quantify safety risks.
Related work Reachability analysis, especially worst-case reachability using sampling- based methods or set propagation, has been studied for decades. The literature on reachability has been extensively studied in many surveys [1, 6, 22, 23]. Here we only discuss few closely related works.
Hamilton Jacobian PDE has been used to derive the exact reachable sets in [24, 1, 5]. However, the HJ-PDE does not provide the density information. Many data-driven approaches can com- pute probabilistic reachable sets using scenario optimization [7, 25]. convex shapes [8, 9, 26), support vector machines [27, 28]. kernel embedding [29], active learning [30], and Gaussian process [2]. However, the probabilistic guarantees they provide are usually in the form that Prob(z(t) E estimated set) > 1 --6 with enough samples, instead of the state density dis- tribution. [31] estimates human state distribution using a probabilistic model but requires state discretization. [32] uses the Liouville equation to maximize the backward reachable set for only polynomial system dynamics. In [33] the authors compute the stochastic reachability by discretiz- ing stochastic hybrid systems to Markov Chains (MC), then perform probabilistic analysis on the discretized MC. The closed-form expression of the probability requires integrals over the whole state space hence is computation-heayy and cannot be used for online safety check. The closest to ours
2

---

is [34] where Perron-Frood Frobenius and Koopman operators are learned from samples of trajectories. Then, the learned operators can be used to transform the moments of the distribution over time. The distribution at time I. is then recovered from the transformed moments. However, as they use mo- ments up to a specific order to represent a distribution, even the learned operators are perfect, the estimation error of the distribution might not be zero. Also, such a moment-based method is hard to scale to aage-dimennsionl systems, Recently, there is also a growing interest in studying the (worst- case) reachability of NN [35, 36, 37, 38, 21] or systems with NN controllers [39, 4, 40, 41, 42]. In this paper, we use [21] to parse our learned NN as a set of linear mappings between polyhedrons for online reachability computation, but this step can be replaced with other NN reachability tools.
To measure the probability distribution in the reachable sets, the most naive approach is to use his- tograms or kernel density estimation [19], similar to the Monte-Cate method used in dispersion analysis [43]. But this could lead to poor accuracy and computational scalability [44]. Another ap- proach propagates the uncertainty by approximating the solution of the probability density function (PDF) transport equation [45], which could still be time-consuming due to the optimization process performed at each time step. Our approach finds the PDF transport equation by solving the Liou- ville PDE 2 using NN. Similar ideas have been explored in [18]. However, the density NN in [18] was learned solely for fixed initial states distribution and therefore, cannot be used for online pre- diction. Instead, we jointly learn the reachable states and density changes, and can perform online reachability density computation for any initial states distribution.
The idea of using deep learning solve PDEs can be traced back to the 1990s [46, 47, 48], where the solutions of the PDE on a priori fixed mesh are approximated by NN. Recently, there a growing interest in the related sub-fields including: mesh-free methods in strong form [20, 49, 50] and weak form [51, 52], solving high-dimension PDEs via BSDE method [53], solving parametric PDE [54, 55] and stochastic differential equations [56, 57], learning differential operators[58,, 59, 60], and developing more advanced toolboxes [61, 62, 63]. Our idea for solving Liouville PDE along trajectories is similar to [53] but without the stochastic term.
2 Preliminaries
We denote by IR, R220 R", R=*= the sets of real numbers, non-negative real numbers, n--dmmnnsional real vectors, and n x n real matrices. We consider autonomous dynamical systems of the form :(t) = fx((t), where for all t R, xxt) E X S E is the state and X' is a compact set that represents the state space. We assume that f Rdd + Ied is locally Lipschitz continuous. The solution of the above differential equation exists and is unique for a given initial condition To- We define the flow map & X x ER Y X such that (<ro, t) (also written as x(t) for brevity) is the state at time starting from Zo at time O. Note that system parameters can be easily incorporated as additional state variables with time derivative to be O.
We analyze the evolution of the dynamical system by equipping it with a density function p X x R + R22, which measures how states distribute in the state space at a specific time instant. A larger density p(z, t) means the state more likely to reside around x at time t, and vice versa. The density function is completely determined by the underlying dynamics (i.e., function f) and the initial density map Po X r> RR200 Specifically, given a Po. the density function p solves the following boundary value problem of the Liouville PDE [17]:
0p Ot
+ (p. 00 p(x, 0) = Po(x),
where V f) Dd Dizl alp-Al: is the divergence for the vector field p f. and [. takes the 0[z).
i-th coordinate of a vector.3 Intuitively, as shown in [17], Liouville PDE is analogous to the mass conservation in fluid mechanics, where the change of density गे at one point is balanced by the total flux traversing the surface of a small volume surrounding that point. It is hard to solve the Liouville PDE for closed from of the density function p. However, it is relatively easy to evaluate
@e
IIn the absence of process noise, this PDF transport equation reduces to stochastic Liouville equation [17] 3 "A general form of Liouville PDE can have a non-zero term on the right- hand side indicating how many (new) states appear or exit from the system during the run time [19]. In all the systems we discuss here, there is no state entering (other than the initial states) or leaving the system, so the right-hand side of Eq. (1) is zero. The total density of the systems consider is invariant over time.
3

---

the density along trajectory of the system. To do that, we first convert the PDE into an ODE as follows. Considering a trajectory & ((z,,t), the density along this trajectory an univariate functionn of I, i.e., p(t) = p (zo, 1), l). If we consider the augmented system with states (x, pl. from Eq. (1) we can easily get the dynamics of the augmented system [19]:
i
f(x) -V (VV-((x))
(2)
P
Therefore, to compute the density an arbitrary point 2t E X' at time T, one can simply proceed as follows: First, find the initial state on = ((TT, -7) using the inverse dynamics -f- Then, solve the Eq. (2) with initial condition [ (x The solution time T just gives the desired density value. However, such a procedure only gives the density at a single point, therefore, cannot be used in reachability analysis. Instead, we need to compute the solution of Eq. (2) for set of initial conditions and use that to compute the reachable sets. To achieve this, we use an NN with ReLU (Rectified Linear Unit) activation functions [64] to jointly approximate the flow map and the density concentration function, as will be shown in the next section.
3 Density learning and online reachability density computation
Let us take closer look Eq. (2). For an initial condition [mo, Po(xr)), the closed form of the solution P is p t) (o(20) exp (-Jo'
J v ( Interestingly, the solution of is linear in the initial condition Po- The gained part denoted by G( t) == exp V- of Po- This mapping is completely determined by the underlying dynamics and we call the density concentration function. With G' in hand, the density at an arbitrary time and state can be quickly computed from any Po. However, G is obviously hard to compute. Therefore, we use NN to approximate the density concentration function. In addition to G, we also use NN to learn the flow map 2 which will be a necessity for computing the reachable set distribution shown in Sec. 3.2.
r))dT) is a function of the initial state zo and time t. but is independent
3.1 System dynamics and density learning framework
Let the parameterized versions of the flow map and the density concentration function be P,, and Gg respectively, where w and 0 are parameters. To train the neural network, we construct a dataset by randomly sampling N trajectories of the system: D Dtrain = N- in T time steps (with time interval At): & { 0), (zi, At), ( 1, (T 1))))} where z) 4 (tt) Then, the goal of the learning is to find parameters w and satisfying
L (ab kat) of. = 0, Vi,k, 8G (G). + Go kDt) (V f(x(( = Vi,k, dt
(3)
where the first constraint is for the flow map estimation, and the second constraint enforces the Liouville equation for all the data points.
As for the implementation, with ReI.U activations. To ensure numerical stability as G is an exponential function, we add a nonlinear transform from the NN output to the density concentration function:
model b.,, and Ge jointly a fully colly-connected neural network NN(:.) :)
Go( t) exp(t NN No (Zo. ((oott) t)) = exp(t 2(20. t)),
(4)
T (1o. t) NN NN] (zo, t),
where NN [+j+i1] ij is to choose the i, + 1, l.....jjjj........ j-th dimensions from the output of the NN, and a is the intermediate density estimation from the NN. In this way, we guarantee that the density concentra- tion function at = 0 is always 1. We optimize our NN via back propagation with the loss function:
NN
C- Σ (s 2,K
1211 + - C(Go k(t) + Go kAt) (V f(xt)))
(5) where the first term denotes the state estimator square error [65], the second term indicates how far (in the sense of L2-norm) the solution deviates from the Liouville Equation, and À balances
1.k

---

these two loss terms. We approximate the time derivative of the density concentration function by Go k&t) = [Go (k + 1)At) - Go kat)] / At. Our method can also work for black- box systems if we approximate V f() numerically. With some tools from statistical learning theory, we can show that with a large enough number N of samples, the learned flow map and the density concentration, function can be arbitrarily accurate. A formal proof is shown in Appendix A.
3.2 System reachable set distribution computation via NN Reachability Analaysis
In the last section, we have learned an NN that can estimate the density at single point given the initial density. In this section, we will boost this single-point estimation to set-based estimation by analyzing the reachability of the learned NN. To compute the set of all reachable states and the corresponding density from a set of initial conditions, we use the Reachable Polyhedron Marching (RPM) method [21] to further process the learned NN for G and Φ. RPM is seyyhhhddbbbsee approach for exact reachability analysis for ReLU NNs. It partitions the input space into polyhedral cells that in each cell the ReLU activation map does not change and the NN becomes fixed affine mapping. With the cells and the corresponding affine mapping on each cell, the exact reachable set for set of input can be quickly evaluated. Backward reachability analysis can also be performed by computing the intersection of the pre-image of a query output set with the input polyhedron cells.
System forward reachable set with density. Recall that the input of the NN in Sec. 3.1 consists of the initial state zig and time t. For the simplicity of comparison with other methods, we only estimate the reachable set and the density at given fixed time instances t and fix the last element of the input of the NN. Thus for each time step I, the input polyhedral cells generated from the RPM will be set of linear inequality constraints on 2, and those input cells together with the set of the affine mappings and output polyhedral cells can be represented as: {(Ax, by. Ck. dx, Es. f where in each input cell H₃ := {u E | S ). the NN becomes an affine mapping y = C&U + du. and thus the image of the input cell is also a polyhedron Mr = fy E T S fie]
Recall that the first dimension of our NN output estimates the density concentration function 2, and the rest dimensions estimate the state F, thus the output cell can be written Mr {(z, æ:) Ez 2/² I fk). By projecting it to the state space, we get the reachable set of the system, i.e., Rg := {I E x2, E Mr.. Then, in each cell Mg. we evaluate the lower and upper bounds of 2, and denote them by zk, mén and zk.maz. The density bound for cell Mr is then computed as
Pk.miin = po(x) exp(t- =k,min),
(pk.mar = (T(x exp(1: =k,mar);
where Fy is the center of Hy. Finally the system forward reachable set is a union of projected output N
polyhedral cells: U R2 where each cell Rg is associated with a density bound [pi,min Pit maz] 1=1
System reachable set probability computation. Given an initial state distribution, want to figure out the probability distribution of the system forward reachable sets, as well as the probability for the states land into query set (c.g., the query set could be the unsafe region).
For an arbitrary initial probability density function (whose support is bounded), we can apply RPM to partition its support into cells* as in the last section. Finally, we obtain the reachable sets with bounded state densities {(Hp, Pk man maz) NN and the probability bound in each cell is:
pk
(
pmmin = Vol(H*)Pk Vol( H*)Pp"r
P.max
where Vol(-) computes the volume for a polyhedron. By computing for all input cells {H& H can derive the system forward reachable set and the corresponding probability bound as {HE, pmmin pmma)) ) } The backward reachable set probability can be computed in a similar fashion, by checking the intersection between the query output region and the output cells derived by RPM, computing each intersection's probability range by its volume and density bound and finally aggregating the probability of all intersections. Detailed computation is shown in Appendix B.
1k=1"
"We can always further divide those input cells to make the bound tighter/ more precise, while still guaran- teeing the Neural Network on each cell can be seen as an affine transformation.
'The RPM method cannot handle systems with higher than 4 dimensional state space in our experiments. But the technique discussed in See. 3.2 can work with any set-based NN reachability tools with slight modifica-
5

---

<L-divereence inormelized)
10*
Systete
Dum.
Consol
Vin der Pet Oscillake (slp)
KDD
Kraichaun Orizag system (hop) [66]
imegrancr (din))
100
Flstogram
Histogram (1U8P)
Inverted gendulum (jend) [67]
LQR
SGPD
Gound navigation (mphos)
JACTEST car taaking system (car) [68]
Trocking
10a
Qualkotur costrol syseet (qual) [42]
NN
Adaptive sruise contrl
[4)
F-16 Ground coliaioo avridunce (gcas) (60]
Hybrid
100
8-Cae plation sysrm (toon) (70)
16
NN
vGp
kop
pend
robot
car
quad
gcas
toon
Table 1: Benchmarks from low dimen- sion academic models complex systems
Figure 2: KI. divergence with ODE-simulated density. Our method consistently outperforms other baselines. The histogram method cannot estimate the density for "acc", "geas" and "toon". ".
with handcrafted or NN controllers.
4
Experimental evaluation in simulation
Here we show the benefits of learning density concentration function from Liouville PDE in state density and reachable set distribution estimation. More experiments are provided in Appendix C-H.
4.1 Implementation details
Datasets: We investigate 10 benchmark dynamical systems as reported in Table 1. These bench- mark systems range from low dimension academic models (2~3 dimensions) to complex and even black-box systems (13~ 16 dimensions) controlled by handcrafted or NN controllers. All controllers (except for the ground robot navigation example) are from the original references. More details (the model description and initial distributions) will be provided in Appendix C.
Training: For each system, we generate 10k trajectories through simulation with varied trajectory lengths from 10 to 100 time-steps, depending on different configurations of the simulation envi- ronment. We use 80% of the samples to train the NN as described in Sec. 3.1 and use the rest for validation. For the trajectory data, we collect the system states and compute for the system diver- gence term. For black-box systems, we use the gradient perturbation method to approximate the derivatives. We use feed-forward NN which has 3 hidden layers with 64 hidden units in each layer. We use PyTorch [71] to train the NN and the training takes 1~2 hours on an RTX2080 Ti GPU.
4.2 Density estimation verification
We first test the density estimation accuracy of our learned NN. We compare our approach with other baselines including kernel density estimation (KDE), Sigmoidal Gaussian Process Density (SGPD) [72] and the histogram approach. For each simulation scenario, we first solve Eq. (2) to generate 20k 100k trajectories of (state, density) pairs, and treat this density value the ground truth. For the KDE method, we choose an Epanechnikov kernel. We then measure the KL divergence between the density estimate of each method and the ground truth. As shown in Fig. our approach has consistently outperformed KDE, SGPD and histogram approaches, with the largest reduction of 99.69% in KI divergence when compared with the histogram approach for the Kraichnan-Orsee system, while our method doesn't 't use any ODE generated density values during training. Also in high-dimension systems (dimension 2 7). the histogram approach fails to predict the density due to the curse of the dimensionality, whereas our approach can always predict the density, with a 30.13% to 99 99.87% % decrease KL divergence comparing to KDE. More plots will be given in Appendix D.
4.3 Reachable set distribution analysis
Forward reachable set distribution analysis. Being confident that our approach is able to pro- vide accurate state density estimation, we extend our NN to do distribution analysis. which is valuable technique in safety-related applications like autonomous driving. Here we use an existing reachability tool RPM [21] to compute the forward reachable sets with probability bounds. Details about how to derive the density and probability bound for the reachable sets are presented in Sec. 3.2 and in Appendix B. Also, RPM was only able to parse the NN for Van der Pol, Double integrator, ground robot navigation, and the FACTEST car model. It fails in handling other high-dimension complex systems due numerical issues when partitioning for the input set. Thus, only report the results on those 4 models in this section. The main purpose is to show that for some systems the
tion based on the set presentation of the tool. Developing better NN reachability tool that can scale to higher dimensional systems is another topic and is out of the scope of paper.

---

CateRed PUZZA
SAPT MARTH
SATURRC PRIEE
Sanpled POWEE
Dryol
Dryak. V-
V-3.6K
5.30X
(500) Ved 172
090 Vad 600
090 Mer
25 50
5.0 -25
-25 25
-5.0 -2.5
(a) The system forward reach set computed by the ConvexHull method((buu) 19), DryVR(red) [68] and GSC[green) [73]. further show the relative volume of each kind reachable set, comparing to the vol- ume the convex hull of the sampled points (gray). As time evolves, the conventional reachability methods often result in a larger over-approximation of the reachable sets with an increasing estimation error.
Oura 54X (pmó 50)
vuts 000005-00.000
Ours voto 48x 0.85.00
voló 44X -6. -0,80)
616
66
100
000
500
12
0.0
-50 -23
00 25
30
30 25 60 25 32
30 -25
to 25 50
-50 -25
20 25 3M
(b) The system forward reach set derived by RPM. The reachable sets are represented in polyhedral cells. The color ranged from dark purple to light yellow indicates the density inside the polyhedral cells. The edges colored in green indicate the boundaries of the RPM polyhedral cells with density below threshold.
Figure 3: The Van der Pol oscillator forward reachable set comparison between (a) the worst-case reachability methods [68, 73] and (b) our probabilistic approach. Our approach clearly identifies reachable sets in high and low densities, and shows that as time evolves, the states will concentrate limit cycle, which expected.
Dryli
10
10
10
20
36
0.00
6.00
0.00
Timastap
Timedeep
Timestep
Van der Pol system
(b) Double integrator
(c) Ground rohot model
(d) FACTEST model
Figure 4: The (relative) volume of reachable set with different probability level using our method, and com- parison other methods. The top part of each figure is in logarithm scale.
density tends to concentrate on certain states, where a small portion of the reachable sets contains the majority of states that are more likely to be reached. Therefore, method can better quantify risks than worst-case reachability, by providing a flexible threshold for the probability of reachability.
We start with the Van der Pol oscillator. The initial states are uniformly sampled from a square region: 2-5,5, 2.3] x 2.5, 2.5). As illustrated in Fig. 1, the system states will gradually converge to a limit cycle. Using worst-case reachability analysis for this system will result in very conservative over-approximation, and this over-approximation will propagate over time and lead to increasing conservativeness of the reachable set estimation. As shown in Fig. 3(a), for the worst-case methods like (red)( [68, 74] and GSG(Green) [73], the volume their estimated reachable set relative to the volume of the convex hull of the system states keeps increasing over time, from 1.9670X to 5.3027X for the VR [68] approach, and from 1.7636X to 2. 2,8086X for GSG [73]. However, our method can give the probability bound for every reachable set in the state space as shown in the heatmap in Fig. 3(b), clearly identifying the region around the limit cycle in high density (bright color), and the rest space in low density (dark color).
We can also compute the relative volume of the reachable sets (comparing to the convex hull) pre- serving different levels of reachable probability, whose evolution over time reflects the system's 's tendency for concentration. As shown in Fig. 4, we use the above systems and study the volume of the reachable set with probability threshold 0.50, 0.70, 0.80, 0.90, and 0.99. As expected, the rel- ative volume will increase as probability threshold increases. In all cases, there exist some time instances where small volume of the reachable set actually preserves high probability, which shows
7

---

the state concentration exists in many existing systems. While as shown in Fig. 4, the worst-case reachability tools can only generate one curve which presents the (relative) volume of the reach- able set that covers all possible states. Not surprisingly, the worst-case reachability tools give very conservative results. We believe using our proposed method to do reachable set distribution analy- sis will benefit future study for systems with uncertainty, and for systems where the failure case is inevitable but happens with a low probability. More comparisons with state-of-the-a reachability methods (Verisig [39]. Sherlock [75] and ReachNN [40]) are shown in Appendix H.
Probability of collssion
Online safety verification under different ini- tial state distributions. Since our approach learns the density concentration function in- stead of absolute density value, it has the flex- ibility to estimate online reach sets distribution with any possible (bounded support) initial dis- tributions. Consider the safety verification for the ground robot navigation problem (shown in Fig. 5(a)): The probability of colliding with a obstacle in the center of the map is deter- mined by the initial state distribution 6, which is parametrized as a truncated Gaussian distri- bution N NW, o² where il and a measure the expectation and uncertainty of the initial robot state. Our method can estimate the upper and lower bounds for the probability of colliding with the obstacle, and this safety evaluation process can run in faster than 50Hz with parallel computation and heuristic searching used (see details in Appendix E). As shown in Fig. 5(b), when the initial state uncertainty decreases from 1.0 to 0.02, the up- per and lower bounds for the probability of collision decrease to close zero (from 0.04236 to 3.5344 x 10-06, where other worst-case reachability analysis methods can only report a collision is inevitable, without quantifying the corresponding risks. This also shows the advantage of our ap- proach in adapting to different density conditions in computation without retraining or fine-tuning.
0.04
Obstacle
Goal
pmax
Start
pmin
0.0
0.2
0.4
0.6
0.8
1.0
Covariance
(a)
(b)
Figure 5: (a) A robot tries to reach the goal while avoiding the obstacle. (b) Probabilistic safety verifica- tion under different initial distributions. As the robot has less uncertainty in its initial position and velocity, it will have lower probability to collide.
Limitations and trade-offs of performing reachable set distribution analysis. Experimental results show that our method can compute much less conservative probabilistic reachable sets than most worst-case reachability methods. This less conservative result benefits from RPM which can provide exact NN reachability analysis by sacrificing scalability. Therefore RPM also constrains us from performing online reachability analysis for high-dimensional systems or systems with large ini- tial sets. Technically, we are solving a harder problem than worst-case reachability approximation, as we need not only the reachable set, but also the density over those reachable states. Like worst- case analysis, this is the reason why it can only scale to initial sets when we want to perform accurate reachable computation. We believe that our approach can better quantify risks under different conditions, especially when unsafe is inevitable (similar to Fig. 5(a)). Our method can also give worst-case reachability by taking all output reachable cells produced by RPM regardless of their density. This worst-case reachability using RPM is less con- servative than other NN reachability tools, at the cost of not scaling to nigh-dimensional systems.
systems and smaller
5 Conclusion and discussion
In this paper, we propose a Neural Network (NN)-based probabilistic safety verification frame- work that can estimate state density, compute reachable sets and corresponding probability. Our Liouville-based NN can accurately estimate the state density even for high-dimension systems. Our probabilistic reachable set framework can handle nonlinear (and potentially black-box) systems with varying initial state distributions and can be used for fast online safety verification. We recognize that the task of computing probabilistic reachable sets is very useful, and our method is more help- ful than worst-case reachability particularly when the system states are more likely to concentrate. One limitation of our approach is that the NN reachability tool we used (RPM) cannot handle high- dimension systems or cases where the initial set is very large, due to the numerical issues when partitioning for polyhedral cells. This limitation is due to the scalability and accuracy trf of NN reachability, which is an independent problem from our paper. We plan to explore other NN reachability methods and more complicated hybrid systems in real-world applications.
6 'How to compute the probability of a reachable set is discussed in Sec.3.2
8

---

Acknowledgments
The NASA University Leadership initiative (grant #80NSSC20M0163) and Ford Motor Company provided funds to assist the authors with their research, but this article solely reflects the opinions and conclusions of its authors and not any NASA or Ford entity,
References
[1] M. Chen and C. J. Tomlin. Hamilton-accbii reachability: Some recent theoretical advances and applications in unmanned airspace management. Annual Review of Control, Robotics, and Autonomous Systems, 1:333-358, 2018.
[2] A. Devonport and M. Arcak. Data-drived reachable set computation using adaptive gaussian process classification and monte carlo methods. In 2020 American Control Conference (ACC), pages 2629-266222222222222222222222222 2634. IEEE, 2020.
[3] X. Chen, E. Ábrahám, and S. Sankaranarayanan. Flow*; An analyzer for non- non-linear hybrid sys- tems. In International Conference on Computer Aided Verification, pages 258-263. Springer, 2013.
[4] H.-D. Tran, X. Yang, D. M. Lopez, P. Musau, L. V. Nguyen, W. Xiang, S. Bak, and T. T. Johnson. Nnv: The neural network verification tool for deep neural networks and learning- enabled cyber-physica systems. In International Conference on Computer Aided Verification, pages 3-17. Springer, 2020.
[5] S. Bansal and C. Tomlin. Deepreach: A deep learning approach to high-dimensional reacha- bility. arXiv preprint (rXiv:2011.02085 2020.
[6] C. Liu, T. Arnon, C. Lazarus, C. Strong, C. Barrett, and M. J. Kochenderfer. Algorithms for verifying deep neural networks. arXiv preprint arXiv: /903.0675, 2019.
[7] A. Devonport and M. Arcak. Estimating reachable sets with scenario optimization. Learning for Dynamics and Control, pages 75-84, PMLR, 2020.
[8] L.. Liehenwein, C. Baykal, I. Gilitschenski, S. Karaman, and Rus. Sampling-based approx- imation algorithms for reachability analysis with provable guarantees. RSS, 2018.
[9] Lew and M. Pavone. Sampling-based reachability analysis: random set theory approach with adversarial sampling. arXiv preprint arXiv:2008. 10180, 2020.
[10] A. B. Kurzhanski and P. Varaiya. Ellipsoidal techniques for reachability analysis. In Inter- national Workshop on Hybrid Systems: Computation and Control, pages 202-214. Springer, 2000.
[11] A. Girard. Reachability of uncertain linear systems using zonotopes. In International Work- shop on Hybrid Systems: Computation and Control, pages 291-305. Springer, 2005.
[12] S. Duggirala and M. Viswanathan. Parsimonious, simulation based verification of linear sys- tems. In International Conference on Computer Aided Verification, pages 477-494. Springer, 2016.
[13] P.J. Meyer, A. Devonport, and M. Arcak. Tira: Toolbox for interval reachability analysis. In Proceedings of the 22nd ACM International Conference on Hybrid Systems: Computation and Control, pages 229, 2019.
[14] Bak. Reducing the wrapping effect flowpipe construction using pseudo-invariant. In Proceedings the 4th ACM SIGBED International Workshop on Design, Modeling, and Eval- uation of Cyber-Physien Systems, pages 40-43, 2014.
[15] M. Althoff, O. Stursberg, and M. Buss. Reachability analysis of nonlinear systems with uncer- tain parameters using conservative linearization. In 2008 47th IEEE Conference on Decision and Control, pages
IEEF, 2008

---

[16] J. Cyranka, M. A. Islam. G. Byrne, P. Jones, S. A. Smolka, and R. Grosu. Lagrangian reacha- bililty. In International Conference on Computer Aided Verification, pages 379-400. Springer, 2017.
[17] M. Ehrendorfer. The liouville equation and prediction of forecast skill. In Predictability and Nonlinear Modelling in Natural Sciences and Economics, pages 29 44. Springer, 1994.
[18] T. Zikmerr,-Zimmerer D. Venturi, Q. Gong, and W. Kang. Density propagation with characteristics-basee deep learning. arXiv preprint arXiv: 1911,09311 2019.
[19] Y. Chen, M. Ahmadi, and A. D. Ames. Optimal safe controller synthesis: density function approach. In 2020 American Control Conference (ACC), pages 5407-5412. IEEE, 2020.
[20] M. Raissi, P. Perdikaris, and G. E. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial dif- ferential equations. Journal of Computational Physics, 378:686-707, 2019.
[21] J. A. Vincent and M. Schwager. Reachable polyhedral marching (rpm): A safety verifica- tion algorithm for robotic systems with deep neural network components. arXiv preprint arXiv:2011. 1609, 2020.
[22] International competition on verifying continuous and hybrid systems. https://cps-vo. Accessed: 2021-06-18 18.
[23] G. Agha and K. Palmskog. A survey of statistical model checking. ACM Transactions on Modeling and Computer Simulation (TOMACS), 28(1):1-39, 2018.
[24] I. M. Mitchell, A. M. Bayen, and C. J. Tomlin. A time-dependent hamilton-jacobi formulation of reachable sets for continuous dynamic games. IEEE Transactions automatic control, 50 (7):947-955 01447-97,,22 2005.
[25] B. Xue, M. Zhang, A. Faswaran, and Q. Li. PAC model checking of black-box continuous- time dynamical systems. IEEE Transactions Computer-Aided Design of Integrated Circuits and Systems, 39(11):3944339) 2020.
[26] A. Berndt, A. Alanwar. K. H. Johansson, and H. Sandberg. Data-driven set-based estimation using matrix zonotopes with set containment guarantees. arXiv preprint arXiv:2101. 10784, 2021.
[27] R. E. Allen, A. A. Clark, J. A. Starek, and M. Pavone. A machine learning approach for real- time reachability analysis. In 2014 IEEE/RSJ international conference on intelligent robots IEEE, 2014.
and systems, pages
[28] M. Rasmussen, J. Rieger, and K. N. Webster. Approximation of reachable sets using optimal control and support vector machines. Journal of Computational and Applied Mathematics, 311: 111:68-33 2017.
[29] A. Thorpe, K. R. Ortiz, and M. M. Oishi. Data-driven stochastic reachability using hilbert space embeddings. arXiv preprint arXiv:2010.08036 2020.
[30] A. Chakrabarty, C. Danielson, S. Di Cairano, and A. Raghunathan. Active learning for estimat- ing reachable sets for systems with unknown dynamics. IEEE Transactions on Cybernetics, 2020.
[31] D. Fridovich-Keil, A. Bajesy, J. F. Fisac, S. L. Herbert, S. Wang, A. D. Dragan, and J. Tom- lin. Confidence-aware motion prediction for real-time collision avoidance The International Journal of Robotics Research, 39(2-3):250-22)
2020.
[32] A. Majumdar, R. Vasudevan, M. M. Tobenkin, and R. Tedrake. Convex optimization of non- linear feedback controllers via occupation measures. The International Journal of Robotics Research, 3309):120- 2014.
[33] A. Abate. Probabilistie reachability for stochastic hybrid systems: theory, computations, and applications. University of California, Berkeley, 2007.
10

---

[34] A. R. R. Matavalam, U. Vaidya, and V. Ajjarapu. Data-driven approach for uncertainty propa- gation and reachability analysis in dynamical systems. In 2020 American Control Conference (ACC), pages 1993-3398. IEEE, 2020.
[35] G. Katz, C. Barrett, D. L. Dill, K. Julian, and M. J. Kochenderfer. Reluplex: An efficient smt solver for verifying deep neural networks. In International Conference on Computer Aided Verification, pages 97-117. Springer, 2017.
[36] G. Katz, D. A. Huang, D. Ibeling, K. Julian, C. Lazarus, R. Lim, P. Shah, S. Thakoor, H. Wu, A. Zeljić, et al. The marabou framework for verification and analysis of deep neural networks. In International Conference Computer Aided Verification, pages 443-452. Springer, 2019.
[37] W. Xiang, H.-D. Tran, and T. T. Johnson. Output reachable set estimation and verification for multilayer neural networks. IEEE transactions on neural networks and learning systems, 29 (11):5777-5783) 2018
[38] X. Yang, H.-D. Tran, W. Xiang, and T. Johnson. Reachability analysis for feed forwad neural networks using face lattices. arXiv preprint arXiv:2003.01226, 2020.
[39] R. Ivanov, J. Weimer, R. Alur, G. J. Pappas, and I. Lee. Verisig: verifying safety properties of hybrid systems with neural network controllers. In Proceedings of the 22nd ACM International Conference on Hybrid Systems: Computation and Control, pages 169-178, 2019.
[40] J. Fan, C. Huang, X. Chen, W. Li, and Q. Zhu. Reachnn*: A tool for reachability analysis of neural-network controlled systems. In International Symposium on Automated Technology for Verification and Analysis, pages 537-542. Springer, 2020.
[41] H. Hu, M. Fazlyab, M. Morari, and G. J. Pappas. Reach-sdp: Reachability analysis of closed- loop systems with neural network controllers via semidefinite programming. 2020 59th IEEE Conference on Decision and Control (CDC), pages 5929-59999 IEEE, 2020.
[42] M. Everett, G. Habibi, and J. P. How. Efficient reachability analysis of closed-loop systems with neural network controllers. arXiv preprint arXiv:2101.01815, 2021.
[43] D. A. Spencer and R. Braun. Mars pathfinder atmospheric entry-trajectory design and dispersion analysis. Journal of Spacecraft and Rockets, 33(5):670-67) 1996.
[44] H. Niederreiter. Random number generation and quasi-Monte Carlo methods. SIAM, 1992.
[45] C. Pantano and B. Shotorban. Least-squares dynamic approximation method for evolution of uncertainty in initial conditions of dynamical systems. Physical Review E, 76(6):66775, 2007.
[46] H. Lee and I. S. Kang. Neural algorithm for solving differential equations. Journal of Compu- tational Physics, 91(1): 10-131, 11011110-13,, 1990.
[47] I. E. Lagaris, A. Likas, and D. Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. IEEE transactions on neural networks, 95))987-11000, 1998.
[48] T. Uchiyama and N. Sonehara. Solving inverse problems in nonlinear pdes by recurrent neural networks. In IEEE International Conference on Neural Networks, pages 99-102. IEEE, 1993.
[49] Sirignano and K. Spiliopoulos. Dgm: A deep learning algorithm for solving partial differ- ential equations, Journal of computational physics, 375: 775:1339-1364 2018
[50] J. Berg and K. Nyström. A unified deep artificial neural network approach to partial differential equations in complex geometries. Neurocomputing. 317:28-44 41, 2018.
[51] E. Weinan and B. Yu. The deep ritz method: deep baaming-based numerical algorithm for solving variational problems. Communications in Mathematics and Staristics, 6(1):1-12, 2018.
[52] J. He, L. J. Xu, and C. Zheng. Relu deep neural networks and linear finite elements. arXiv preprint arXiv: 1rXXiv1807.005 03973, 2018.
[53] Han, A. Jentzen, al. Algorithms for solving high dimensional pdes: From nonlinear monte carlo to machine learning, arXiv preprint arXiv:2088.1 13333, 2020.

---

[54] Y. Khoo, J. Lu, and L. Ying. Solving parametric pde problems with artificial neural networks European Journal of Applied Mathematics, 32(3):421-433 2021.
[55] G. Kutyniok, Petersen, M. Raslan, and R. Schneider. A theoretical analysis of deep neural networks and parametric pdes. Constructive Approximation, pages 1-53, 2021.
[56] D. Zhang, L. Lu, L. Guo, and G. E. Karniadakis. Quantifying total uncertainty in physics- informed neural networks for solving forward and inverse stochastic problems. Journal of Computational Physics, 397:10890 2019.
[57] L.. Yang, D. Zhang, and G. E. Karniadakis. Physics-informed generative adversarial networks for stochastic differential equations. arXiv preprint arXiv:1811.0203: 2018.
[58) Z. Li, N. Kovachki, K. Arizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anand- kumar. Fourier neural operator for parametric partial differential equations. arXiv preprint arXiv:2010.08897 2020.
[59] Z. Long, Y. Lu, X. Ma, and B. Dong. Pde-net: Learning pdes from data. In International Conference Machine Learning, pages 3208-3216. PMLR, 2018.
[60] L. Lu, P. Jin, and G. E. Karniadakis. Deeponet: Learning nonlinear operators for identify- ing differential equations based on the universal approximation theorem of operators. arXiv preprint arXiv:1910.033
2019.
[61] A. Koryagin, R. Khudorozkov, and S. Tsimfer. Pydens: A python framework for solving differential equations with neural networks. arXiv preprint
2019.
[62] L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis. Deepxde: A deep learning library for solving differential equations. SIAM Review, 63(1):208-228, 2021.
[63] F. Chen, D. Sondak, P. Protopapas, M. Mattheakis, S. Liu, D. Agarwal, and M. Di Giovanni. Neurodiffeq: A python package for solving differential equations with neural networks. Jour- nal of Open Source Software, 5(46): 1931, 2020.
[64] V. Nair and G. E. Hinton. Rectified linear units improve restricted boltzmann machines. In leml, 2010.
[65] T. Chai and R. R. Draxler. Root mean square error (rmse) or mean absolute error (mac)?- arguments against avoiding rmse in the literature. Geoscientific model development, 7(3): 1247-1250, 2014.
[66] S. A. Orszag and L. Bissonnette. Dynamical properties of truncated wiener-hermite expan- sions. The Physics of Fluids, 10( (0(12):20032261, 1967.
[67] Y.-C. Chang, Roohi, and S. Gao. Neural Yyарunov control. arXiv:2005.00611 2020.
arXiv preprint
[68] C. Fan, K. Miller, and S. Mitra. Fast and guaranteed safe controller synthesis for nonlinear vehicle models. In International Conference on Computer Aided Verification, pages 629-652. Springer, 2020.
[69] Heidlauf, A. Collins, M. Bolender, and Bak. Verification challenges in f-16 ground collision avoidance and other automated maneuvers. In ARCH@ ADHS, pages 208-217 2018.
[70] H. Zhu, Z. Xiong, S. Magill, and S. Jagannathan. An inductive synthesis framework for ver- ifiable reinforcement learning. In Proceedings of the 40th ACM SIGPLAN Conference on Programming Language Design and Implementation, pages 686-701, 2019.
[71] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, et al. Pytorch: An imperative style, high-performance deep learning library. arXiv preprint arXiv:1912.0177 2019.
[72] C. Donner and Opper. Efficient bayesian inference of sigmoidal gaussian cox processes. 10. 10144279//epooiinne--8888 2018
12

---

[73] M. Everett, G. Habibi, and J. P. How. Robustness analysis of neural networks via efficient partitioning with applications in control systems. IEEE Control Systems Letters, 2020.
[74] Du, Z. Huang, Т. Liu, T. Ji, K. Xu, Q. Gao, Sibai, K. Driggs-Campbol and S. Mitra. Online monitoring for safe vechice interactions. In 2020 IEEE 23rd International Conference on Intelligent Transportation Systems (ITSC), pages 1-8. IEEE, 2020.
[75] S. Dutta, S. Jha, S. Sanakaranarayanan, and A. Tiwari. Output range analysis for deep neural networks, arXiv preprint arXiv: 1700:11709..991300 2017.
13


Reviews:
[{"review_id": "kpUS3qzke5z", "reviewer": "Reviewer_bVFw", "summary": "This paper presents a data-driven method to approximate the probability density of reachable states for a nonlinear system. The proposed method learns the state density concentration function and the flow map (forward dynamics model that predicts the state at t given initial state x_0) using a neural network. The learned models are used to predict reachable states and corresponding density distribution. \nThe model is kind of a physics-informed NN, which is constructed upon the Liouville theorem. The proposed method formulates learning objectives (or apply constraints) based on the theorem that asserts that the state density follows the Liouville PDE (Eq. 1). The theorem characterizes the evolution of the state density along the trajectory.\n\n\nThe paper presents three experiments. \nThe training data is generated by simulating the benchmark systems. and evaluation is done by the density values obtained by solving the ODE (Eq. 2).\n\n(1) The proposed method is evaluated 10 different systems including low-dimensional academic models such as Van der Pol oscillator and high-dimensional systems. In all examples, the proposed method could predict the state density better than the traditional KDE method in terms of the KL divergence to the ground truth. \n(2) The proposed model is used to compute reachable sets of some of the benchmark systems. RPM is used to represent & compute reachable regions. In all the experiments, the proposed method could predict the forward reachable set close to the ground truth while the existing worst-case methods often fail in the systems where the states concentrate in small regions.\n(3) Online safety verification experiment of the proposed method. Here the authors show that the reachability analysis can be done efficiently (20 ms per step) online using the learned model given arbitrary initial states.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "- The problem definition is clear and well-presented. It was easy to understand what problem this work is addressing.\n- The model is constructed such that the density concentration function (G) is independent of the initial state, which can be used to quickly compute the state density from any given initial density. This smart formulation enables the online reachability analysis shown in 4.3.\n- The benefit of predicting the state density is that we can leverage different density conditions as demonstrated in supplementary E. However, this advantage is not clearly mentioned in the main text.\n- The online reachability analysis has the potential for building a practical real-time control framework for safety-critical systems. Unfortunately, this paper does not provide a scalable solution for computing the reachability set.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "I enjoyed reading this paper. The problem is well defined and the approach is very logical. \nHowever, the experimental results and analyses are weak. The density values computed by solving Liouville PDE are used as ground truth to compare the proposed method against existing methods, which seems to be an unfair comparison because the proposed model is constructed to follow the Liouville PDE. The computation of the reachability sets presented in the last part is very interesting, but it is limited to low-dimensional systems.", "technical_quality": "Very Good:  Convincing.", "clarity_of_presentation": "Good: Minor revision needed to improve organization or clarity.", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"impact": "3: The work is incremental and unlikely to have much impact, even though it may be technically correct and well executed.", "reviewer_expertise": "Fair: Some knowledge of the area"}, "extra_text": {"issues": "- In the first experiment, the ground truth labels are generated by solving the Liouville PDE. The fairness of this comparison is questionable. The proposed method is constructed to follow the Liouville PDE while the baselines (KDE & histogram) are not.  The experiment seems to be designed such that the proposed method outperforms others.\n- It is unclear how the initial density is defined for each experiment.\n- In the conclusion, there is a statement saying that the reachability analysis can be scaled up by taking other tools. I wonder If there is a more scalable method than the RPM that the authors are considering.\n- Line 32, The author claims that \"almost the entire space in the figure is reachable\" if an existing reachability algorithm is used. This claim contradicts the result in Figure 2 (1). It is true for GSG or DryVR, but it seems like the Convex Hull method [9] does not exhibit such behavior.\n- In Figure 2 and Figure 37-48 in supplementary, the gray points are not given in the legend. It seems like they are Monte Carlo samples that approximate the ground truth.\n- The result & takeaway of the last experiment is not clearly presented. I had to look up the supplementary to understand the experiment."}}, {"review_id": "bRLi6HkmQ8s", "reviewer": "Reviewer_mxEJ", "summary": "This paper addresses the problem of estimating the density\ndistribution of reachable states of autonomous systems. The authors\ntrain two neural networks (NNs): one to approximate the system\ndynamics and one to approximate the density evolution over time. Once\nthe NNs are trained, a reachability tool, RPM, is used to provide\ndensity estimates over sets of initial conditions and\ndensities. Evaluation is provided over 10 benchmarks selected from the\nliterature.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The idea of estimating density distributions, as opposed to standard\nreachable sets, is interesting and relevant, especially in the context\nof evaluating the risk associated with a certain task. The paper is\nwell written and the problem is clearly defined.\n\nMy main comment is that it is not clear how the proposed method\ncompares with state-of-the-art worst-case reachability tools such as\nVerisig, Sherlock, ReachNN, etc. While I agree that estimating the\ndensity distribution would be very useful, this problem is strictly\nharder than estimating the reachable sets without taking into account\nthe probability distribution. In particular, the authors discuss\nnumeric issues as a challenge to their approach but I would argue that\nthis is not a numeric issue as much as it is an issue of whether to\ninclude all reachable states, even those that have a very low\nprobability of being reached. On the other hand, if states with very\nlow probability are being ignored, then this needs to be clearly\nstated in the problem statement and in the solution.\n\nRelated to the above, another interesting comparison would be with a\nsampling-based method, which should perform well in low-dimensional\nsystems. In particular, given a sufficient number of simulations, one can\nobtain a good estimate of the density distribution that can be\ncompared with the one produced in the paper.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "I am recommending a weak reject because I think there is a mismatch between the author claims and the actual contributions. In particular, the authors seem to imply that the proposed method is more efficient than worst-case reachability methods (\"...can be used for fast online safety verification\"). I think the opposite is true, as discussed in my detailed comments -- this problem is strictly harder unless there are simplifying assumptions (which are not clearly stated).", "technical_quality": "Good: Generally convincing, but some additional experiments or analyses needed.", "clarity_of_presentation": "Good: Minor revision needed to improve organization or clarity.", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"impact": "2: Unsure, difficult to assess.", "reviewer_expertise": "Excellent: Expert knowledge on the topic of the paper"}, "extra_text": {"issues": "I would like to see a discussion on the limitations of the proposed approach as compared with existing worst-case reachability methods. Furthermore, since neural networks are used to approximate the system dynamics and density, I am not sure how the reachability problem changes since it's not clear that the neural network can well approximate low-probability spaces."}}, {"review_id": "AwqG4IgDWu5", "reviewer": "Reviewer_HgHK", "summary": "In this paper the authors present a method for approximating the probabilistic forward evolution of an autonomous dynamical system by finding an approximate solution to the Liouville PDE. They train a fully connected neural network which learns the flowmap (i.e. trajectory) and the density function. The loss function is composed of two terms: prediction error for the dynamics, and PDE error. Moving from single-point estimation to set-based estimation, the authors additionally use RPM, a method which propagates the polyhedral partitions of the input space of a neural network as an affine transformation, in order to compute the forward reachable set and it's associated probability density. They showcase their algorithm on a variety of examples with improved KL divergence results than other methods.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n* It deals with an important problem in safety by analyzing probabilities rather than yes or no safety queries.\n* The paper is clear and easy to read.\n* They demonstrate the framework on a large variety of experiments.\n* The results indicate that their approach obtains better KL performance than other methods. \n\nWeaknesses: \n* While RPM is a useful tool for verification, it is only used for a small number of examples in the paper given that polyhedral propagation does not scale well with input dimension. I would have preferred to see some of the information in the appendix replacing section 3.2.\n* Some explanation of why they jointly model the flow and density would have been useful.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "I think this paper proposes a method which is useful for safety applications. Taking a probabilistic perspective rather than a binary safe/unsafe approach can be of good use in robotics.", "technical_quality": "Good: Generally convincing, but some additional experiments or analyses needed.", "clarity_of_presentation": "Very Good:  Well organized and clear.", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"impact": "4: The work contains interesting new ideas that will potentially have major impact in robotics or machine learning.", "reviewer_expertise": "Good: General knowledge of the area"}, "extra_text": {"issues": "I would like these questions/issues to be addressed:\n1.) Why are there two sets of parameters omega and theta different when only one single neural network is trained? Why not use two separate networks.\n2.) There should be some more in-depth discussion of how the dynamics and dimensionality of the dynamical systems affect the number of samples required to closely approximate the Liouville PDE and the flow."}}]

Author rebuttal:

Meta-review: ### Final Meta-Review

After the revision and discussion process, reviewers agree that the authors have addressed their concerns to a satisfactory degree, both through clarifications during the discussion and through improvements to the paper's rigor and completeness. The final manuscript is significantly improved, and I can confidently recommend its inclusion in the CoRL technical program.


### Original Meta-Review

Reviewers agree that the paper presents interesting ideas around the use of density distributions to provide a richer risk quantification relative to the more commonly used reachable sets. However, they also express various concerns regarding the strength of the results and the accuracy of contribution claims.

In particular, Reviewer mxEJ points out a “mismatch between the author claims and the actual contributions”. This Area Chair agrees that statements about the better scalability and accuracy of state density estimation relative to worst-case (nondeterministic) analysis are not substantiated by the results. Reviewers HgHK and bVFw also express concerns related to scalability with state dimension not being clearly established. In addition, Reviewer bVFw points out multiple issues with clarity and rigor in the results and conclusions.

Adding to the above Reviewer concerns, the paper’s claim that this work is “the first to characterize and compute the reachability density/probability distribution for safety analysis” is incorrect as stated. There are in fact multiple works, in robotics alone, that compute probabilistic state distributions: for example, the work on confidence-aware predictions from Claire Tomlin and Anca Dragan’s groups, and on occupation measures from Ram Vasudevan’s group (also based on the Liouville equation).

The authors need to qualify the novelty of their contribution with respect to prior work, more rigorously state—and support—the advantages of their proposed approach, and adequately address the issues raised around the reported experimental results.
Decision: Accept (Poster)