Model-based Validation as Probabilistic Inference

We build upon the notation of Corso et al. (2021) to describe the safety validation problem. Consider a system under test that takes actions in an environment after receiving observations. We define a safety property $\psi$ that we aim to evaluate for the system under test. The property $\psi$ is defined over the state trajectories of the environment, $\mathbf{s}=[s_{1},\dots,s_{t}]$, where $s_{t}$ is the state of the environment at time $t$. If the trajectory results in a failure, we say that $\mathbf{s}\notin\psi$.

We perturb the environment with disturbances with the goal of inducing the system under test to violate the safety specification. A disturbance trajectory $\mathbf{x}$ has probability density $p(\mathbf{x})$. We assume that disturbances determine all sources of stochasticity in the environment. For example, sensor noise and stochastic dynamics are potential disturbances that could lead to failure.

We denote the simulation of the system under test with environment disturbances by a dynamics function $f$. Note that $f$ represents the dynamics of both the system under test and the environment. The resulting state trajectory under disturbances is written as $\mathbf{s}=f(\mathbf{x})$. Since all stochasticity in the simulation is determined by $\mathbf{x}$, the disturbance model $p(\mathbf{x})$ induces a distribution over state trajectories, $p(\mathbf{s})$. We seek the distribution over state trajectories that violate the safety property, the conditional distribution $p(\mathbf{s}\mid\mathbf{s}\notin\psi)$.

Environment dynamics or characteristics of the system under test often make it impossible to express the failure distribution analytically. However, it is usually easy to sample from these systems by simulating the system under test with disturbances. Bayesian inference methods such as MCMC approximate a target distribution by drawing samples. In this work, we frame the safety validation problem as Bayesian inference of the distribution over failures.

Sampling-based Bayesian inference requires a way to draw samples from a target distribution. We represent the state trajectory distribution as a probabilistic model using probabilistic programming (Gordon et al., 2014). Under this paradigm, we can describe the distribution over state trajectories as a function of random variables that represent the disturbance model. Modern probabilistic programming languages, such as Turing (Ge et al., 2018), let us wrap existing simulations in a probabilistic program with little modification of the underlying simulation.

This probabilistic model represents the distribution over state trajectories induced by the disturbance model. We represent the distribution over failures by conditioning the probabilistic model on the desired outcome of the simulation — that the simulated trajectory is a failure. Any safe trajectory has zero likelihood under the conditional distribution. We represent this condition using the Dirac delta distribution, whose value is zero for all state trajectories except when $\mathbf{s}\notin\psi$. Under this formulation, the log-likelihood of a state trajectory $\ell(\mathbf{s}\mid\mathbf{s}\notin\psi)$ is

The main challenges associated with sampling from $p(\mathbf{s}\mid\mathbf{s}\notin\psi)$ are the distribution’s high dimensionality, discontinuity, and multimodality. The following sections show how we address each of these challenges.

State trajectories in sequential systems of interest tend to be high dimensional. It is well established that black-box, gradient-free sampling from high-dimensional distributions becomes difficult with increasing dimensionality (Neal, 2011). Gradient-based sampling algorithms such as Hamiltonian Monte Carlo (HMC) (Duane et al., 1987; Girolami and Calderhead, 2011) use information about the shape of the underlying distribution to take large steps in parameter space, allowing them to scale to higher dimensions.

To address the challenge of sampling high-dimensional failure trajectories in sequential systems, we take a model-based approach. Gradients of the dynamics function $f$ with respect to disturbances are computed using automatic differentiation. The simulation environment and system under test are modeled in a framework compatible with automatic differentiation (AD). AD enables algorithms to compute exact gradients of the posterior likelihood density with respect to the input disturbances. We use Zygote (Innes et al., 2019), which performs source-to-source AD of Julia code.

In contrast with some previous approaches that assume the system is Markov (Corso et al., 2020), our approach can handle non-Markov systems. We estimate the distribution over failures using rollouts of the system dynamics under disturbances, which fully characterize the system behavior and allows us to validate policies that depend on state estimates rather than true states.

Gradients of the posterior density will be undefined in regions where the state trajectory is safe due to the Dirac delta in Eq. 1. Additionally, the boundary between safe and failing state trajectories is discontinuous. As a result, gradient-based sampling algorithms have no information about how to improve a proposal. We relax the discontinuity by approximating the Dirac delta distribution as a one-dimensional Gaussian with a small variance $\epsilon$. We also assume that we can calculate a distance to failure $\Delta(\mathbf{s})$ for any safe trajectory. The likelihood of the distance to failure under the Gaussian approximation of the discontinuity is added to the likelihood of the sampled trajectory. The log-likelihood of a trajectory $\ell(\mathbf{s})$ under this smoothing approximation is

The role of AD is to calculate $\nabla_{\mathbf{x}}\ell(\mathbf{s}\mid\mathbf{s}\notin\psi)$, which is now defined for safe and failure trajectories. The Gaussian smooths the discontinuity between the safe and failure regions of the state trajectory space while changing the density in the failure region by a constant. This smoothing approximation along with AD enables us to use gradient-based sampling methods.

The variance $\epsilon$ acts as a hyperparameter, with smaller values better approximating the discontinuity and larger values yielding smoother posterior densities. Consider a validation problem for illustrative purposes only, where the state distribution is the unit Gaussian and failure occurs when $|s|>5$. Fig. 2 shows the log-likelihood of single-step state trajectories for various values of $\epsilon$. As $\epsilon$ decreases towards zero, the resulting distribution is a better approximation of the true failure distribution, shown in red. Larger values of $\epsilon$ smooth the gap between the separated failure modes on the far left and right tails of the distribution. Empirically, we find values in the range $\epsilon\in[0.01,0.1]$ work well for the problems investigated.

The generality of the approach makes available a wide variety of off-the-shelf Bayesian inference methods. We use HMC with No-U-Turn Sampling (NUTS) (Hoffman et al., 2014), which is a general HMC variant that excels at sampling from complex, high-dimensional posterior distributions.

Multimodality in the distribution over failure trajectories poses a challenge for many MCMC algorithms (Rudoy and Wolfe, 2006). Sampling algorithms may converge to a single mode rather than mixing across the modes in parameter space. In our illustrative example, a single MCMC chain may only sample from the left or right mode of the failure distribution. One way to explore a multimodal space is to run many chains of MCMC from different starting points. In our sampling approach, we use multiple NUTS chains with different starting disturbance trajectories. The number of chains selected depends on the difficulty of the problem, with higher dimension and more modes generally requiring more chains.

The first system we validate is an inverted pendulum with a neural network controller. While this is not a safety-critical system, it demonstrates the ability of the proposed approach to validate policies with machine learning components. In the inverted pendulum, a policy exerts a control torque to keep the pendulum balanced upright. The state of the system is $s=(\theta,\dot{\theta})$ where $\theta$ is the angle of the pendulum from vertical and $\dot{\theta}$ is the angular velocity. The control policy takes continuous-valued actions $a\in[-T_{\max},T_{\max}]$, where $T_{\max}$ is the maximum control effort of the actuator. The pendulum is underactuated, meaning that the available torque is too weak to result in arbitrary state trajectories. If the pendulum reaches beyond about $0.5\text{\,}\mathrm{rad}\text{/}$ from vertical, the available torque will not be able to overcome the torque due to the pendulum’s mass, and the pendulum will fall. The pendulum and failure condition is illustrated in Fig. 3.

We optimize a neural network policy using the Proximal Policy Optimization (Schulman et al., 2017) reinforcement learning algorithm, which has been shown to perform well on the inverted pendulum. The policy is represented by a two-layer neural network with $32$ neurons per layer and tanh activations between hidden layers. The policy is trained for ${10}^{6}\text{\,}\mathrm{t}$ime steps.

We sample failures of control policy by considering external torque disturbances from the environment. Adversarial torque disturbances are assumed to be distributed according to a zero-mean Gaussian with variance $\sigma_{n}$. This experiment uses $\sigma_{n}=$0.1\text{\,}\mathrm{N}\text{\,}\mathrm{m}\text{/}$$ and $T_{\max}=$2\text{\,}\mathrm{N}\text{\,}\mathrm{m}\text{/}$$ so that disturbances are likely to be much smaller than the maximum control torque. A failure occurs when the torque due to gravity overcomes the maximum control torque, which occurs when $|\theta|\geq$0.5\text{\,}\mathrm{rad}\text{/}$$. The distance to failure is the minimum distance from the failure threshold. Simulations span a $5\text{\,}\mathrm{s}\text{/}$ horizon with a $0.1\text{\,}\mathrm{s}\text{/}$ time step. We sample $10000$ samples for each method, with $10$ MCMC chains and $1000$ samples per chain.

The second experiment involves an autonomous vehicle (AV) approaching a single pedestrian at a crosswalk and is based on the work of Koren et al. (2018). This scenario is illustrated in Fig. 3. The system under test is the autonomous vehicle’s control policy, which is defined by the Intelligent Driver Model (IDM) (Treiber et al., 2000). The IDM is a model-based controller that keeps the vehicle in its lane, following the traffic or other obstacles ahead while maintaining a safe distance. The goal of the AV is to identify the pedestrian crossing the road and apply adequate braking to avoid coming too close. Each agent’s state is represented by a 4-tuple $s=(r_{x},r_{y},v_{x},v_{y})$ where $(r_{x},r_{y})$ is the agent’s position and $(v_{x},v_{y})$ is the agent’s velocity. In the noise-free case, the AV comes to a stop and the pedestrian safely crosses.

The adversary has control over pedestrian acceleration and noise in the AV’s sensors. The disturbances to the system are represented by a 6-tuple $x=(\delta r_{x},\delta r_{y},\delta v_{x},\delta v_{y},a_{x},a_{y})$ where $(\delta r_{x},\delta r_{y})$ is the noise in the pedestrian’s observed position, $(\delta v_{x},\delta v_{y})$ is the noise in the pedestrian’s observed velocity, and $(a_{x},a_{y})$ is the pedestrian’s acceleration. The noise components are processed at each time step and added to the vehicle’s observation of each pedestrian’s position and velocity. The disturbances are assumed to be distributed according to a multivariate Gaussian distribution with zero mean and diagonal covariance matrix. The diagonal components of covariance are $\sigma^{2}_{r}=0.1$ for position measurement noise, $\sigma^{2}_{v}=0.1$ for velocity measurement noise, $\sigma^{2}_{lat}=0.01$ over the pedestrian’s lateral acceleration, and $\sigma^{2}_{lon}=0.1$ over the pedestrian’s longitudinal acceleration. These values are chosen to reflect relatively small sensor noise, and that the pedestrian is less likely to accelerate longitudinally while crossing the street.

Failures of this scenario are defined by any collision between the pedestrian and vehicle. The distance to failure is the minimum Euclidean distance between the vehicle and the pedestrian. We sample $20000$ samples for each method, with $20$ MCMC chains and $1000$ samples per chain.

The final validation problem is a simulation of a lunar lander with partially observed state information. The objective of the lunar lander’s policy is to guide the vehicle in a target area with a soft landing. The vehicle state is represented by the 6-tuple $(x,y,\theta,\dot{x},\dot{y},\dot{\theta})$, where $x$ is the horizontal position, $y$ is the vertical position, and $\theta$ is the vehicle’s orientation. The vehicle’s dynamics are deterministic. The vehicle makes noisy observations of its angular rate $\omega$, horizontal speed $v$, and above ground level altitude ($AGL$). $AGL$ is a noisy measurement of the distance from the vehicle’s center of mass to the ground along the vehicle’s longitudinal axis. The continuous actions are defined by $(T,F_{x},d)$ where $T$ is the main thrust acting along the longitudinal axis, $F_{x}$ is lateral corrective thrust acting at an offset of $d$ from the center of mass. The lander is illustrated in Fig. 3.

We optimize a control policy by considering a discrete, fully observable version of the continuous, partially observable problem. We use approximate value iteration to solve for a table-based policy in the discrete problem. Next, we create a differentiable representation of the table-based policy by performing behavior cloning with a neural network approximator. The network is trained to minimize the mean squared error between the continuous output and the table-based action.

We assume that observation noise is identically and independently distributed at each time step according to a zero-mean Gaussian with a diagonal covariance matrix. We use $\sigma^{2}_{\omega}=0.02$, $\sigma^{2}_{v}=0.1$, and $\sigma^{2}_{AGL}=1.0$. An extended Kalman filter maintains a multivariate Gaussian belief over the true state in the partially observable problem. The policy uses the mean to calculate an action.

Our approach samples sequences of observation noise that lead to a hard landing, defined as the landing velocity being greater than $6\text{\,}\mathrm{m}\text{/}\mathrm{s}$. Empirically, we discovered the lunar lander contained multiple isolated failure modes. We select initial conditions for the MCMC samplers by performing black-box optimization of the probabilistic model’s posterior density. Specifically, we use particle swarm optimization with five particles and $100$ iterations to select initial conditions for each MCMC chain. We sample $30$ MCMC chains with $1000$ samples each to explore the different failure modes.