Probabilistic Reach-Avoid for Bayesian Neural Networks
Model-based reinforcement learning seeks to simultaneously learn the dynamics
of an unknown stochastic environment and synthesise an optimal policy for
acting in it. Ensuring the safety and robustness of sequential decisions made
through a policy in such an environment is a key challenge for policies
intended for safety-critical scenarios. In this work, we investigate two
complementary problems: first, computing reach-avoid probabilities for
iterative predictions made with dynamical models, with dynamics described by
Bayesian neural network (BNN); second, synthesising control policies that are
optimal with respect to a given reach-avoid specification (reaching a "target"
state, while avoiding a set of "unsafe" states) and a learned BNN model. Our
solution leverages interval propagation and backward recursion techniques to
compute lower bounds for the probability that a policy's sequence of actions
leads to satisfying the reach-avoid specification. Such computed lower bounds
provide safety certification for the given policy and BNN model. We then
introduce control synthesis algorithms to derive policies maximizing said lower
bounds on the safety probability. We demonstrate the effectiveness of our
method on a series of control benchmarks characterized by learned BNN dynamics
models. On our most challenging benchmark, compared to purely data-driven
policies the optimal synthesis algorithm is able to provide more than a
four-fold increase in the number of certifiable states and more than a
three-fold increase in the average guaranteed reach-avoid probability.
Matthew Wicker, Luca Laurenti, Andrea Patane, Nicola Paoletti, Alessandro Abate, Marta Kwiatkowska
# Probabilistic Reach-Avoid for

###### Abstract

Model-based reinforcement learning seeks to simultaneously learn the dynamics of an unknown stochastic environment and synthesise an optimal policy for acting in it. Ensuring the safety and robustness of sequential decisions made through a policy in such an environment is a key challenge for policies intended for safety-critical scenarios. In this work, we investigate two complementary problems: first, computing reach-avoid probabilities for iterative predictions made with dynamical models, with dynamics described by Bayesian neural network (BNN); second, synthesising control policies that are optimal with respect to a given reach-avoid specification (reaching a "target" state, while avoiding a set of "unsafe" states) and a learned BNN model. Our solution leverages interval propagation and backward recursion techniques to compute lower bounds for the probability that a policy's sequence of actions leads to satisfying the reach-avoid specification. Such computed lower bounds provide safety certification for the given policy and BNN model. We then introduce control synthesis algorithms to derive policies maximizing said lower bounds on the safety probability. We demonstrate the effectiveness of our method on a series of control benchmarks characterized by learned BNN dynamics models. On our most challenging benchmark, compared to purely data-driven policies the optimal synthesis algorithm is able to provide more than a four-fold increase in the number of certifiable states and more than a three-fold increase in the average guaranteed reach-avoid probability.

keywords: Reinforcement Learning, Formal Verification, Certified Control Synthesis, Bayesian Neural Networks, Safety, Reach-while-avoid +
## 1 Introduction

The capacity of deep learning to approximate complex functions makes it particularly attractive for inferring process dynamics in control and reinforcement learning problems (Schrittwieser et al., 2019). In safety-critical scenarios where the environment and system state are only partially known or observable (e.g., a robot with noisy actuators/sensors), Bayesian models have recently been investigated as a safer alternative to standard, deterministic, Neural Networks (NNs): the uncertainty estimates of Bayesian models can be propagated through the system decision pipeline to enable safe decision making despite unknown system conditions (McAllister and Rasmussen, 2017; Carbone et al., 2020; Depeweg et al., 2017). In particular, _Bayesian Neural Networks_ (BNNs) retain the same advantages of NNs (relative to their approximation capabilities) and also enable reasoning about uncertainty in a principled probabilistic manner (Neal, 2012; Murphy, 2012), making them very well-suited to tackle safety-critical problems.

In problems of sequential planning, time-series forecasting, and model-based reinforcement learning, evaluating a model with respect to a control policy (or strategy) requires making several predictions that are mutually dependent across time (Liang, 2005; Deisenroth and Rasmussen, 2011). While multiple models can be learned for each time step, a common setting is for these predictions to be made iteratively by the same machine learning model (Huang and Rosendo, 2020), where the state of the predicted model at each step is a function of the model state at the previous step and possibly of an action (from the policy). We refer to this setting as _iterative predictions_.

Unfortunately, performing iterative predictions with BNN models poses several practical issues. In facts, BNN models output probability distributions, so that at each successive timestep the BNN needs to be evaluated over a probability distribution, rather than a fixed input point - thus posing the problem of successive predictions over a stochastic input. Even when the posterior distribution of the BNN weights is inferred using analytical approximations, the deep and non-linear structure of the network makes the resulting predictive distribution analytically intractable (Neal, 2012). In iterative prediction settings, the problem is compounded and exacerbated by the fact that one would have to evaluate the BNN, sequentially, over a distribution that cannot be computed analytically (Depeweg et al., 2017). Hence, computing sound, formal bounds on the probability of BNN-based iterative predictions remains an open problem. Such bounds would enable one to pro
vide safety guarantees over a given (or learned) control policy, which is a necessary precondition before deploying the policy in a real-world environment (Polymenakos et al., 2020; Vinogradska et al., 2016).

In this paper, we develop a new method for the computation of probabilistic guarantees for iterative predictions with BNNs over _reach-avoid_ specifications. A reach-avoid specification, also known as constrained reachability (Soudjani and Abate, 2013), requires that the trajectories of a dynamical system reach a goal/target region over a given (finite) time horizon, whilst avoiding a given set of states that are deemed "unsafe". Probabilistic reach-avoid is a key property for the formal analysis of stochastic processes (Abate et al., 2008), underpinning richer temporal logic specifications: its computation is the key component for probabilistic model checking algorithms for various temporal logics such as PCTL, csLTL, or BLTL (Kwiatkowska et al., 2007; Cauchi et al., 2019).

Even though the exact computation of reach-avoid probabilities for iterative prediction with BNNs is in general not analytically possible, with our method, we can derive a guaranteed (conservative) lower bound by solving a backward iterative problem obtained via a discretisation of the state space. In particular, starting from the final time step and the goal region, we back-propagate the probability lower bounds for each discretised portion of the state space. This backwards reachability approach leverages recently developed bound propagation techniques for BNNs (Wicker et al., 2020). In addition to providing guarantees for a given policy, we also devise methods to synthesise policies that are maximally certifiable, i.e., that maximize the lower bound of the reach-avoid probability. We first describe a numerical solution that, by using dynamic programming, can synthesize policies that are maximally safe. Then, in order to improve the scalability of our approach, we present a method for synthesizing approximately optimal strategies parametrised as a neural network. While our method does not yet scale to state-of-the-art reinforcement learning environments, we are able to verify and synthesise challenging non-linear control case studies.

We validate the effectiveness of our certification and synthesis algorithms on a series of control benchmarks. Our certification algorithm is able to produce non-trivial safety guarantees for each system that we test. On each proposed benchmark, we also show how our synthesis algorithm results in actions whose safety is significantly more certifiable than policies derived via deep reinforcement learning. Specifically, in a challenging planar navigation benchmark, our synthesis method results in policies whose certified safety
probabilities are eight to nine times higher than those for learned policies.

We further investigate how factors like the choice of approximate inference method, BNN architecture, and training methodology affect the quality of the synthesised policy. In summary, this paper makes the following contributions:

* We show how probabilistic reach-avoid for iterative predictions with BNNs can be formulated as the solution of a backward recursion.
* We present an efficient certification framework that produces a lower bound on probabilistic reach-avoid by relying on convex relaxations of the BNN model and said recursive problem definition.
* We present schemes for deriving a maximally certified policy (i.e., maximizing the lower bound on safety probability) with respect to a BNN and given reach-avoid specification.
* We evaluate our methodology on a set of control case studies to provide guarantees for learned and synthesized policies and conduct an empirical investigation of model-selection choices and their effect on the quality of policies synthesised by our method.

A previous version of this work (Wicker et al., 2021) has been presented at the thirty-seventh Conference on Uncertainty in Artificial Intelligence. Compared to the conference paper, in this work, we introduce several new contributions. Specifically, compared to Wicker et al. (2021) we present novel algorithms for the synthesis of control strategies based on both a numerical method and a neural network-based approach. Moreover, the experimental evaluation has been consistently extended, to include, among others, an analysis of the role of approximate inference and NN architecture on safety certification and synthesis, as well as an in-depth analysis of the scalability of our methods. Further discussion of related works can be found in Section 7.

## 2 Bayesian Neural Networks

In this work, we consider fully-connected neural network (NN) architectures \(f^{w}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}\) parametrised by a vector \(w\in\mathbb{R}^{n_{w}}\) containing all the weights and biases of the network. Given a NN \(f^{w}\) composed by \(L\) layers, we denote by \(f^{w,1},...,f^{w,L}\) the layers of \(f^{w}\) and we have that \(w=\left(\{W_{i}\}_{i=1}^{L}\right)\cup\left(\{b_{i}\}_{i=1}^{L}\right)\), where \(W_{i}\) and \(b_{i}\) represent weights and biases of the \(i-\)th layer of \(f^{w}.\) For \(x\in\mathbb{R}^{n}\) the output of layer \(i\in\{1,...,L\}\) can be
explicitly written as \(f^{w,i}(x)=a(W_{i}f^{w,i-1}(x)+b_{i})\) with \(f^{w,1}(x)=a(W_{1}x+b_{1})\), where \(a:\mathbb{R}\rightarrow\mathbb{R}\) is the activation function. We assume that \(a\) is a continuous monotonic function, which holds for the vast majority of activation functions used in practice such as sigmoid, ReLu, and tanh (Goodfellow et al., 2016). This guarantees that \(f^{w}\) is a continuous function.

Bayesian Neural Networks (BNNs), denoted by \(f^{\mathbf{w}}\), extend NNs by placing a prior distribution over the network parameters, \(p_{\mathbf{w}}(w)\), with \(\mathbf{w}\) being the vector of random variables associated to the parameter vector \(w\). Given a dataset \(\mathcal{D}\), training a BNN on \(\mathcal{D}\) requires to compute posterior distribution, \(p_{\mathbf{w}}(w|\mathcal{D})\), which can be computed via Bayes' rule (Neal, 2012). Unfortunately, because of the non-linearity introduced by the neural network architecture, the computation of the posterior is generally intractable. Hence, various approximation methods have been studied to perform inference with BNNs in practice. Among these methods, we consider Hamiltonian Monte Carlo (HMC) (Neal, 2012), and Variational Inference (VI) (Blundell et al., 2015). In our experimental evaluation in Section 6.3 we employ both HMC and VI.

Hamiltonian Monte Carlo (HMC).HMC proceeds by defining a Markov chain whose invariant distribution is \(p_{\mathbf{w}}(w|\mathcal{D})\), and relies on Hamiltonian dynamics to speed up the exploration of the space. Differently from VI discussed below, HMC does not make any parametric assumptions on the form of the posterior distribution and is asymptotically correct. The result of HMC is a set of samples that approximates \(p_{\mathbf{w}}(w|\mathcal{D})\). We refer interested readers to (Neal et al., 2011; Izmailov et al., 2021) for further details.

Variational Inference (VI).VI proceeds by finding a Gaussian approximating distribution \(q(w)\sim p_{\mathbf{w}}(w|\mathcal{D})\) in a trade-off between approximation accuracy and scalability. The core idea is that \(q(w)\) depends on some hyperparameters that are then iteratively optimized by minimizing a divergence measure between \(q(w)\) and \(p_{\mathbf{w}}(w|\mathcal{D})\). Samples can then be efficiently extracted from \(q(w)\). See (Khan and Rue, 2021; Blundell et al., 2015) for recent developments in variational inference in deep learning.
## 3 Problem Formulation

Given a trained BNN \(f^{\mathbf{w}}\) we consider the following discrete-time stochastic process given by iterative predictions of the BNN:

\[\mathbf{x}_{k}=f^{\mathbf{w}}(\mathbf{x}_{k-1},\mathbf{u}_{k-1})+\mathbf{v}_{k}, \quad\mathbf{u}_{k}=\pi_{k}(\mathbf{x}_{k}),\quad k\in\mathbb{N}_{>0}, \tag{1}\]

where \(\mathbf{x}_{k}\) is a random variable taking values in \(\mathbb{R}^{n}\) modelling the state of System (1) at time \(k\), \(\mathbf{v}_{k}\) is a random variable modelling an additive noise term with stationary, zero-mean Gaussian distribution \(\mathcal{N}(0,\sigma^{2}\cdot I)\), where \(I\) is the identity matrix of size \(n\times n\). \(\mathbf{u}_{k}\) represents the action applied at time \(k\), selected from a compact set \(\mathcal{U}\subset\mathbb{R}^{c}\) by a (deterministic) feedback Markov strategy (a.k.a. policy, or controller) \(\pi:\mathbb{R}^{n}\times\mathbb{N}\to\mathcal{U}\).1

Footnote 1: We can limit ourselves to consider deterministic Markov strategies as they are optimal in our setting (Bertsekas and Shreve, 2004; Abate et al., 2008). Also, in the following, we denote with \(\pi\) the time-varying policy described, at each step \(k\), by policy \(\pi_{k}:\mathbb{R}^{n}\to\mathcal{U}\).

The model in Eqn. (1) is commonly employed to represent noisy dynamical models driven by a BNN and controlled by the policy \(\pi\)(Depeweg et al., 2017). In this setting, \(f^{\mathbf{w}}\) defines the transition probabilities of the model and, correspondingly, \(p(\bar{x}|(x,u),\mathcal{D})\) is employed to describe the _posterior predictive distribution_, namely the probability density of the model state at the next time step being \(\bar{x}\), given that the current state and action are \((x,u)\), as:

\[p(\bar{x}|(x,u),\mathcal{D})=\int_{\mathbb{R}^{n_{w}}}\mathcal{N}(\bar{x}\mid f ^{w}(x,u),\sigma^{2}\cdot I)p_{\mathbf{w}}(w|\mathcal{D})dw, \tag{2}\]

where \(\mathcal{N}(\cdot\mid f^{w}(x,u),\sigma^{2}\cdot I)\) is the Gaussian likelihood induced by noise \(\mathbf{v}_{k}\) and centered at the NN output (Neal, 2012).

Observe that the posterior predictive distribution induces a probability density function over the state space. In iterative prediction settings, this implies that at each step the state vector \(\mathbf{x}_{k}\) fed into the BNN is a random variable. Hence, a \(N\)-step _trajectory_ of the dynamic model in Eqn (1) is a sequence of states \(x_{0},...,x_{N}\in\mathbb{R}^{n}\) sampled from the predictive distribution. As a consequence, a principled propagation of the BNN uncertainty through consecutive time steps poses the problem of predictions over stochastic inputs. In Section 4.1 we will tackle this problem for the particular case of reach-avoid properties, by designing a backward computation scheme that starts its calculations from the goal region, and proceeds according to Bellman iterations (Bertsekas and Shreve, 2004).
We remark that \(p(\bar{x}|(x,u),\mathcal{D})\) is defined by marginalizing over \(p_{\mathbf{w}}(w|\mathcal{D})\), hence, the particular \(p(\bar{x}|(x,u),\mathcal{D})\) depends on the specific approximate inference method employed to estimate the posterior distribution. As such, the results that we derive are valid w.r.t. a specific BNN posterior.

Probability MeasureFor an action \(u\in\mathbb{R}^{c}\), a subset of states \(X\subseteq\mathbb{R}^{n}\) and a starting state \(x\in\mathbb{R}^{n}\), we call \(T(X|x,u)\) the _stochastic kernel_ associated (and equivalent (Abate et al., 2008)) to the dynamical model of Equation (1). Namely, \(T(X|x,u)\) describes the one-step transition probability of the model of Eqn. (1) and is defined by integrating the predictive posterior distribution with input \((x,u)\) over \(X\), as:

\[T(X|x,u)= \int_{X}p(\bar{x}|(x,u),\mathcal{D})d\bar{x}. \tag{3}\]

In what follows, it will be convenient at times to work over the space of parameters of the BNN. To do so, we can re-write the stochastic kernel by combining Equations (2) and (3) and applying Fubini's theorem (Fubini, 1907) to switch the integration order, thus obtaining:

\[T(X|x,u)=\int_{\mathbb{R}^{n_{w}}}\left[\int_{X}\mathcal{N}(\bar{x}|f^{w}(x,u ),\sigma^{2}\cdot I)d\bar{x}\right]p_{\mathbf{w}}(w|\mathcal{D})dw. \tag{4}\]

From this definition of \(T\) it follows that, under a strategy \(\pi\) and for a given initial condition \(x_{0}\), \(\mathbf{x}_{k}\) is a Markov process with a well-defined probability measure \(\Pr\) uniquely generated by the stochastic kernel \(T\)(Bertsekas and Shreve, 2004, Proposition 7.45) and such that for \(X_{0},X_{k}\subseteq\mathbb{R}^{n}\):

\[\Pr[\mathbf{x}_{0}\in X_{0}]=\mathbf{1}_{X_{0}}(x_{0}),\] \[\Pr[\mathbf{x}_{k}\in X_{k}|\mathbf{x}_{k-1}=x,\pi]=T(X_{k}|x, \pi_{k-1}(x)),\]

where \(\mathbf{1}_{X_{0}}\) is the indicator function (that is, 1 if \(x\subseteq X_{0}\) and 0 otherwise). Having a definition of \(\Pr\) allows one to make probabilistic statements over the stochastic model in Eqn (1).

**Remark 1**.: _Note that, as is common in the literature (Depeweg et al., 2017), according to the definition of the probability measure \(\Pr\) we marginalise over the posterior distribution at each time step. Consequently, according to our modelling framework, the weights of the BNN are not kept fixed during each trajectory, but we re-sample from \(\mathbf{w}\) at each time step._
### Problem Statements

We consider two problems concerning, respectively, the certification and the control of dynamical systems modelled by BNNs. We first consider safety certification with respect to probabilistic reach-avoid specifications. That is, we seek to compute the probability that from a given state, under a selected control policy, an agent navigates to the goal region without encountering any unsafe states. Next, we consider the formal synthesis of policies that maximise this probability and thus attain maximal certifiable safety.

**Problem 1** (Computation of Probabilistic Reach-Avoid): _Given a strategy \(\pi\), a goal region \(\mathrm{G}\subseteq\mathbb{R}^{n}\), a finite-time horizon \([0,N]\subseteq\mathbb{N},\) and a safe set \(\mathrm{S}\subseteq\mathbb{R}^{n}\) such that \(\mathrm{G}\cap\mathrm{S}=\emptyset\), compute for any \(x_{0}\in\mathrm{G}\cup\mathrm{S}\)_

\[P_{reach}(\mathrm{G},\mathrm{S},x_{0},[0,N]|\pi)=\] \[\Pr\big{[}\exists k\in[0,N],\mathbf{x}_{k}\in\mathrm{G}\,\wedge \forall 0\leq k^{\prime}<k,\mathbf{x}_{k^{\prime}}\in\mathrm{S}\mid\mathbf{x}_{0}= x_{0},\pi\big{]}. \tag{5}\]

Outline of the ApproachIn Section 4.1 we show how \(P_{reach}(\mathrm{G},\mathrm{S},x_{0},[0,N]|\pi)\) can be formulated as the solution of a backward iterative computational procedure, where the uncertainty of the BNN is propagated backward in time, starting from the goal region. Our approach allows us to compute a sound lower bound on \(P_{reach}\), thus guaranteeing that \(\mathbf{x}_{k}\), as defined in Eqn (1), satisfies the specification with a given probability. This is achieved by extending existing lower bounding techniques developed to certify BNNs (Wicker et al., 2020) and applying these at each propagation step through the BNN.

Note that, in Problem 1, the strategy \(\pi\) is provided, and the goal is to quantify the probability with which the trajectories of \(\mathbf{x}_{k}\) satisfy the given specification. In Problem 2 below, we expand the previous problem and seek to synthesise a controller \(\pi\) that maximizes \(P_{reach}\). The general formulation of this optimization is given below.

**Problem 2** (Strategy Synthesis for Probabilistic Reach-Avoid): _For an initial state \(x_{0}\in\mathrm{G}\cup\mathrm{S}\), and a finite time horizon \(N\), find a strategy \(\pi^{*}:\mathbb{R}^{n}\times\mathbb{R}_{\geq 0}\to\mathbb{R}^{c}\) such that_

\[\pi^{*}=\operatorname*{arg\,max}_{\pi}P_{reach}(\mathrm{G},\mathrm{S},x_{0},[0,N]\mid\pi). \tag{6}\]

In Section 5, we will provide specific schemes for synthesizing optimal strategies when \(\pi\) is either a look-up table or a deterministic neural network.
Outline of the ApproachTo solve this problem, we notice that the backward iterative procedure outlined to solve Problem 1 has a substructure such that dynamic programming will allow us to compute optimal actions for each state that we verify, thus producing an optimal policy with respect to the given posterior and reach-avoid specification. With low-dimensional or discrete action spaces, we can then derive a tabular policy by solving the resulting dynamic programming problem. For higher-dimensional action spaces instead, in Section 5.1 we consider (generalising) policies represented as neural networks.

## 4 Methodology

In this section, we illustrate the methodology used to compute lower bounds on the reach-avoid probability, as described in Problem 1. We begin by encoding the reach-avoid probability through a sequence of value functions.

### Certifying Reach-Avoid Specifications

We begin by showing that \(P_{reach}(\mathrm{G},\mathrm{S},x,[k,N]|\pi)\) can be obtained as the solution of a backward iterative procedure, which allows to compute a lower bound on its value. In particular, given a time \(0\leq k<N\) and a strategy \(\pi\), consider the value functions \(V_{k}^{\pi}:\mathbb{R}^{n}\to[0,1]\), recursively defined as

\[V_{N}^{\pi}(x) =\mathbf{1}_{\mathrm{G}}(x),\] \[V_{k}^{\pi}(x) =\mathbf{1}_{\mathrm{G}}(x)+\mathbf{1}_{\mathrm{S}}(x)\!\!\! \int\!\!V_{k+1}^{\pi}(\bar{x})p\big{(}\bar{x}|(x,\pi_{k}(x)),\mathcal{D}\big{)} d\bar{x}. \tag{7}\]

Intuitively, \(V_{k}^{\pi}\) is computed starting from the goal region \(\mathrm{G}\) at \(k=N\), where it is initialised at value \(1\). The computation proceeds backwards at each state \(x\), by combining the current values with the transition probabilities from Eqn (1). The following proposition, proved inductively over time in the Supplementary Material, guarantees that \(V_{0}^{\pi}(x)\) is indeed equal to \(P_{reach}(\mathrm{G},\mathrm{S},x,[0,N]|\pi)\).

**Proposition 1**.: _For \(0\leq k\leq N\) and \(x_{0}\in\mathrm{G}\cup\mathrm{S},\) it holds that_

\[P_{reach}(\mathrm{G},\mathrm{S},x_{0},[k,N]|\pi)=V_{k}^{\pi}(x).\]
The backward recursion in Eqn (7) does not generally admit a solution in closed-form, as it would require integrating over the BNN posterior predictive distribution, which is in general analytically intractable. In the following section, we present a computational scheme utilizing convex relaxations to lower bound \(P_{reach}\).

### Lower Bound on \(P_{reach}\)

We develop a computational approach based on the discretisation of the state space, which allows convex relaxation methods such as (Wicker et al., 2020) to be used. The proposed computational approach is illustrated in Figure 1 and formalized in Section 4.3. Let \(Q=\{q_{1},...,q_{n_{q}}\}\) be a partition of \(\mathrm{S}\cup\mathrm{G}\) in \(n_{q}\) regions and denote with \(z:\mathbb{R}^{n}\to Q\) the function that associates to a state in \(\mathbb{R}^{n}\) the corresponding partitioned state in \(Q\). For each \(0\leq k\leq N\) we iteratively build a set of functions \(K_{k}^{\pi}:Q\rightarrow[0,1]\) such that for all \(x\in\mathrm{G}\cup\mathrm{S}\) we have that \(K_{k}^{\pi}(z(x))\leq V_{k}^{\pi}(x)\). Intuitively, \(K_{k}^{\pi}\) provides a lower bound for the value functions on the computation of \(P_{reach}\).

Figure 1: Examples of functions \(K_{N-1}^{\pi}\) (left) and \(K_{N-2}^{\pi}\) (right), which are lower bounds of \(V_{k}^{\pi}\). On the left, we consider the first step of our backward algorithm, where we compute \(K_{N-2}^{\pi}(q)\) by computing the probability that \(\mathbf{x}_{N}\in\mathrm{G}\) given that \(\mathbf{x}_{N-1}\in q.\) On the right, we consider the subsequent step. We outline the state we want to verify in red and the goal region in green. With the orange arrow, we represent the \(0.95\) transition probability of the BNN dynamical model, and in pink we represent the worst-case probabilities spanned by the BNN output. On top, we show where each of these key terms comes into play in Eqn (9).
The functions \(K_{k}^{\pi}\) are obtained by propagating backward the BNN predictions from time \(N\), where we set \(K_{N}^{\pi}(q)=\mathbf{1}_{\mathrm{G}}(q)\), with \(\mathbf{1}_{\mathrm{G}}(q)\) being the indicator function (that is, \(1\) if \(q\subseteq\mathrm{G}\) and \(0\) otherwise). Then, for each \(k<N\), we first discretize the set of possible probabilities in \(n_{p}\) sub-intervals \(0=v_{0}\leq v_{1}\leq...\leq v_{n_{p}}=1\). Hence, for any \(q\in Q\) and probability interval \([v_{i},v_{i+1}]\), one can compute a lower bound, \(R(q,k,\pi,i)\), on the probability that, starting from any state in \(q\) at time \(k\), we reach in the next step a region that has probability \(\in[v_{i},v_{i+1}]\) of safely reaching the goal region. The resulting values are used to build \(K_{k}^{\pi}\) (as we will detail in Eqn (9)). For a given \(q\subset\mathrm{S}\), \(K_{k}^{\pi}(q)\) is obtained as the sum over \(i\) of \(R(q,k,\pi,i)\) multiplied by \(v_{i-1}\), i.e., the lower value that \(K_{k+1}^{\pi}\) obtains in all the states of the \(i-th\) region. Note that the discretisation of the probability values does not have to be uniform, but can be adaptive for each \(q\in Q\). A heuristic for picking the value of thresholds \(v_{i}\) will be given in Algorithm 1. In what follows, we formalise the intuition behind this computational procedure.

### Lower Bounding of the Value Functions

For a given strategy \(\pi\), we consider a constant \(\eta\in(0,1)\) and \(\epsilon=\sqrt{2\sigma^{2}}\mathrm{erf}^{-1}(\eta)\), which are used to bound the value of the noise, \(\mathbf{v}_{k}\), at any given time. Intuitively, \(\eta\) represents the proportion of observational error we consider.2 Then, for \(0\leq k<N\), \(K_{k}^{\pi}:Q\to[0,1]\) are defined recursively as follows:

Footnote 2: The threshold is such that it holds that \(Pr(|\mathbf{v}_{k}^{(i)}|\leq\epsilon)=\eta\). In the experiments of Section 6 we select \(\eta=0.99\).

\[K_{N}^{\pi}(q) =\mathbf{1}_{\mathrm{G}}(q), \tag{8}\] \[K_{k}^{\pi}(q) =\mathbf{1}_{\mathrm{G}}(q)+\mathbf{1}_{\mathrm{S}}(q)\sum_{i=1 }^{n_{p}}v_{i-1}R(q,k,\pi(q),i), \tag{9}\]

where

\[R(q,k,\pi(q),i)=\eta^{n}\int_{H_{k,i}^{q,\pi,\epsilon}}p_{ \mathbf{w}}(w|\mathcal{D})dw, \tag{10}\] \[H_{k,i}^{q,\pi,\epsilon} =\{w\in\mathbb{R}^{n_{w}}|\,\forall x\in q,\forall\gamma\in[- \epsilon,\epsilon]^{n},\text{ it holds that:}\] \[v_{i-1} \leq K_{k+1}^{\pi}(q^{\prime})\leq v_{i},\text{ with }q^{ \prime}=z(f^{w}(x,\pi_{k}(x))+\gamma)\}.\]
The key component for the above backward recursion is \(R(q,k,\pi,i)\), which bounds the probability that, starting from \(q\) at time \(k\), we have that \(\mathbf{x}_{k+1}\) will be in a region \(q^{\prime}\) such that \(K_{k+1}^{\pi}(q^{\prime})\in[v_{i},v_{i+1}]\). By definition, the set \(H_{k,i}^{q,\pi,\epsilon}\) defines the weights for which the BNN maps all states covered by \(q\) into the goal states given action \(\pi(q)\). Given this, it is clear that integration of the posterior \(p_{\mathbf{w}}(w|\mathcal{D})\) over the \(H_{k,i}^{q,\pi,\epsilon}\) will return the probability mass of system (1) transitioning from \(q\) to \(q^{\prime}\) with probability in \([v_{i},v_{i+1}]\) in one time step. The computation of Eqn (9) then reduces to computing the set of weights \(H_{k,i}^{q,\pi,\epsilon}\), which we call the _projecting weight set_. A method to compute a safe under-approximation \(\bar{H}\subseteq H_{k,i}^{q,\pi,\epsilon}\) is discussed below. Before describing that, we analyze the correctness of the above recursion.

**Theorem 1**.: _Given \(x\in\mathbb{R}^{n}\), for any \(k\in\{0,...,N\}\) and \(q=z(x)\), assume that \(H_{k,i}^{q,\pi,\epsilon}\cap H_{k,j}^{q,\pi,\epsilon}=\emptyset\) for \(i\neq j.\) Then:_

\[\inf_{x\in q}V_{k}^{\pi}(x)\geq K_{k}^{\pi}(q).\]

A proof of Theorem 1 is given in the Supplementary Material. Note that the assumption on the null intersection between different projecting weight sets required in Theorem 1 can always be enforced by taking their intersection and complement.

### Computation of Projecting Weight Set

Theorem 1 allows us to compute a safe lower bound to Problem 1, by relying on an abstraction of the state space, that is, through the computation of \(K_{0}^{\pi}(q)\). This can be evaluated once the projecting set of weight values \(H_{k,i}^{q,\pi,\epsilon}\) associated to \([v_{i-1},v_{i}]\) is known.3 Unfortunately, direct computation of \(H_{k,i}^{q,\pi,\epsilon}\) is intractable. Nevertheless, a method for its lower bounding was developed by Wicker et al. (2020) in the context of adversarial perturbations for one-step BNN predictions, and can be directly adapted to our settings.

Footnote 3: In the case of Gaussian VI the integral of Equation (10) can be computed in terms of the _erf_ function, whereas more generally Monte Carlo or numerical integration techniques can be used.

The idea is that an under approximation \(\bar{H}\subseteq H_{k,i}^{q,\pi,\epsilon}\) is built by sampling weight boxes of the shape \(\hat{H}=[w^{L},w^{U}]\), according to the posterior, and
checking whether:

\[v_{i-1}\leq K_{k+1}^{\pi}(z(f^{w}(x,\pi_{k}(x))+\gamma))\leq v_{i},\] \[\forall x\in q,\,\forall w\in\hat{H},\,\forall\gamma\in\left[- \epsilon,\epsilon\right]^{n}. \tag{11}\]

Finally, \(\bar{H}\) is built as a disjoint union of boxes \(\hat{H}\) satisfying the above condition. For a full discussion of the details of this method we refer interested readers to (Wicker et al., 2020). In order to apply this method to our setting, we propagate the abstract state \(q\) through the policy function \(\pi_{k}(x)\), so as to obtain a bounding box \(\widehat{\Pi}=[\pi^{L},\pi^{U}]\) such that \(\pi^{L}\leq\pi_{k}(x)\leq\pi^{U}\) for all \(x\in q\). In the experiments, this bounding is only necessary when \(\pi_{k}(x)\) is given by an NN controller, for which bound propagation of NNs can be used for the computation of \(\widehat{\Pi}\)(Gowal et al., 2018; Gehr et al., 2018).

The results of Proposition 2 and Proposition 3 from Wicker et al. (2020) can then be used to propagate \(q\), \(\widehat{\Pi}\) and \(\hat{H}\) through the BNN. For discrete posteriors (e.g., those resulting from HMC) one can use the method described by Gowal et al. (2018) (Equations 6 and 7). Propagation of \(q\), \(\widehat{\Pi}\) amounts to using these method to compute values \(f^{L}_{q,\epsilon,k}\) and \(f^{U}_{q,\epsilon,k}\) such that, for all \(x\in q,\gamma\in[-\epsilon,\epsilon]^{n},w\in\hat{H}\), it holds that:

\[f^{L}_{q,\epsilon,k}\leq f^{w}(x,\pi_{k}(x))+\gamma\leq f^{U}_{q,\epsilon,k}. \tag{12}\]

Furthermore, \(f^{L}_{q,\epsilon,k}\) and \(f^{U}_{q,\epsilon,k}\) are differentiable w.r.t. to the input vector (Gowal et al., 2018; Wicker et al., 2021a).

Finally, the two bounding values can be used to check whether or not the condition in Eqn (11) is satisfied, by simply checking whether \([f^{L}_{q,\epsilon,k},f^{U}_{q,\epsilon,k}]\) propagated through \(K_{k+1}^{\pi}\) is within \([v_{i},v_{i+1}]\). We highlight that computing this probability is equivalent to a conservative estimate of \(R(q,k,\pi,i)\).

### Probabilistic Reach-Avoid Algorithm

In Algorithm 1 we summarize our approach for computing a lower bound for Problem 1. For simplicity of presentation, we consider the case \(n_{p}=2\), (i.e., we partition the range of probabilities in just two intervals \([0,v_{1}]\), \([v_{1},1]\) - the case \(n_{p}>2\) follows similarly). The algorithm proceeds by first initializing the reach-avoid probability for the partitioned states \(q\) inside the goal region \(G\) to \(1\), as per Eqn (8). Then, for each of the \(N\) time steps and for each one of the remaining abstract states \(q\), in line 4 we set the threshold probability \(v_{1}\) equal to the maximum value that \(K^{\pi}\) attains at the next time step over the
states in the neighbourhood of \(q\) (which we capture with a hyper-parameter \(\rho_{x}>0\)). We found this heuristic for the choice of \(v_{1}\) to work well in practice (notice that the obtained bound is formal irrespective of the choice of \(v_{1}\), and different choices could potentially be explored). We then proceed in the computation of Eqn (9). This computation is performed in lines 5-14. First, we initialise to the null set the current under-approximation of the projecting weight set, \(\bar{H}\). We then sample \(n_{s}\) weights boxes \(\hat{H}\) by sampling weights from the posterior, and expanding them with a margin \(\rho_{w}\) heuristically selected (lines 6-8). Then, for each of these sets, we first propagate the state \(q\), policy function, and weight set \(\bar{H}\) to build a box \(\bar{X}\) according to Eqn (12) (line 9), which is then accepted or rejected based on the value that \(K^{\pi}\) at the next time step attains in states in \(\bar{X}\) (lines 10-12). \(K^{\pi}_{N-i}(q)\) is then computed in line 14 by integrating \(p_{\mathbf{w}}(w|\mathcal{D})\) over the union of the accepted sets of weights.

**Input:** BNN model \(f^{\mathbf{w}}\), safe region S, goal region G, discretization \(Q\) of \(\mathrm{S}\cup\mathrm{G}\), time horizon \(N\), neural controller \(\pi\), number of BNN samples \(n_{s}\), weight margin \(\rho_{w}\), state space margin \(\rho_{x}\)

**Output:** Lower bound on \(V^{\pi}\)

```
1:For all \(0\leq k\leq N\) set \(K^{\pi}_{k}(q)=1\) iff \(q\subseteq\mathrm{G}\) and \(0\) otherwise
2:for\(k\gets N\) to \(1\)do
3:for\(q\in Q\setminus\mathrm{G}\)do
4:\(v_{1}\leftarrow\max_{x\in[q-\rho_{x},q+\rho_{x}]}K^{\pi}_{k+1}(z(x))\)
5:\(\bar{H}\leftarrow\emptyset\quad\)# \(\bar{H}\) is the set of safe weights
6:for desired number of samples, \(n_{s}\)do
7:\(w^{\prime}\sim P(w|\mathcal{D})\)
8:\(\hat{H}\leftarrow[w^{\prime}-\rho_{w},w^{\prime}+\rho_{w}]\)
9:\(\bar{X}\leftarrow[f^{L}_{q,\epsilon,k},f^{U}_{q,\epsilon,k}]\quad\)# Computed according to Eqn (12)
10:if\(\min_{x\in\bar{X}}K^{\pi}_{k+1}(z(x))\geq v_{1}\)then
11:\(\bar{H}\leftarrow\bar{H}\bigcup\hat{H}\)
12:endif
13:endfor
14: Ensure \(H_{i}\cap H_{j}=\emptyset\quad\forall H_{i},H_{j}\in\bar{H}\)
15:\(K^{\pi}_{k}(q)=v_{1}\cdot\eta^{n}\int_{\bar{H}}p_{\mathbf{w}}(w|\mathcal{D})dw\quad\)(Eqn (9))
16:endfor
17:endfor
18:return\(K^{\pi}\)
```

**Algorithm 1** Probabilistic Reach-Avoid for BNNs
## 5 Strategy Synthesis

We now focus on synthesising a strategy that maximizes our lower bound on \(P_{reach}\), thus solving Problem 2. Notice that, while no global optimality claim can be made about the strategy that we obtain, maximising the lower bound guarantees that the true reach-avoid probability will still be greater than the improved bound obtained after the maximisation.

**Definition 1**.: _A strategy \(\pi^{*}\) is called maximally certified (max-cert), w.r.t. the discretised value function \(K^{\pi}\), if and only if, for all \(x\in\mathrm{G}\cup\mathrm{S}\), it satisfies_

\[K_{0}^{\pi^{*}}(z(x))=\sup_{\pi}K_{0}^{\pi}(z(x)),\]

_that is, the strategy \(\pi^{*}\) maximises the lower bound of \(P_{reach}\)._

It follows that, if \(K_{0}^{\pi^{*}}(z(x))>1-\delta\) for all \(x\in\mathrm{G}\cup\mathrm{S}\), then the max-cert strategy \(\pi^{*}\) is a solution of Problem 2. Note that a max-cert strategy is guaranteed to exist when the set of admissible controls \(\mathcal{U}\) is compact (Bertsekas and Shreve, 2004, Lemma 3.1), as we assume in this work. In the next theorem, we show that a max-cert strategy can be computed via dynamic programming with a backward recursion similar to that of Eqn (9).

**Theorem 2**.: _For \(0\leq k<N\) and \(0=v_{0}<...<v_{n_{p}}=1,\) define the functions \(K_{k}^{*}:\mathbb{R}^{n}\to[0,1]\) recursively as follows_

\[K_{k}^{*}(q)=\sup_{u\in\mathcal{U}}\big{(}\mathbf{1}_{\mathrm{G}}(q)+\mathbf{ 1}_{\mathrm{S}}(q)\sum_{i=1}^{n_{p}}v_{i}R(q,k,u,i)\big{)},\]

_where \(R(q,k,u,i)\) and \(H_{k,i}^{q,u,\epsilon}\) are defined as in Eqn (10). If \(\pi^{*}\) is s.t. \(K_{0}^{*}=K_{0}^{\pi^{*}}\), then \(\pi^{*}\) is a max-cert strategy. Furthermore, for any \(x\), it holds that \(K_{0}^{\pi^{*}}(z(x))\leq P_{reach}(\mathrm{G},\mathrm{S},[0,N],x|\pi^{*})\)._

Theorem 2 is a direct consequence of the Bellman principle of optimality (Abate et al., 2008, Theorem 2) and it guarantees that for each state \(q\in\mathrm{S}\) and time \(k\), we have that \(\pi^{*}(q,k)=\arg\max_{u\in\mathcal{U}}\sum_{i=1}^{n_{p}}v_{i}R(q,k,u,i)\).

In Algorithm 2 we present a numerical scheme based on Theorem 2 to find a max-cert policy \(\pi^{*}\). Note that the optimization problem required to be solved at each time step state, i.e. \(\arg\max_{u\in\mathcal{U}}\sum_{i=1}^{n_{p}}v_{i}R(q,k,u,i)\), is non-convex. Hence, in Algorithm 2, in Line 1, we start by partitioning the action space \(\mathcal{U}\). Then, in Lines 4-15 for each action in the partition we estimate
the expectation of \(K_{k+1}^{*}\) starting from \(q\) via \(n_{s}\) samples taken from the BNN posterior (250 in all our experiments). Finally, in Lines 11-14 we keep track of the action maximising \(K_{k+1}^{*}\).

The described approach for synthesis, while optimal in the limit of an infinitesimal discretization of \(\mathcal{U}\), may become infeasible for large state and action spaces. As a consequence, in the next subsection, we also consider when \(\pi\) is parametrised by a neural network and thus can serve as a function over a larger (even infinite) state space. Specifically, we show how a set of neural controllers, one for each time step, can be trained in order to maximize probabilistic reach-avoid via Theorem 2. In Section 6 we empirically investigate both controller strategies.

```
1:\(\Upsilon\leftarrow\) middle points of each region in a partition of \(\mathcal{U}\)
2:\(\kappa^{*}\gets 0\)
3:\(u^{*}\gets 0\)
4:for\(u\in\Upsilon\)do
5:\(\hat{\kappa}\gets 0\)
6:for\(j\) from \(0\) to \(n_{s}\)do
7:\(w^{\prime}\sim P(w|\mathcal{D})\)
8:\(\bar{X}\leftarrow[f_{q,\epsilon,k}^{L},f_{q,\epsilon,k}^{U}]\)# Computed for \(q\) and \(u\) via (Eqn (12))
9:\(\hat{\kappa}=\hat{\kappa}+\frac{\min_{x\in\bar{X}}K_{k+1}^{*}(z(x))}{n_{s}}\)
10:endfor
11:if\(\hat{\kappa}>\kappa^{*}\)then
12:\(\kappa^{*}\leftarrow\hat{\kappa}\)
13:\(u^{*}\gets u\)
14:endif
15:endfor
16:return\(u^{*}\)
```

**Algorithm 2** Numerical Synthesis of Action for region \(q\) at time \(k\)

### An Approach for Strategy Synthesis Based on Neural Networks

In this subsection we show how we can train a set of NN policies \(\pi_{0},...,\pi_{N-1}:\mathbb{R}^{n}\rightarrow\mathcal{U}\) such that at each time step \(k\), \(\pi_{k}\) approximately solves the dynamic programming equation in Theorem 2. Note that, because of the approximate nature of the NN training, the resulting neural policies will necessarily
be sub-optimal, but have the potential to scale to larger and more complex systems, compared to the approach presented in Algorithm 2.

At time \(k\) we start with an initial set of parameters (weights and biases) \(\theta_{k}\) for policy \(\pi_{k}\). These parameters can either be initialized to \(\theta_{k+1}\), the parameters synthesised at the previous time step of the value iteration for \(\pi_{k+1}\), or to a policy employed to collect the data to train the BNN as in Gal et al. (2016b), or simply selected at random. In our implementation where no previous policy is available, we start with a randomly initialized NN, and then at time \(k\) we set our initial neural policy with that obtained at time \(k+1\). We then employ a scheme to learn a "safer" set of parameters via backpropagation. In particular, we first define the following loss function penalizing policy parameters that lead to an unsafe behavior for an ensemble of NNs sampled from the BNN posterior distribution:

\[\mathcal{L}(x,\theta_{k})=-\alpha||\sum_{w\in\bar{W}}f^{w}(x,\pi_ {k}(x))-\mathbf{A}_{k}||_{2}+(1-\alpha)||\sum_{w\in\bar{W}}f^{w}(x,\pi_{k}(x)) -\mathbf{R}_{k}||_{2}, \tag{13}\]

where \(\bar{W}\) are a set of parameters independently sampled from the BNN posterior \(p_{\mathbf{w}}(w|\mathcal{D})\), for a probability threshold \(0\leq p_{t}\leq 1\), \(\mathbf{A}_{k}=\{x:K_{k+1}^{\pi_{k+1}}(x)\geq p_{t}\}\) and \(\mathbf{R}_{k}=\{x:K_{k+1}^{\pi_{k+1}}(x)\leq 1-p_{t}\}\) are the sets of states for which the probability of satisfying the specification at time \(k+1\) is respectively greater than \(p_{t}\) and smaller than \(1-p_{t}\). For \(X\subset\mathbb{R}^{n}\), \(||x-X||_{2}=\inf_{\bar{x}\in X}||x-\bar{x}||_{2}\) is the standard \(L_{2}\) distance of a point from a set, and \(0\leq\alpha\leq 1\) is a parameter taken to be \(0.25\) in our experiments, that weights between reaching the goal and staying away from "bad" states. Intuitively, the first term in \(\mathcal{L}(x,\theta_{k})\) enforces \(\theta_{k}\) that leads to high values of \(K_{k+1}^{\pi_{k+1}}\), while the second term penalizes parameter sets that lead to small values of this quantity.

\(\mathcal{L}(x,\theta_{k})\) only considers the behaviour of the dynamical system of Equation (1) starting from initial state \(x\). Then, in order to also enforce robustness in a neighbourhood of initial states around \(x\), similarly to the adversarial training case (Madry et al., 2017), we consider the robust loss

\[\bar{\mathcal{L}}(x,\theta_{k})=\max_{x^{\prime}:||x-x^{\prime}|| _{2}\leq\epsilon}\mathcal{L}(x,\theta_{k}), \tag{14}\]

Note that by employing Eqn (12) we obtain a differentiable upper bound of \(\bar{\mathcal{L}}(x,\theta_{k})\), which can be employed for training \(\theta_{k}\).
### Discussion on the Algorithms

In this section we provide further discussion of our proposed algorithms including the complexity and the various sources of approximation that may lead to looser guarantees. To frame this discussion, we start by highlighting the complexity and approximation introduced by the chosen bound propagation technique shared by both of the algorithms. We then proceed to discuss how discretisation choices made with respect to the state-space, the weight-space, and the observational noise, practically affect the tightness of our probability bounds for both algorithms, and finally how the action-space discretisation affects our synthesis algorithm.

Bound Propagation TechniquesGiven that there are currently no methods for BNN certification that are both sound and complete (Wicker et al., 2020, 2021b; Berrada et al., 2021), the evaluation of the \(R\) function will always introduce some approximation error. While it is difficult to characterize this error in general, it is known that for deeper networks, BNN certification methods introduce more approximation than shallow networks (Wicker, 2021). The recently developed bounds from Berrada et al. (2021), have been shown to be tighter than the IBP and LBP approaches from Wicker et al. (2020) at the cost of computation complexity that is exponential in the number of dimensions of the state-space. In contrast, each iteration of the interval bound propagation method proposed in and Wicker et al. (2020) requires the computational complexity of four forward passes through the neural network architecture.

Discretization Error and ComplexityWhile our formulation supports any form of state space discretisation, we can assume for simplicity that each dimension of the \(n\)-dimensional state-space is broken into \(m\) equal-sized abstract states. This implies that certification of the system requires us to evaluate the \(R\) function \(\mathcal{O}\big{(}N(m^{n})\big{)}\) many times, where \(N\) is the time horizon we would like to verify. Given that \(n\) is fixed, the user has control over \(m\), the size of each abstract state. For large abstract states, small \(m\), one introduces more approximation as the \(R\) function must account for all possible behaviors in the abstract state. For small abstract states, large \(m\), there is much less approximation, but considerably larger runtime. Assume the \(c\)-dimensional action space is broken into \(t\) equal portions at each dimension, then the computational complexity of the algorithm becomes \(\mathcal{O}\big{(}t^{c}N(m^{n})\big{)}\) as each of the \(m^{n}\) states must be evaluated \(t^{c}\)-many times to determine the approximately
optimal action. As with the state-space discretization, larger \(t\) will lead to a more-optimal action choice but requires greater computational time.

## 6 Experiments

We provide an empirical analysis of our BNN certification and policy synthesis methods. We begin by providing details on the experimental setting in Section 6.1. We then analyse the performance of our certification approach on synthesized policies in Section 6.2. Next, in Section 6.3, we discuss how the choice of the BNN inference algorithm affects synthesis and certification results. Finally, in Section 6.4 we study how our method scales with larger neural network architectures and in higher-dimensional control settings.

### Experimental Setting

We consider a planar control task consisting of a point-mass agent navigating through various obstacle layouts. The point-mass agent is described by four dimensions, two encoding position information and two encoding velocity (Astrom and Murray, 2008). To control the agent there are two continuous action dimensions, which represent the force applied on the point-mass in each of the two planar directions. The task of the agent is to navigate to a goal region while avoiding various obstacle layouts. The knowledge supplied to the agent about the environment is the locations of the goal and obstacles. The full set of equations describing the agent dynamics is given in Appendix A.1. In our experiments, we analyse three obstacle layouts of varying difficulty, which we name v1, v2 and Zigzag - visualized in the left column of Figure 3. Obstacle layout v1 places an obstacle directly between the agent's initial position and the goal, forcing the agent to navigate its way around it. Obstacle layout v2 extends this setting by adding two further obstacles that block off one side of the state space. Finally, scenario Zigzag has 5 interleaving triangles and requires the agent to navigate around them in order to reach the goal.

In order to learn an initial policy to solve the task, we employ the episodic learning framework described in Gal et al. (2016a). This consists of iteratively collecting data from deploying our learned policy, updating the BNN dynamics model to the new observations, and updating our policy. When collecting data, we start by randomly sampling state-action pairs and observing their resulting state according to the ground-truth dynamics. After this initial sample, all future observations from the ground-truth environment
are obtained from deploying our learned policy. The initial policy is set by assigning a random action to each abstract state. This is equivalent to tabular policy representations in standard reinforcement learning (Sutton and Barto, 1998). We additionally discuss neural network policies in Section 6.4. Actions in the policy are updated by performing gradient descent on a sum of discounted rewards over a pre-specified finite horizon. The reward of an action is taken to be the \(\ell_{2}\) distance moved towards the goal region penalized by the \(\ell_{2}\) proximity to obstacles as is done in (Sutton and Barto, 1998; Gal et al., 2016a). For the learning of the BNN, we perform approximate Bayesian inference over the neural network parameters. For our primary investigation, we select an NN architecture with a single fully connected hidden layer comprising 50 hidden units, and learn the parameters via Hamiltonian Monte Carlo (HMC). Larger neural network architectures are considered in Section 6.4, while results for variational approximate inference are given in Section 6.3.

Unless otherwise specified, in performing certification and synthesis we employ abstract states spanning a width of 0.02 around each position dimension and 0.08 around each velocity dimension. Velocity values are clipped to the range \([-0.5,0.1]\). When performing optimal synthesis, we discretise the two action dimensions for the point-mass problem into 100 possible vectors which uniformly cover the continuous space of actions \([-1,1]\). When running our backward reachability scheme, at each state, we test all 100 action vectors and take the action that maximizes our lower bound to be the policy action at that state, thus giving us the locally optimal action within the given discretisation. Further experimental details are presented in Appendix A and code to reproduce all results in this paper can be found at [https://github.com/matthewwicker/BNNReachAvoid](https://github.com/matthewwicker/BNNReachAvoid).

The computational times for each system were roughly equivalent. This is to be expected given that each has the same state space. The following average times are reported for a parallel implementation of our algorithm run on 90 logical CPU cores across 4 Intel Core Xeon 6230 clocked at 2.10GHz. Training of the initial policy and BNN model takes in the order of 10 minutes, 6 hours for the certification with a horizon of 50 time steps, and 8 hours for synthesis.
### Comparing Certification of Learned and Max-Cert Policies

In Figure 2 and Figure 3 we visualize systems from both learned and synthesized policies. Each row represents one of our control environments and is comprised of four figures. These figures show, respectively, simulations from the dynamical system, BNN uncertainty, the control policy plotted as a gradient field, and the certified safety probabilities. The first column of the Figures depicts 200 simulated trajectories of the learned (Figure 2) or synthesized (Figure 3) control policies on the BNN (whose uncertainty is plotted along the second column). Notice how in both cases we visually obtain the

Figure 2: **Left Column: 200 simulated trajectories for the learned policy starting from the initial state. Center Left Column: A 2D visualization of the learned policy. Each arrow represents the direction of the applied force. Center Right Column: The epistemic uncertainty for the learned dynamics model. Right Column: Certified lower-bounds of probabilistic reach avoid for each abstract state according to BNN and final learned policy.**
Figure 3: A version of each learned system after a new policy has been synthesized from the reach-avoid specification. **First column:** 200 simulations of the synthesized policy in the real environment. **Second column:** BNN epistemic uncertainty given as the variance of the BNN predictive distribution. **Third column:** A visualization of the maximally certifiable policies, which demonstrate a clearer tendency to avoid obstacles throughout the state space compared to the policies in Figure 2. **Fourth column:** Synthesized policies have remarkably higher lower-bounds than learned policies, corresponding plots for learned systems in Figure 2.
behaviour expected, with the overwhelming majority of the simulated trajectories safely avoiding the obstacles (red regions in the figure) and terminating in the goal (green region). A vector field associated with the policy is depicted in the third column of the figures. Notice that, the actions returned by our synthesis method intuitively align with the reach-avoid specification, that is, synthesized actions near the obstacle and out-of-bounds are aligned with moving away from such unsafe states. Exceptions to this are represented by locations where the agent is already unsafe and that, as such, are not fully explored during the BNN learning phase (e.g., the lower triangles in the Zigzag scenario), locations where two directions are equally optimal (e.g., in the top right corner of the v1 environment) and locations which are not along any feasibly optimal path (e.g., the lower right corner of v2) and as such are not accounted by the BNN learning.

In Table 1 we compare the certification results of the synthesized policy against the initial learned policy. As the synthesized policy is computed by improving on the latter, we expect the former to outperform the learned policy in terms of the guarantees obtained. This is in fact confirmed and quantified by the results of Table 1, which lists, for each of the three environments, the average reach-avoid probability estimated over 500 trajectories, the average certification lower bound across the state space, and the certification coverage (i.e., the proportion of states where our algorithm returns a non-zero probability lower bound). This notion of coverage only requires a state to be certified with a probability above 0, and so it is most informative when evaluated together with the average lower-bound and visual inspection of Figure 2 and Figure 3.

Indeed, the synthesized policy significantly improves on the certification guarantees given by the learned policy, and consistently so across the three environments analysed, with the lower bound improving by a factor of roughly 3.5. This considerable improvement is to be expected as worst-case guarantees can be poor for deep learning systems that are not trained with specific safety objectives(Mirman et al., 2018; Gowal et al., 2018; Wicker et al., 2021a). In particular, for both the V1 and Zigzag case studies, we observe that the average lower bound jumps from roughly 0.2 to greater than 0.7. Moreover, the most significant improvements are obtained in the most challenging case, i.e., the Zigzag environment, with the certification coverage increasing of a 4.75 factor. Interestingly, also the average model performance increases for the synthesized models. Intuitively this occurs because while in the learning of the initial policy passing through the obstacle is only pe
nalised by a continuous factor, the synthesized policy strives to rigorously enforce safety across the BNN posterior. A visual representation of these results is provided in the last column of Figure 2 for the learned policy and in Figure 3 for the synthesized policy. We note that the uncertainty maps in these figures are identical as the BNN model is not changed, only the policy.

### On the Effect of Approximate Inference

The results provided so far have been generated with BNN dynamical models learned via HMC training. However, different inference methods produce different approximations of the BNN posterior thus leading to different dynamics approximations and hence synthesized policies.

Table 2 and the plots in Figures 4 and 5 analyse the effect of approximate inference on both learned and synthesized policies, comparing results obtained by HMC with those obtained by VI training on the v1 scenario. We notice that also in the case of VI the synthesized policy significantly improves on the initial policy. Interestingly, the certification results over the learned policy for VI are higher than those obtained for HMC, but the results are comparable for the synthesized policies. In fact, it is known in the literature that VI tends to under-estimate uncertainty (Myshkov and Julier, 2016;

\begin{table}
\begin{tabular}{|l|c|c|c|} \multicolumn{4}{c}{_Learned Policy_}} \\ \cline{2-4}  & Performance & Avg. Lower Bound & Cert. Coverage \\ \hline V1 & 0.789 & 0.212 & 0.639 \\ \hline V2 & 0.805 & 0.192 & 0.484 \\ \hline Zigzag & 0.815 & 0.189 & 0.193 \\ \hline \end{tabular}
\end{table}
Table 1: Certification comparisons for learned and synthesized policy across the three environments. **Performance** indicates the proportion of simulated trajectory that respect the reach-avoid specification. **Avg. Lower Bound** is the mean certification probability across all states. **Cert. Coverage** is the proportion of states that we are able to certify (i.e., with a non-zero lower bound for the reach-avoid probability).
Michelmore et al., 2020) and is more susceptible to model misspecification (Masegosa, 2020). As such, being probabilistic, the bound obtained is tighter for VI where the uncertainty is lower than that of HMC which provides a more conservative representation of the agent dynamics. For example, we see in the first two rows of Table 2 that the average lower bound achieved for the variational inference posterior is 1.88 times higher than the bounds for HMC posterior. However, our synthesis method reduces this gap between HMC and VI, while still accounting for the higher uncertainty of the former, and hence the more conservative guarantees.

While HMC approximates the posterior by relying on a Monte Carlo estimate of it, VI is a gradient-based technique, where the number of training epochs (i.e., the number of full sweeps through the dataset) is a key hyperparameter. We thus analyse the effect of training epochs in the quality of the dynamics obtained in Figure 6, along with the effect on the synthesized policies and the certificates obtained for such policies. The left plot of the figure shows a set of predicted trajectories over a 10 time-step horizon for a varying number of training epochs, with the ground truth behaviour highlighted in red. The BNN trajectories are color-coded based on the number of epochs each dynamics model was trained for. In yellow, we see that the BNN which has only been trained for 10 epochs displays considerable error in its iterative predictions. This is reduced considerably for a model trained for 50 epochs,

\begin{table}
\begin{tabular}{|c|c|c|c|} \multicolumn{4}{c}{_Learned Policy_}} \\ \cline{2-4} \multicolumn{1}{c|}{} & Performance & Avg. Lower Bound & Coverage \\ \hline Var. Inference & 0.832 & 0.399 & 0.696 \\ \hline Ham. Monte Carlo & 0.789 & 0.212 & 0.639 \\ \hline \end{tabular}
\end{table}
Table 2: Certification comparisons for learned and synthesized policy between VI and HMC BNN learning on obstacle layout V1. **Performance** indicates the proportion of simulated trajectory that respect the reach-avoid specification. **Avg. Lower Bound** is the mean certification probability across all states. **Cert. Coverage** is the proportion of states that we are able to certify with non-zero probability.
but the cumulative error after 10 epochs is still considerable. Finally, as expected, for models trained for 250 and 1500 epochs we empirically observe a trend toward convergence to the ground truth.

We notice that the policy and certifications directly reflect the quality of the approximation. In fact, as we increase the model fit, we see that there are significant improvements in both the intuitive behavior of the synthesized policy as well as the resulting guarantees we are able to compute.

### Depth Experiments

In this section, we evaluate how our method performs when we vary the depth of the BNN dynamics model considered. In Figure 7, we plot the certified reach-avoid probabilities for a learned policy and a one, two, and three-layer BNN dynamics model where each layer has a width of 12 neurons. Similar architectures are found in the BNNs studied in recent related work (Lechner et al., 2021). Other than the depth of the BNN, all the other variables in the experiment are held equal (e.g., number of episodes during

Figure 4: **Top Row:** Visualization of the learned system using HMC to approximately infer dynamics. **Bottom Row:** Visualization of the learned system using VI to approximately infer dynamics. We highlight that the VI approximation displays a 5 to 10 times reduction in epistemic uncertainty.
learning, discretization of state-space, and number of BNN samples considered for the lower bound). The learning procedure results in BNN models with roughly equivalent losses and in policies that are qualitatively similar, see Appendix A.3 for further visualizations. Given that the key factors of the system have been held equal, we notice a decrease in our certified lower bound as depth increases. Specifically, the average lower bound for the one-layer model is 0.811, for the two-layer model it is 0.763, and for the three-layer model it decreases further to 0.621. Figure 7 clearly demonstrates that as the BNN dynamics model becomes deeper, then our lower bound becomes more conservative. This finding is consistent with existing results in certification of BNNs (Wicker et al., 2020; Berrada et al., 2021) and DNNs (Gowal et al., 2018; Mirman et al., 2018). We note, however, that when the verification parameters are refined, i.e., more samples from the BNN are taken, we are able to certify the three-layer system with an average certified lower-bound of 0.867 (see Appendix A.3). These additional BNN samples increase the runtime of our certification procedure by 1.5 times.

Figure 5: **Top Row:** Visualization of the synthesized policy and its performance based on the HMC dynamics model. **Bottom Row:** Visualization of the synthesized policy and its performance based on the VI dynamics model.
## 7 Related Work

Certification of machine learning models is a rapidly growing area (Gehr et al., 2018; Katz et al., 2017; Gowal et al., 2018; Wicker et al., 2022). While most of these methods have been designed for deterministic NNs, recently safety analysis of Bayesian machine learning models has been studied both for Gaussian processes (GPs) (Grosse et al., 2017; Cardelli et al., 2019; Blaas et al., 2020) and BNNs (Athalye et al., 2018; Cardelli et al., 2019; Wicker et al., 2020), including methods for adversarial training (Liu et al., 2019; Wicker et al., 2021). The above works, however, focus exclusively on the input-output behaviour of the models, that is, can only reason about static properties. Conversely, the problem we tackle in this work has additional complexity, as we aim to formally reason about iterative predictions, i.e., trajectory-level behaviour of a BNN interacting in a closed loop with a controller.

Iterative predictions have been widely studied for Gaussian processes (Gi

Figure 6: Analysis for number of training epochs used in performing VI training on the V1 environment. Left: sample of 10-step agent trajectories obtained with BNNs trained with VI and different number of epochs (red: ground truth trajectory). Right: synthesis and certification results for a selection of training epochs.
rard et al., 2003) and safety guarantees have been proposed in this setting in the context of model-based RL with GPs (Jackson et al., 2020; Polymenakos et al., 2019; Berkenkamp et al., 2017, 2016). However, all these works are specific to GPs and cannot be extended to BNNs, whose posterior predictive distribution is intractable and non-Gaussian even for the more commonly employed approximate Bayesian inference methods (Neal, 2012). Recently, iterative prediction of _neural network dynamic models_ have been studied (Wei and Liu, 2021; Adams et al., 2022) and methods to certify these models against temporal logic formulae have been derived (Adams et al., 2022). However, these works only focus on standard (i.e., non-Bayesian) neural networks with additive Gaussian noise. Closed-loop systems with known (deterministic) models and control policies modelled as BNNs are considered in (Lechner et al., 2021). In contrast with our work, Lechner et al. (2021) can only support deterministic models without noisy dynamics, only focus on the safety verification problem, and are limited to BNN posterior with unimodal weight distribution.

Various recent works consider verification or synthesis of RL schemes against reachability specifications (Sun et al., 2019; Konighofer et al., 2020; Bacci and Parker, 2020). None of these approaches, however, support both continuous state-action spaces and probabilistic models, as in this work. Continuous action spaces are supported in (Hasanbeig et al., 2019), where the authors provide RL schemes for the synthesis of policies maximising given temporal requirements, which is also extended to continuous state- and

Figure 7: We vary the depth of the BNN used to learn the dynamics and observe its effects on our certified safety probabilities over the first half of the Puck-V1 state-space. From left to right we plot the lower-bound reach-avoid probabilities for a one-layer BNN dynamics model, a two-layer BNN dynamics model, and a three-layer BNN dynamics model.
action-spaces in (Hasanbeig et al., 2020). However, the guarantees resulting from these model-free algorithms are asymptotic and thus of a different nature than those in this work. The work of Haesaert et al. (2017) integrates Bayesian inference and formal verification over control models, additionally proposing strategy synthesis approaches for active learning (Wijesuriya and Abate, 2019). In contrast to our paper these works do not support unknown noisy models learned via BNNs.

A related line of work concerns the synthesis of runtime monitors for predicting the safety of the policy's actions and, if necessary, correct them with fail-safe actions (Alshiekh et al., 2018; Avni et al., 2019; Bouton et al., 2019; Fulton and Platzer, 2019; Phan et al., 2020). These approaches, however, do not support continuous state-action spaces or require some form of ground-truth mechanistic model for safety verification (as opposed to our data-driven BNN models).

## 8 Conclusions

In this paper, we considered the problem of computing the probability of time-bounded reach-avoid specifications for dynamic models described by iterative predictions of BNNs. We developed methods and algorithms to compute a lower bound of this reach-avoid probability. Additionally, relying on techniques from dynamic programming and non-convex optimization, we synthesized certified controllers that maximize probabilistic reach-avoid. In a set of experiments, we showed that our framework enables certification of strategies on non-trivial control tasks. A future research direction will be to investigate techniques to enhance the scalability of our methods so that they can be applied to state-of-the-art reinforcement learning environments. However, we emphasise that the benchmark considered in this work remains a challenging one for certification purposes, due to both the non-linearity and stochasticity of the models, and the sequential, multi-step dependency of the predictions. Thus, this paper makes an important step toward the application of BNNs in safety-critical scenarios.

## Acknowledgements

This project received funding from the ERC under the European Union's Horizon 2020 research and innovation programme (FUN2MODEL, grant agreement No. 834115).