# Stochastic Control Barrier Functions with Bayesian Inference for Unknown
Stochastic Differential Equations

Chuanzheng Wang
Department of Applied Mathematics
University of Waterloo, Canada
&Yiming Meng
Grainger College of Engineering
University of Illinois Urbana-Champaign
&Jun Liu
Department of Applied Mathematics
University of Waterloo, Canada
&Stephen Smith
Electrical and Computer Engineering
University of Waterloo, Canada

###### Abstract

Control barrier functions are widely used to synthesize safety-critical
controls. However, the presence of Gaussian-type noise in dynamical systems
can generate unbounded signals and potentially result in severe consequences.
Although research has been conducted in the field of safety-critical control
for stochastic systems, in many real-world scenarios, we do not have precise
knowledge about the stochastic dynamics. In this paper, we delve into the
safety-critical control for stochastic systems where both the drift and
diffusion components are unknown. We employ Bayesian inference as a data-
driven approach to approximate the system. To be more specific, we utilize
Bayesian linear regression along with the central limit theorem to estimate
the drift term, and employ Bayesian inference to approximate the diffusion
term. Through simulations, we verify our findings by applying them to a
nonlinear dynamical model and an adaptive cruise control model.

For some real-world control problems, safety-critical control must be used in
order to prevent severe consequences. It requires not only achieving control
objectives, but also providing control actions with guaranteed safety [1].
Hence incorporating safety criteria is of great importance in practice when
designing controllers. These requirements need to be satisfied in practice
both with or without noise and disturbance. Since the notion of safety-
critical control was first introduced in [2], there has been extensive
research in safety verification problems, e.g., using discrete approximations
[3] and computation of reachable sets [4].

Recently, control barrier functions (CBFs) have been widely used to deal with
safety-critical control [5]. Quadratic programming (QP) problems are used to
solve safety-critical control with constraints from CBFs, together with
control Lyapunov functions (CLFs) for achieving stability objectives [6]. The
authors show that the safety criteria can be transformed into linear
constraints of the QP problems. By taking the derivative of the CBF, the
control inputs can be treated as the decision variables of the QP problem so
that we can find a sequence of actions that can guarantee safe trajectories.
The authors in [7] then show that finding safe control inputs by solving QP
problems can be extended to an arbitrary number of constraints and any nominal
control law. Consequently, safety-critical control using QP with CBF
conditions have been applied in a wide range of applications such as lane
keeping [8] and obstacle avoidance [9]. However, for many applications in
robotics [10], the derivative of the CBFs is not dependent on the control
input, and thus the control action can not be directly solved using QP. To
address this issue, CBFs have been extended to exponential control barrier
functions (ECBF) to handle high relative degree constraints using input-output
linearization [11] and in [12] the authors propose a more general form of
high-order control barrier functions (HOCBFs).

However, in practice, models that are used to design controllers are imperfect
and this may lead to unsafe or even dangerous behavior of the systems.
Consequently, designing controllers considering uncertainty is important in
practical applications. In [13] and [14], a bounded disturbance is considered
for the model, in which the time derivative of the barrier function is
separated into the time derivative of the nominal barrier function and a
remainder that can be approximated using neural networks. For systems driven
by Gaussian-type noise, stochastic differential equations (SDEs) are usually
used to characterize the effect of randomness. Studies for stochastic
stability of diffusion-type stochastic differential equations have seen a
variety of applications in verifying probabilistic quantification of safe set
invariance [15]. Control barrier functions for stochastic systems have also
been studied in recent years. The authors in [16] applied the strong set-
invariance certificate from [17] to high-order stochastic control systems
using stochastic reciprocal control barrier functions (SRCBFs) and [18]
investigates the worst-case safety verification utilizing stochastic zeroing
control barrier functions (SZCBFs) regardless of the magnitude of noise. In
[19], the authors proposed the stochastic control barrier function (SCBF) with
milder conditions at the cost of sacrificing the almost sure safety.

However, in some practical scenarios, we do not have precise information about
the system and in this case, we cannot directly apply (stochastic) CBFs to
control the system safely. One way of handling these unknown dynamics is to
identify the model first and then applied control barrier certificates to
guarantee safety. The Matrix-Variate Gaussian process model is used to learn
nonlinear control-affined discrete systems and the learned model is
incorporated into a multi-agent CBF framework as in [20]. A two-stage solution
is proposed in [21] for unknown systems, where Gaussian learning is used to
learn the model first and then use CBF for safety guarantee. In [22] and [23],
Bayesian learning is used to approximate the unknown model with barrier
certificates for safety-critical control. However, all of these methods only
focus on deterministic systems. Data-driven methods with scenario convex
programming (SCP) is also studied in guaranteeing safety of unknown dynamical
systems. In [24] and [25], random trajectories are sampled with barrier
certificates as finite number of constraints in the optimization problems for
unknown deterministic system and stochastic systems, respectively. However,
such methods only apply to discrete systems. Formal methods with Gaussian
process learning is also studied in [26] for discrete systems as well. In
[27], the authors assume that the diffusion term of the SDE is unknown and use
sampling method with central limit theorem to estimate the generator the such
SDEs for safety-critical control. However, such sampling-based method is very
time-inefficient when the relative degree is higher than one. On the other
hand, data-driven methods with Bayesian learning has been widely used to
identify unknown stochastic differential equations recently as in [28] and
[29]. And recently, the authors in [30] proposed sparse Bayesian learning
architecture to handle insufficient data resources or prior information for
unknown SDEs. Motivated by such results, we consider safety-critical control
problems for unknown stochastic systems. To the best of our knowledge, this is
the first paper to address safety-control problems for SDEs with fully unknown
drift and diffusion terms.

In this paper, we extend our work in [19] for stochastic control barrier
functions (SCBFs) to unknown SDEs. We use Bayesian learning to estimate the
system and then use SCBFs on the identified system. The proposed method relies
on data-driven method such that we will estimate both drift and diffusion
terms based our observations of the system. We also compare the simulation
results both in safety ratio and running time against our previous study in
[19] and [27]. As a result, the main contribution of this paper is as follows.

  * •

We propose the stochastic control barrier function with high relative degree
for safety-critical control of SDEs.

  * •

We propose a Bayesian learning framework for safety-critical control of
unknown SDEs

  * •

We validate our results and compare with our previous work in [19] and [27]
for safety ratio and running time.

##  1 Preliminary and Problem Definition

###  1.1 System Description

Given a state space $\mathcal{X}\subseteq\mathbb{R}^{n}$, a (compact) set of
control values $\mathcal{U}\subset\mathbb{R}^{p}$, consider the following form
of perturbed SDE:

$dX_{t}=f(X_{t})dt+g(X_{t})\mathfrak{u}_{t}dt+b(X_{t})dW_{t},$ (1)

where $W$ is an $m$-dimensional standard Wiener process;
$f:\mathcal{X}\rightarrow\mathbb{R}^{n}$ is a nonlinear local Lipschitz
continuous vector field; $g,b:\mathcal{X}\rightarrow\mathbb{R}^{n\times m}$
are smooth mappings. The input $\mathfrak{u}:\mathbb{R}_{\geq
0}\rightarrow\mathcal{U}$ is a locally bounded control signal.

In brief, without loss of generality, we assume that a stochastic process
$X:=\\{X_{t}\\}_{t\in[0,\infty)}$ and a process of control values
$\mathfrak{u}:=\\{\mathfrak{u}_{t}\\}_{t\in[0,\infty)}$ are defined on some
(unknown) probability space where the noise is generated. Given any measurable
process $\mathfrak{u}$, the probability law of the joint process
$(X,\mathfrak{u}):=\\{(X_{t},\mathfrak{u}_{t})\\}_{t\in[0,\infty)}$ can be
determined on the measurable space
$((\mathcal{X}\times\mathcal{U})^{\otimes[0,\infty)},\mathcal{F},\mathbf{P})$,
where $\mathcal{F}$ is the Borel $\sigma$-algebra of
$(\mathcal{X}\times\mathcal{U})^{\otimes[0,\infty)}$. We also denote
$X^{\mathfrak{u}}$ by the controlled process if we emphasize on the state-
space marginal of $(X,\mathfrak{u})$.

######  Definition 1 ( Control strategy)

A control strategy is a set-valued measurable function

$\kappa:\mathcal{X}\rightarrow 2^{\mathcal{U}}.$ (2)

We consider the following memoryless set of non-randomized control signals
throughout this paper.

######  Definition 2 (State-dependent control)

We say that a control signal $\mathfrak{u}$ conforms to a control strategy
$\kappa$ for (4) if, given $\\{X_{t}=x\\}$ for any $x\in\mathcal{X}$, we have
that

$\mathfrak{u}_{t}\in\kappa(x),\quad\forall t\geq 0.$ (3)

We further denote the set of all control signals that confirm to $\kappa$ as
$\mathbf{u}_{\kappa}$.

In this view, for system (1), we only consider the weak sense of solutions.
The system admits a weak solution if there exists a filtered probability space
$(\Omega^{\dagger},\mathscr{F}^{\dagger},\\{\mathscr{F}^{\dagger}_{t}\\},\mathbb{P}^{\dagger})$
with a pair $(X,\mathfrak{u},W)$ of adapted stochastic processes, such that
$W$ is a Wiener process and $(X,\mathfrak{u})$ solves the controlled SDE (1).

Given a filtered probability space
$(\Omega,\mathcal{F},\\{\mathcal{F}_{t}\\},\mathbb{P})$ with a natural
filtration, a state space $\mathcal{X}\subseteq\mathbb{R}^{n}$, a (compact)
set of control values $\mathcal{U}\subset\mathbb{R}^{p}$, consider a
continuous-time stochastic process
$X:[0,\infty)\times\Omega\rightarrow\mathcal{X}$ that solves the SDE

$dX_{t}=(f(X_{t})+g(X_{t})u(t))dt+\sigma(X_{t})dW_{t},$ (4)

where $u:\mathbb{R}_{\geq 0}\rightarrow\mathcal{U}$ is a bounded measurable
control signal; $W$ is a $d$-dimensional standard
$\\{\mathcal{F}\\}_{t}$-Brownian motion;
$f:\mathcal{X}\rightarrow\mathbb{R}^{n}$ is a nonlinear vector field;
$g:\mathcal{X}\rightarrow\mathbb{R}^{n\times p}$ and
$\sigma:\mathcal{X}\rightarrow\mathbb{R}^{n\times d}$ are smooth mappings.

###### Assumption 1

We make the following assumptions on system (4) for the rest of this paper:

  * (i)

There is a $\xi\in\mathcal{X}$ such that $\mathbb{P}[X_{0}=\xi]=1$;

  * (ii)

The mappings $f,g,\sigma$ satisfy local Lipschitz continuity and a linear
growth condition.

######  Definition 3 (Strong solutions)

A stochastic process $X$ is said to be a strong solution to (4) if it
satisfies the following integral equation

$X_{t}=\xi+\int_{0}^{t}(f(X_{s})+g(X_{s})u(s))ds+\int_{0}^{t}\sigma(X_{s})dW_{s},$
(5)

where the stochastic integral is constructed based on the given Brownian
motion $W$.

######  Definition 4 (Infinitesimal generator of $X_{t}$)

Let $X$ be the strong solution to (4), the infinitesimal generator
$\mathcal{A}$ of $X_{t}$ is defined by

$\mathcal{A}h(x)=\lim\limits_{t\downarrow
0}\frac{\mathbb{E}^{x}[h(X_{t})]-h(x)}{t};\;\;x\in\mathbb{R}^{n},$ (6)

where $h:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is in a set
$\mathcal{D}(\mathcal{A})$ (called the domain of the operator $\mathcal{A}$)
of functions such that the limit exists at $x$.

######  Proposition 1 (Dynkin)

Let $X$ solve (4). If $h\in C_{0}^{2}(\mathbb{R}^{n})$ then
$h\in\mathcal{D}(\mathcal{A})$ and

$\mathcal{A}h(x)=\frac{\partial h}{\partial
x}(f(x)+g(x)u(t))+\frac{1}{2}\sum\limits_{i,j}\left(\sigma\sigma^{T}\right)_{i,j}(x)\frac{\partial^{2}h}{\partial
x_{i}\partial x_{j}}.$ (7)

###### Remark 1

The solution $X$ to (4) is right continuous and satisfies strong Markov
properties, and for any finite stopping time $\tau$ and $h\in
C_{0}^{2}(\mathbb{R}^{n})$, we have the following Dynkin’s formula

$\mathbb{E}^{\xi}[h(X_{\tau})]=h(\xi)+\mathbb{E}^{\xi}\left[\int_{0}^{\tau}\mathcal{A}h(X_{s})ds\right],$

and therefore

$h(X_{\tau})=h(\xi)+\int_{0}^{t}\mathcal{A}h(X_{s})ds+\int_{0}^{\tau}\nabla_{x}h(X_{s})dW_{s}.$

The above is an analogue of the evolution of $h$ along trajectories

$h(x(t))=h(\xi)+\int_{0}^{t}[L_{f}h(x(s))+L_{g}h(x(s))u(s)]ds$

driven by deterministic dynamics $\dot{x}=f(x)+g(x)u$, where
$L_{f}h=\nabla_{x}h(x)\cdot f(x)$ and $L_{g}h=\nabla_{x}h(x)\cdot g(x)$.

###  1.2 Set Invariance and Control

In deterministic settings, a set $\mathcal{C}\subseteq\mathcal{X}$ is said to
be invariant for a dynamical system $\dot{x}=f(x)$ if, for all
$x(0)\in\mathcal{C}$, the solution $x(t)$ is well defined and
$x(t)\in\mathcal{C}$ for all $t\geq 0$. As for stochastic analogies, we have
the following probabilistic characterization of set invariance.

######  Definition 5 (Probabilistic set invariance)

Let $X$ be a stochastic process. A set $\mathcal{C}\subset\mathcal{X}$ is said
to be invariant w.r.t. a tuple $(x,T,p)$ for $X$, where $x\in\mathcal{C}$,
$T\geq 0$, and $p\in[0,1]$, if $X_{0}=x$ a.s. implies

$\mathbb{P}_{x}[X_{t}\in\mathcal{C},\;0\leq t\leq T]\geq p.$ (8)

Moreover, if $\mathcal{C}\subset\mathcal{X}$ is invariant w.r.t. $(x,T,1)$ for
all $x\in\mathcal{C}$ and $T\geq 0$, then $\mathcal{C}$ is strongly invariant
for $X$.

For stochastic dynamical systems with controls such as system (4), we would
like to define similar probabilistic set invariance property for the
controlled processes. Before that, we first define the following concepts.

######  Definition 6 (Control strategy)

A control strategy is a set-valued function

$\kappa:\mathcal{X}\rightarrow 2^{\mathcal{U}}.$ (9)

######  Definition 7 (Controlled probabilistic invariance)

Given system (4) and a set of control signals $\mathbf{u}$, a set
$\mathcal{C}\subset\mathcal{X}$ is said to be controlled invariant under
$\mathbf{u}$ w.r.t. a tuple $(x,T,p)$ for system (4), if for all
$u\in\mathbf{u}$, $\mathcal{C}$ is invariant w.r.t. $(x,T,p)$ for $X$, where
$X$ is the solution to (4) with $u$ as input.

Similarly, $\mathcal{C}\subset\mathcal{X}$ is strongly controlled invariant
under $\mathbf{u}$ if $\mathcal{C}\subset\mathcal{X}$ is controlled invariant
under $\mathbf{u}$ w.r.t. $(x,T,1)$ for all $x\in\mathcal{C}$ and $T\geq 0$.

###  1.3 Problem Formulation

For the rest of this paper, we consider a safe set of the form

$\mathcal{C}:=\\{x\in\mathcal{X}:h(x)\geq 0\\},$ (10)

where $h\in C^{2}(\mathbb{R}^{n})$. We also define the boundary and interior
of $\mathcal{C}$ explicitly as below

$\partial\mathcal{C}:=\\{x\in\mathcal{X}:h(x)=0\\},$ (11)
$\mathcal{C}^{\circ}:=\\{x\in\mathcal{X}:h(x)>0\\}.$ (12)

The objective of this paper is to control the stochastic system (4) with
unknown drift and diffusion terms to stay inside the safe set. The problem is
defined as follows.

###### Problem 1

Given system (4) with $f$, $g$ and $\sigma$ unknown, a compact set
$\mathcal{C}\subseteq\mathcal{X}$ defined by (10), a point
$x\in\mathcal{C}^{\circ}$, and a $p\in[0,1]$, design a (deterministic) control
strategy $\kappa$ such that under $\mathbf{u}_{\kappa}$, the interior
$\mathcal{C}^{\circ}$ is controlled $p$-invariant for the resulting solutions
to (4).

###### Remark 2

In our previous paper [27], we assume that we know the drift term of the
system and only the diffusion term $\sigma$ is unknown. In this paper, we
extend our assumption that the system is fully unknown.

##  2 Safe-critical Control Design via Barrier Functions

Before proposing our stochastic barrier certificates that can be used to
design a control strategy $\kappa$ for Problem 1, it is necessary to review
(stochastic) control barrier functions to interpret (probabilistic) set
invariance. Note that we consider the safe set as constructed in (10), where
the function $h$ is given a priori.

###  2.1 Stochastic Reciprocal and Zeroing Barrier Functions

Similar to the terminology for deterministic cases [8], we introduce the
construction of stochastic control barrier functions as follows.

######  Definition 8 (SRCBF)

A function $B:\mathcal{C}^{\circ}\rightarrow\mathbb{R}$ is called a stochastic
reciprocal control barrier function (SRCBF) for system (4) if
$B\in\mathcal{D}(\mathcal{A})$ and satisfies the following properties:

  * (i)

there exist class-$\mathcal{K}$ functions $\alpha_{1},\alpha_{2}$ such that
for all $x\in\mathcal{X}$ we have

$\frac{1}{\alpha_{1}(h(x))}\leq B(x)\leq\frac{1}{\alpha_{2}(h(x))};$ (13)

  * (ii)

there exists a class-$\mathcal{K}$ function $\alpha_{3}$ such that

$\inf\limits_{u\in\mathcal{U}}[\mathcal{A}B(x)-\alpha_{3}(h(x))]\leq 0.$ (14)

We refer to the control strategy generated by (14) as

$\varrho(x):=\\{u\in\mathcal{U}:\mathcal{A}B(x)-\alpha_{3}(h(x))\leq 0\\}$
(15)

and the corresponding control constraint as $\mathbf{u}_{\varrho}$ (see in
Definition 2).

######  Proposition 2 ([17])

Suppose that there exists an SRCBF for system (4). If
$u(t)\in\mathbf{u}_{\upsilon}$, then for all $t\geq 0$ and
$X_{0}=\xi\in\mathcal{C}^{\circ}$, we have
$\mathbb{P}_{\xi}[X_{t}\in\mathcal{C}^{\circ}]=1$ for all $t\geq 0$.

######  Definition 9 (SZCBF)

A function $B:\mathcal{C}\rightarrow\mathbb{R}$ is called a stochastic zeroing
control barrier function (SZCBF) for system (4) if
$B\in\mathcal{D}(\mathcal{A})$ and

  1. (i)

$B(x)\geq 0$ for all $x\in\mathcal{C}$;

  2. (ii)

$B(x)<0$ for all $x\notin\mathcal{C}$;

  3. (iii)

there exists an extended $\mathcal{K}_{\infty}$ function $\alpha$ such that

$\sup_{u\in\mathcal{U}}{[\mathcal{A}B(x)}+\alpha(B(x))]\geq 0.$ (16)

We refer the control strategy generated by (16) as

$\varkappa(x):=\\{u\in\mathcal{U}:{\mathcal{A}B(x)}+\alpha(B(x))\geq 0\\}$
(17)

and the corresponding set of constrained control signals as
$\mathbf{u}_{\varkappa}$.

######  Proposition 3 (Worst-case probabilistic quantification)

Suppose the mapping $h$ is an SZCBF with linear function $kx$ as the
class-$\mathcal{K}$ function (where $k>0$), and the control strategy as
$\varkappa(x)=\\{u\in\mathcal{U}:\mathcal{A}h(x)+kh(x)\geq 0\\}$. Let
$c=\sup_{x\in\mathcal{C}}h(x)$ and $X_{0}=\xi\in\mathcal{C}^{\circ}$, then
under any $u\in\mathbf{u}_{\varkappa}$ we have the following worst-case
probability estimation:

$\mathbb{P}_{\xi}\left[X_{t}\in\mathcal{C}^{\circ},\;0\leq t\leq
T\right]\geq\left(\frac{h(\xi)}{c}\right)e^{-cT}.$ (18)

###### Proof 1

Let $s=c-h(\xi)$ and $V(x)=c-h(x)$, then $V(x)\in[0,c]$ for all
$x\in\mathcal{C}$. It is clear that $\mathcal{A}V(x)=-\mathcal{A}h(x)$. For
$u(t)\in\mathbf{u}_{\varkappa}$ for all $t\in[0,T]$, we have

$\mathcal{A}V(x)\leq-kV(x)+kc.$

By [15, Theorem 3.1],

$\begin{split}\mathbb{P}_{\xi}\left[\sup\limits_{t\in[0,T]}V(X_{t})\geq
c\right]\leq 1-\left(1-\frac{s}{c}\right)e^{-cT}.\end{split}$ (19)

The result follows directly after this.

###  2.2 Motivation of Stochastic Control Barrier Functions

The authors in [16] constructed high-order SRCBF and have found the sufficient
conditions to guarantee pathwise set invariance with probability 1. While the
results seem strong, they come with significant costs. At the safety boundary,
the control inputs need to be unbounded in order to guarantee safety. (A
motivated example is illustrated in [19]). On the other hand, the synthesis of
controller for a high-order system via SZCBF is with mild constraints. The
trade-off is that the probability estimation of set invariance is of low
quality. As a result, we propose high-order stochastic control barrier
functions in this paper in order to reduce the high control efforts and
improve the worst-case quantification.

###  2.3 High-order Stochastic Control Barrier Functions

To obtain non-vanishing worst-case probability estimation (compared to SZCBF),
we propose a safety certificate via a stochastic Lyapunov-like control barrier
function [15].

######  Definition 10 (Stochastic control barrier functions)

A continuously differentiable function $B:\mathbb{R}^{n}\rightarrow\mathbb{R}$
is said to be a stochastic control barrier function (SCBF) if
$B\in\mathcal{D}(\mathcal{A})$ and the following conditions are satisfied:

  1. (i)

$B(x)\geq 0$ for all $x\in\mathcal{C}$;

  2. (ii)

$B(x)<0$ for all $x\notin\mathcal{C}$;

  3. (iii)

$\sup\limits_{u\in\mathcal{U}}\mathcal{A}B(x)\geq 0$.

We refer the control strategy generated by (iii) as

$\upsilon(x):=\\{u\in\mathcal{U}:\mathcal{A}B(x)\geq 0\\}$ (20)

and the corresponding set of constrained control signals as
$\mathbf{u}_{\upsilon}$.

###### Remark 3

Condition (iii) of the above definition is an analogue of
$\sup\limits_{u\in\mathcal{U}}[L_{f}B(x)+L_{g}B(x)u]\geq 0$ for the
deterministic settings. The consequence is such that
$\mathbb{E}^{\xi}[B(x)]\geq B(x)$ for all $x\in\mathcal{C}$. A relaxation of
condition (iii) is given in the deterministic case such that the set
invariance can still be guaranteed [31],

$\sup\limits_{u\in\mathcal{U}}[L_{f}B(x)+L_{g}B(x)u]\geq-\alpha(B(x)),$

where $\alpha$ is a class-$\mathcal{K}$ function. If $\alpha(x)=kx$ where
$k>0$, under the stochastic settings, the condition formulates an SZCBF and
provides a much weaker quantitative estimation of the lower bound of
satisfaction probability. In comparison with SZCBF, we provide the worst-case
quantification in the following proposition.

###### Proposition 4

Suppose the mapping $h$ is an SCBF with the corresponding control strategy
$\upsilon(x)$. Let $c=\sup_{x\in\mathcal{C}}h(x)$ and
$X_{0}=\xi\in\mathcal{C}^{\circ}$, then under the set of constrained control
signals $\mathbf{u}_{\upsilon}$, we have the following worst-case probability
estimation:

$\mathbb{P}_{\xi}\left[X_{t}\in\mathcal{C}^{\circ},\;0\leq
t<\infty\right]\geq\frac{h(\xi)}{c}.$

###### Proof 2

Let $V=c-h(x)$, then $V(x)\geq 0$ for all $x\in\mathcal{X}$ and
$\mathcal{A}V\leq 0$. Let $s=c-h(\xi)$ then $s\leq c$ by definition. The
result for every finite time interval $t\in[0,T]$ is followed by [15, Lemma
2.1],

$\mathbb{P}_{\xi}\left[X_{t}\in\mathcal{C}^{\circ},\;0\leq t<T\right]\geq
1-\frac{s}{c}.$

The result follows by letting $T\rightarrow\infty$.

###### Definition 11

A function $B:\mathcal{X}\rightarrow\mathbb{R}$ is called a stochastic control
barrier function with relative degree $r$ for system (4) if
$B\in\mathcal{D}(\mathcal{A}^{r})$, and
$\mathcal{A}\circ\mathcal{A}^{r-1}h(x)\neq 0$ as well as
$\mathcal{A}\circ\mathcal{A}^{j-1}h(x)=0$ for $j=1,2,\ldots,r-1$ and
$x\in\mathcal{C}$.

If the system (4) is an $r^{\text{th}}$-order stochastic control system, to
steer the process $X$ to satisfy probabilistic set invariance w.r.t.
$\mathcal{C}$, we recast the mapping $h$ as an SCBF with relative degree $r$.
For $h\in\mathcal{D}(\mathcal{A}^{r})$, we define a series of functions
$b_{0},b_{j}:\mathcal{X}\to\mathbb{R}$ such that for each $j=1,2,\dots,r$
$b_{0},b_{j}\in\mathcal{D}(\mathcal{A})$ and

$\begin{split}b_{0}(x)&=h(x),\\\
b_{j}(x)&=\mathcal{A}\circ\mathcal{A}^{j-1}b_{0}(x).\\\ \end{split}$ (21)

We further define the corresponding superlevel sets $\mathcal{C}_{j}$ for
$j=1,2,\dots,r$ as

$\begin{split}\mathcal{C}_{j}&=\\{x\in\mathbb{R}^{n}:b_{j}(x)\geq 0\\}.\\\
\end{split}$ (22)

###### Theorem 1

If the mapping $h$ is an SCBF with relative degree $r$, the corresponding
control strategy is given as
$\upsilon(x)=\\{u\in\mathcal{U}:\mathcal{A}^{r}h(x)\geq 0\\}$. Let
$c_{j}=:\sup_{x\in\mathcal{C}_{j}}b_{j}(x)$ for each $j=0,1,...,r$ and
$X_{0}=\xi\in\bigcap_{j=0}^{r}\mathcal{C}_{j}^{\circ}$. Then under the set of
constrained control signals $\mathbf{u}_{\upsilon}$, we have the following
worst-case probability estimation:

$\mathbb{P}_{\xi}[X_{t}\in\mathcal{C}^{\circ},\;0\leq
t<\infty]\geq\prod\limits_{j=0}^{r-1}\frac{b_{j}(\xi)}{c_{j}}.$

###### Proof 3

We introduce the notations
$p_{j}:=\mathbb{P}_{\xi}[X_{t}\in\mathcal{C}_{j}^{\circ},\;0\leq t<\infty]$ ,
$\hat{p}_{j}:=\mathbb{P}_{\xi}[X_{t}\in\mathcal{C}_{j}^{\circ},\;0\leq
t<\infty\;|\;\mathcal{A}b_{j}\geq 0]$.

The control signal $u(t)\in\\{u\in\mathcal{U}:\mathcal{A}^{r}b(x)\geq 0\\}$
for all $t\geq 0$ provides $\mathcal{A}b_{r-1}\geq 0$. By Proposition 4,

$\begin{split}p_{r-1}&=\mathbb{P}_{\xi}[X_{t}\in\mathcal{C}_{r-1}^{\circ},\;0\leq
t<\infty]\\\ &=\mathbb{P}_{\xi}[X_{t}\in\mathcal{C}_{r-1}^{\circ},\;0\leq
t<\infty\;|\;\mathcal{A}b_{r-1}\geq 0]\\\ &=\hat{p}_{r-1}.\end{split}$ (23)

For $j=0,1,...,r-2$, and $0\leq t<\infty$, we have the following recursion:

$\begin{split}&p_{j}=\mathbb{P}_{\xi}[X_{t}\in\mathcal{C}_{j}^{\circ}]\\\
&=\mathbb{E}^{\xi}[\mathds{1}_{\\{X_{t}\in
C_{j}^{\circ}\\}}\mathds{1}_{\\{X_{t}\in
C_{j+1}^{\circ}\\}}]+\mathbb{E}^{\xi}[\mathds{1}_{\\{X_{t}\in
C_{j}^{\circ}\\}}\mathds{1}_{\\{X_{t}\notin C_{j+1}^{\circ}\\}}]\\\
&=\mathbb{P}_{\xi}[X_{t}\in C_{j}^{\circ}\;|\mathds{1}_{\\{X_{t}\in
C_{j+1}^{\circ}\\}}]\cdot\mathbb{P}_{\xi}[X_{t}\in C_{j+1}^{\circ}]\\\
&+\mathbb{P}_{\xi}[X_{t}\in C_{j}^{\circ}\;|\mathds{1}_{\\{X_{t}\notin
C_{j+1}^{\circ}\\}}]\cdot\mathbb{P}_{\xi}[X_{t}\notin
C_{j+1}^{\circ}],\end{split}$ (24)

where we have used shorthand notations
$\\{X_{t}\in\mathcal{C}_{j}\\}:=\\{X_{t}\in\mathcal{C}_{j},\;0\leq
t<\infty\\}$ and $\\{X_{t}\notin C_{j+1}\\}:=\\{X_{t}\notin
C_{j+1}\;\text{for\;some}\;0\leq t<\infty\\}$. Indeed, we have

$X_{t}=X_{t}\mathds{1}_{\\{X_{t}\in
C_{j+1}^{\circ}\\}}+X_{t}\mathds{1}_{\\{X_{t}\notin C_{j+1}^{\circ}\\}},$

then

$\mathbb{E}[X_{t}]=\mathbb{E}[X_{t}\mathds{1}_{\\{X_{t}\in
C_{j+1}^{\circ}\\}}]+\mathbb{E}[X_{t}\mathds{1}_{\\{X_{t}\notin
C_{j+1}^{\circ}\\}}],$

and (24) follows. Note that

$\begin{split}&\mathbb{P}_{\xi}[X_{t}\in
C_{j}^{\circ}\;|\mathds{1}_{\\{X_{t}\in C_{j+1}^{\circ}\\}}]\\\
&\geq\mathbb{P}_{\xi}[X_{t}\in C_{j}^{\circ}\;|\mathcal{A}b_{j}\geq
0]=\hat{p}_{j},\end{split}$ (25)

and therefore

$\mathbb{P}_{\xi}[X_{t}\in C_{j}^{\circ}\;|\mathds{1}_{\\{X_{t}\in
C_{j+1}^{\circ}\\}}]\cdot\mathbb{P}_{\xi}[X_{t}\in
C_{j+1}^{\circ}]\geq\hat{p}_{j}p_{j+1}.$

Now define stopping times $\tau_{j}=\inf\\{t:b_{j}(X_{t})\leq 0\\}$ for
$j=0,1,...,r-2$, then $b_{j}(X_{t\wedge\tau_{j}})\geq 0$ a.s.. In addition,

$X_{t\wedge\tau_{j}}=\mathds{1}_{\\{\tau_{j}\leq\tau_{j+1}\\}}X_{t\wedge\tau_{j}}+\mathds{1}_{\\{\tau_{j}>\tau_{j+1}\\}}X_{t\wedge\tau_{j}}$

Assume the worst scenario, which is for all $t\geq\tau_{j+1}$, we have
$\mathcal{A}b_{j}\leq 0$. On
$\\{\mathcal{A}b_{j}<0\\}\cap\\{\tau_{j}>\tau_{j+1}\\}$, we have
$b_{j}(X_{\tau_{j+1}})>0$ and
$\mathbb{E}^{\xi}[b_{j}(X_{t\wedge\tau_{j}})]\leq
b_{j}(X_{\tau_{j+1}})-\int_{\tau_{j+1}}^{t}\varepsilon(s)ds$ for some
$\varepsilon:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{>0}$ and
$t\geq\tau_{j+1}$. Therefore, the process $b_{j}(X_{t\wedge\tau_{j}})$ is a
nonnegative supermartingale and $\mathbb{P}_{\xi}[\sup_{T\leq
t<\infty}b_{j}(X_{t\wedge\tau_{j}})\geq\lambda]\leq\frac{b_{j}(X_{\tau_{j+1}})-\int_{\tau_{j+1}}^{T}\varepsilon(s)ds}{\lambda}$
by Doob’s supermartingale inequality. For any $\lambda>0$, we can find a
finite $T\geq\tau_{j+1}$ such that $\mathbb{P}_{\xi}[\sup_{T\leq
t<\infty}b_{j}(X_{t\wedge\tau_{j}})\geq\lambda]=0$. Since $\lambda$ is
arbitrarily selected, we must have $\mathbb{P}_{\xi}[\sup_{T\leq
t<\infty}b_{j}(X_{t\wedge\tau_{j}})>0]=0$, which means $\tau_{j}$ is triggered
within finite time. On the other hand, on
$\\{\mathcal{A}b_{j}<0\\}\cap\\{\tau_{j}\leq\tau_{j+1}\\}$, $\tau_{j}$ has
been already triggered. Therefore,
$\\{X_{t}\notin\mathcal{C}_{j}^{\circ}\;\text{for\;some}\;0\leq t<\infty\\}$
a.s. given $\\{\mathcal{A}b_{j}<0\\}$. Hence,

$\begin{split}&\mathbb{P}_{\xi}[X_{t}\in
C_{j}^{\circ}\;|\mathds{1}_{\\{X_{t}\notin C_{j+1}^{\circ}\\}}]\\\
&\geq\mathbb{P}_{\xi}[X_{t}\in C_{j}^{\circ}\;|\mathds{1}_{\\{X_{t}\notin
C_{j+1}^{\circ},\;\forall t\geq\tau_{j+1}\\}}]\\\
&\geq\mathbb{P}_{\xi}[X_{t}\in
C_{j}^{\circ}\;|\mathds{1}_{\\{\mathcal{A}b_{j}<0,\;\forall
t\geq\tau_{j+1}\\}}]=0\end{split}$ (26)

Combining the above, for $j=0,1,...,r-2$, we have

$p_{j}\geq\hat{p}_{j}p_{j+1},$

and ultimately $p_{0}\geq\prod\limits_{j=0}^{r-1}\hat{p}_{j}$.

###### Remark 4

The above result estimates the lower bound of the safety probability given the
constrained control signals $\mathbf{u}_{\upsilon}$. Based on recursion (24),
we can easily obtain the same result by dropping the last term. However, we
argued that under some extreme conditions the worst case may happen. Indeed,
we have assumed that $t\geq\tau_{j+1}\implies\mathcal{A}b_{j}\leq 0$. This
conservative assumption is made such that within finite time $X$ will cross
the boundary of each $\mathcal{C}_{j}$.

Another implicit condition may cause the worst-case lower bound as well, that
is when $\bigcup_{j=0}^{r-1}\\{\mathcal{A}b_{j}=0,\;0\leq t<\infty\\}$ is a
$\mathbb{P}_{\xi}$-null set. This, however, is practically possible since the
controller indirectly influences the value of $\mathcal{A}b_{j}$ for all
$j<r$, the strong invariance of the level set $\\{\mathcal{A}b_{j}=0\\}$ is
not guaranteed using QP scheme.

##  3 Data-driven method for safety-critical control of unknown SDEs

In this section, we will first address the process of identifying the system
with data-driven method using Bayesian inference and then, propose our QP-
based control framework for safety-critical control of the unknown SDE with
learned dynamics.

###  3.1 Data collection for drift terms

In order to collect training data to estimate $f(x)$ and $g(x)$, we sample an
initial state and let the system evolve under control inputs. More
specifically, given some $T>0$, we divide the time interval $[0,T]$ into a
sufficiently refined partition $\\{t_{0},t_{1},...,t_{N}\\}$ with
$N\in\mathbb{N}$ such that $0=t_{0}<t_{1}<\dots<t_{N}=T$. Define the index set
as $\mathcal{I}=\\{0,1,\dots,N\\}$. At each $i\in\mathcal{I}$, we simulate the
system with $K$ times from $x_{i}$ with two control inputs $u_{i,1}$ and
$u_{i,2}$ such that $u_{i,1}\neq u_{i,2}$. For each $u_{i,j}$ with
$j\in\\{1,2\\}$, we will get $x^{k}_{i,j}$ for $k\in\\{1,2,...,K\\}$. By
applying central limit theorem, we get

$dx_{i,j}=\sum_{k=1}^{K}\frac{[x_{i+1,j}-x_{i}|x_{i}]}{K}.$

Accordingly, we can calculate the target value for $f(x_{i})$ and $g(x_{i})$
as

$\begin{split}y_{f(x_{i})}&=\frac{dx_{i,1}\cdot u_{i,2}-dx_{i,2}\cdot
u_{i,2}}{(u_{i,2}-u_{i,1})(t_{i+1}-t_{i})}\\\
y_{g(x_{i})}&=\frac{dx_{i,1}-dx_{i,2}}{(u_{i,1}-u_{i,2})(t_{i+1}-t_{i})}\\\
\end{split}$ (27)

Finally we will get training data as

$\begin{split}\textbf{X}&=[x_{0},\dots,x_{i},\dots,x_{N}]^{T},\\\
\textbf{Y}_{f}&=[y_{f(x_{0})},\dots,y_{f(x_{i})},\dots,y_{f(x_{N})}]^{T}\\\
\textbf{Y}_{g}&=[y_{g(x_{0})},\dots,y_{g(x_{i})},\dots,y_{g(x_{N})}]^{T}\\\
\end{split}$ (28)

###  3.2 Identification of drift terms

We use Bayesian linear regression model [32] to identify the drift term of the
system. Define $\Phi$ as the base function and $\theta$ as the weight vector
where $\Phi\in\mathbb{R}^{N\times M}$ and $\theta\in\mathbb{R}^{M\times 1}$.
Here, $M$ is the number of base function. We can write the regression equation
as

$Y=\Phi\theta+\varepsilon$ (29)

The noise term $\varepsilon$ is assumed to follow a Gaussian distribution
$\mathcal{N}(0,\Psi)$ where $\Psi$ is a diagonal matrix with the $i$th element
as $\sigma^{2}/K$. Initially, we impose a Gaussian prior $p(\theta)$ with zero
mean and covariance matrix $\Sigma_{0}$ on weight vector $\theta$. According
to Bayes’ theorem, the posterior distribution can be calculated analytically
as

$P(\theta|Y)=\mathcal{N}(\bar{\theta},\bar{\Sigma})$

with $\bar{\theta}$ the maximum a posterior estimation (MAP) as

$\bar{\theta}=(\Phi^{T}\Phi+\frac{\sigma^{2}}{K}\Sigma_{0}^{-1})^{-1}\Phi^{T}Y$
(30)

and $\bar{\Sigma}$ as the posterior covariance matrix as

$\bar{\Sigma}=\frac{\sigma^{2}}{K}(\Phi^{T}\Phi+\frac{\sigma^{2}}{K}\Sigma_{0}^{-1})^{-1}$

As a result, given training data $X$, $Y_{f}$ and $Y_{g}$ as in Equation 28,
we can approximate the drift terms $\hat{f}(x)$ and $\hat{g}(x)$. By properly
choosing the base functions $\Phi_{\hat{f}}$ and $\Phi_{\hat{g}}$, we can
construct models for both $\hat{f}(x)$ and $\hat{g}(x)$ according to Equation
29 as

$\begin{split}\hat{f}(x)&=\Phi_{f}\theta_{f}+\varepsilon_{f}\\\
\hat{g}(x)&=\Phi_{g}\theta_{g}+\varepsilon_{g}\\\ \end{split}$

and calculate weights $\bar{\theta}_{f}$ and $\bar{\theta}_{g}$ based on
Equation 30.

###  3.3 Data collection for diffusion term

After we estimate the drift terms, we can use the estimated $\hat{f}(x)$ and
$\hat{g}(x)$ to identify the diffusion term. Similar as above, we also
randomly sample an initial state and let the system evolve under some control
$u_{i}$ with $i\in\mathcal{I}$ at each time step given the refined partition
$\\{t_{0},t_{1},...,t_{N}\\}$. Specifically, at each time step, we calculate
the $\xi_{i}$ using estimated drift term as

$\xi_{i}=x_{i+1}-x_{i}-\hat{f}(x_{i})\Delta t-\hat{g}(x_{i})u_{i}\Delta t.$

Then we will get the dataset $\mathcal{D}_{\xi}$ containing
$[\xi_{0},\dots,\xi_{N-1}]^{T}$ to estimate $\sigma$.

###  3.4 Identification of diffusion term

The likelihood function given dataset $\mathcal{D}$ is

$L(\mathcal{D}|\sigma)=\Pi_{i=0}^{N-1}P(\xi_{i}|\sigma),$

where

$P(\xi_{i}|\sigma)=\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\frac{1}{2}\frac{\xi_{i}^{2}}{\sigma^{2}}}$

Given a prior distribution $P(\sigma)$, the posterior distribution of $\sigma$
is

$P(\sigma|\mathcal{D})\varpropto L(\mathcal{D}|\sigma)\cdot P(\sigma).$ (31)

Then we use the MAP method to estimate $\sigma$ by maximizing the log-
posterior distribution. According to Eq (31), we can get

$\log(P(\sigma|\mathcal{D}))=\log{L(\sigma|\mathcal{D})+\log(P(\sigma))}$

By using the inverse-gamma distribution

$P(\sigma)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}\frac{1}{\sigma^{{}^{\alpha+1}}}e^{-\frac{\beta}{\sigma}}$
(32)

as the prior distribution, we can get

$\log{P(\sigma)|\mathcal{D}}=-(N+\alpha)\log{\sigma}-\frac{1}{2\sigma^{2}}\sum_{i=0}^{N-1}\xi_{i}^{2}-\frac{\beta}{\sigma}+C$

with

$C=\alpha\log{\beta}-\log{\Gamma(\alpha)}-(N-1)\log{\sqrt{2\pi}}$

is a constant [32]. The best estimation $\hat{\sigma}$ is obtained by taking
the derivative of the log-posterior with respect to $\sigma$ as

$\frac{d\log{P(\hat{\sigma}|\mathcal{D})}}{d\hat{\sigma}}=0.$ (33)

###  3.5 QP-based control framework for learned dynamics with SCBFs

Based on the above identification process, we can make prediction using the
learned model to approximate the dynamics of an unknown SDE. For drift terms,
we can make prediction as $\hat{f}(x)=\Phi_{f}\bar{\theta}_{f}$ and
$\hat{g}(x)=\Phi_{g}\bar{\theta}_{g}$ and for diffusion term, we estimate
noise level using $\hat{\sigma}$ as in Equation 33. As a result, we can
control the unknown SDE using a QP-based control framework with SCBF as

$\begin{split}u(x)&=\operatorname*{arg\,min}_{u\in\mathbb{R}}\frac{1}{2}||u||^{2},\\\
&\text{s.t.}\quad\hat{\mathcal{A}}B(x)\geq 0,\end{split}$ (34)

where

$\hat{\mathcal{A}}B(x)=\frac{\partial B}{\partial
x}(\hat{f}(x)+\hat{g}(x)u+\frac{1}{2}\sum\limits_{i,j}\left(\hat{\sigma}\hat{\sigma}^{T}\right)_{i,j}(x)\frac{\partial^{2}B}{\partial
x_{i}\partial x_{j}}.$

##  4 Simulation results

###  4.1 Example 1

In the first example, we test our result using a nonlinear system given by the
following stochastic differential equation as in [27]:

$d\begin{bmatrix}\dot{x}_{1}\\\
\dot{x}_{2}\end{bmatrix}=\begin{bmatrix}-0.6x_{1}-x_{2}\\\
x_{1}^{3}\end{bmatrix}dt+\begin{bmatrix}0\\\
x_{2}\end{bmatrix}udt+\begin{bmatrix}\sigma_{1}&0\\\
0&\sigma_{2}\end{bmatrix}dW.$

The safe region is defined as

$h(x)=-x_{2}^{2}-x_{1}^{2}+1>0.$

The generator of $h$ is calculated as

$\mathcal{A}(h)=1.2x_{1}^{2}+x_{1}x_{2}-2x_{1}^{3}x_{2}-2x_{2}^{2}u-(\sigma_{2}^{2}+\sigma_{2}^{2}),$

and we select

$\begin{split}b_{0}(x)&=h(x),\\\ b_{1}(x)&=\mathcal{A}(x).\end{split}$

In order to identify $f(x)$ in the drift term, we first randomly sample 10
initial points and at each initial point, we set $u=0$ and simulate the system
with $\Delta t=0.01$. For any state $x_{t}=\xi$, we apply the control to the
system at state $\xi$ for $K$ times to calculate the expectation of $f(x_{t})$
using the central limit theorem as

$f(\xi)+\varepsilon\approx\sum_{j=0}^{K}\frac{[x_{t_{K}+1}-x_{t}|x_{t}=\xi]}{K\Delta
t}.$ (35)

The we will get the training data for $f(x)$ as

$\begin{split}\mathbb{X}&=[\xi_{1},\xi_{2},...,.\xi_{N}]^{T},\\\
\mathbb{Y}&=[f(\xi_{1}),f(\xi_{2}),...,f(\xi_{N})].\end{split}$

and estimate $f(x)$ using linear regression. After $f(x)$ is estimated, we use
the estimated function $\hat{f}(x)$ to identify $g(x)$ in the drift term. We
set control input as $u=0.1$ since the system is control affined and $g(x)$ is
independent of $u$. We collect the training data similarly as in estimating
$f(x)$. Through out this example of the simulation, we use polynomial base
functions as

$\Phi=[1,x_{1},x_{2},x_{1}^{2},x_{2}^{2},x_{1}x_{2},x_{1}^{3},x_{2}^{3}]$

for both $f(x)$ and $g(x)$. The result of estimation for $f(x)$ is shown in
Fig 1. We also calculate the mean square error (MSE) for $f_{1}(x)$ and
$f_{2}(x)$ over 100 randomly sample points and show the result in Table 1.

Table 1: MSE for estimation of $f_{1}(x)$ and $f_{2}(x)$ with $K=10,30,100$ over 100 randomly sample points. | $f_{1}(x)$ | $f_{2}(x)$
---|---|---
K=10 | $e^{-2}$ | $5e^{-2}$
K=30 | $2e^{-3}$ | $8e^{-3}$
K=100 | $6e^{-4}$ | $6e^{-4}$
Figure 1: Estimation of $f_{2}(x)=x_{1}^{3}$ using $K=10$, $K=30$ and $K=100$.
The estimated result is compared with the true value of the function using 30
randomly sampled points. The diffusion of the system for all the cases is
$[0.2,0.2]$.

After we estimate the drift term, we collect data using the estimated
$\hat{f}(x)$ and $\hat{g}(x)$ for 100 initial points and 300 simulation steps
for each initial point. At each time step, we calculate the noise data using

$\phi_{t}=x_{t+1}-x_{t}-\hat{f}(x_{t})\cdot dt-\hat{g}(x_{t})\cdot u\cdot dt$

and add it into the data set $\mathcal{D}$. We use inverse gamma distribution
as in Equation 32 with $\alpha=1$ and $\beta=1$.

We use log-likelihood function as

$\log
L(\sigma|\mathcal{D})=\sum_{i=1}^{n}\left[-\frac{1}{2}\log(2\pi\sigma^{2})-\frac{\phi_{i}^{2}}{2\sigma^{2}}\right],$

where $n$ is the number of data in $\mathcal{D}$. We estimate a $\hat{\sigma}$
such that the log-posterior distribution

$\log P(\hat{\sigma}|\mathcal{D})=\log P(\sigma)+\log
L(\hat{\sigma}|\mathcal{D})$

is maximized. The distribution over 10000 random sample from posterior
distribution of $\sigma_{1}$ and $\sigma_{2}$ is shown in Fig 2.

Figure 2: Posterior distribution of $\sigma_{1}$ and $\sigma_{2}$ over 10000
random samples for Example 2.

In order to analyze the performance of safety control w.r.t the estimated
system, we count the number of safe trajectories over 1000 simulations under
different initial points between the SCBF, DDSCBF [27] and Bayesian SCBF. For
SCBF, we assume that we have knowledge of true system and for DDSCBF, we
assume that we know the drift of the system and only the diffusion is unknown
to us. The statistical results at initial state $[-0.1,0.7]^{T}$ and
$[-0.1,0.8]^{T}$ are shown in Table 2. The analytical result in the table is
the worst-case safe probability calculated based on [19, Theorem III.8]. Since
the relative degree if of $r=1$, the worst-case probability is calculated as
$P=\frac{b_{0}(x_{0})}{c_{0}}$ where $c_{0}=\sup_{x\in\mathcal{C}}h(x)$ with
$\mathcal{C}$ as the safe set. It is easy to find that $c_{0}=1$ so the worst-
case probability is calculated as $P=0.5$ and $P=0.35$ for initial state
$[-0.1,0.7]^{T}$ and $[-0.1,0.8]^{T}$, respectively.

Table 2: Safety ratio over 1000 simulation trajectories at initial state $[-0.1,0.7]^{T}$ and $[-0.1,0.8]^{T}$ with diffusion $\sigma_{1},\sigma_{2}=0.2$. | Safety ratio |
---|---|---
| $(-0.1,0.7)$ | $(-0.1,0.8)$ |
SCBF[19] | 91% | 55% |
DDSCBF[27] | 88% | 45% |
Bayesian SCBF | 90% | 43% |
Analytical[19] | 50% | 35% |

###  4.2 Example 2

In the second example, we use an adaptive cruise control (ACC) model as in
[33] to validate the estimation of the system and safety control. The model is

$d\begin{bmatrix}v\\\ z\\\ \end{bmatrix}=\begin{bmatrix}-F_{r}(x)/M\\\
v_{f}-v\\\ \end{bmatrix}dt+\begin{bmatrix}1/M\\\ 0\\\
\end{bmatrix}udt+\begin{bmatrix}\sigma_{1}&0\\\ 0&\sigma_{2}\\\
\end{bmatrix}dW,$

where $x=[v,z]$ is the state of the system representing the speed of the
behind vehicle and the distance between two vehicles, respectively. $v_{f}$ is
the speed of the front vehicle. The aerodynamic drag is
$F_{r}(x)=f_{0}+f_{1}x_{1}+f_{2}x_{1}^{2}$ with $f_{0}=0.1$, $f_{1}=5$,
$f_{2}=0.25$ and the mass of the vehicle is $M=1650$. We require the behind
vehicle to reach a desired speed $v_{d}=22$ while keeping a minimum distance
with $D=10$ from the front vehicle. Since the second-order generator of $h(x)$
is

$\mathcal{A}\mathcal{A}h(x)=L^{2}_{f}h(x)+L_{g}L_{f}h(x)u+\Sigma^{T}\frac{\partial^{4}h(x)}{\partial
x^{4}}\Sigma$

with

$\Sigma=\begin{bmatrix}\sigma_{1}&0\\\ 0&\sigma_{2}\\\ \end{bmatrix},$

we select $h(x)=(z-D)^{5}$ with a non-zero fourth-order derivative. As a
result, select

$\begin{split}b_{0}(x)&=h(x),\\\
b_{1}(x)&=\mathcal{A}h(x)=5(z-D)^{4}(v_{f}-v),\\\
b_{2}(x)&=\frac{5(z-D)^{4}}{M}\cdot F+20(z-D)^{3}\cdot(v_{f}-v)^{2}\\\
&-\frac{(z-D)^{4}}{M}\cdot u+120(z-D)(\sigma_{1}^{2}+\sigma_{2}^{2}).\\\
\end{split}$

We use the same base functions, prior distribution and likelihood function as
in the first example for system estimation. The result of estimation for
$f(x)$ is displayed as in Fig 3 and the distribution over 10000 random samples
from the posterior distribution of $\sigma_{1}$ and $\sigma_{2}$ is shown in
Fig 4.

(a)

(b)

Figure 3: Estimation of $f(x)$ for Example 2. The result is validated using
100 randomly sampled states. (a): Estimation of $f_{1}(x)$. (b): Estimation of
$f_{2}(x)$. Figure 4: Posterior distribution of $\sigma_{1}$ and $\sigma_{2}$
over 10000 random samples for Example 2. The true values are $\sigma_{1}=0.5$
and $\sigma_{2}=0.5$.

To control the vehicle to reach the desired speed, we use control Lyapunov
function as in [19]. Similar as in example 1, we calculate the safety ratio
over 1000 simulations between SCBF, DDSCBF and Bayesian SCBF. Note that the
SCBF is a non-convex function in this example, we can not find the supreme
value of $c_{0}$ as in the first example. So we only compare simulated safety-
probability and the result is shown in Table 3.

Table 3: Safety ratio over 1000 simulation trajectories with diffusion $\sigma_{1},\sigma_{2}=0.5$. The initial state is $[v,z]^{T}=[10,15]^{T}$. | Safety ratio
---|---
SCBF[19] | 83%
DDSCBF[27] | 65%
Bayesian SCBF | 78%

Also, as stated in [27], the bottleneck of the DDSCBF is that when the
relative degree is higher than 1, we have to recursively apply this method in
each order to calculate higher order generators. As in this example, we have
relative degree $r=2$, we need to learn the first-order generator, and then
use the estimated first-order generator to sample data for estimating the
second order generator. In the process of collecting training data for second
order generator using DDSCBF, we have to sample $N$ points and at each point,
sample $K$ simulations and apply CLT to calculate training data. However, in
each simulation, we need to get the value of first-order generator through
network, which is very time consuming. As a result, we compare the running
time between DDSCBF and Bayesian estimated SCBF. The runtime of both methods
are displayed in Table 4. In the simulation of DDSCBF, we select $N=300$ and
$K=10000$.

Table 4: Runtime of learning process | Runtime
---|---
DDSCBF[27] | over 5 hours
Bayesian SCBF | 15s

##  5 conclusion

This paper is the first one to address safety-critical control problem for
fully unknown stochastic differential equations. In this paper, we use
Bayesian inference to estimate both the drift and diffusion terms of the
system and use SCBF to guarantee safety for the learned system. We validate
safety ratio with statistical results using nonlinear models and compare with
worse-case analytical safety probabilities. Future work would be focused on
safety-critical control with limited observations with sparse learning
inspired by recently work in SDE identifications.

## References

  * [1] J. Garcıa and F. Fernández, “A comprehensive survey on safe reinforcement learning,” Journal of Machine Learning Research, vol. 16, no. 1, pp. 1437–1480, 2015.
  * [2] L. Lamport, “Proving the correctness of multiprocess programs,” IEEE Transactions on Software Engineering, no. 2, pp. 125–143, 1977.
  * [3] S. Ratschan and Z. She, “Safety verification of hybrid systems by constraint propagation-based abstraction refinement,” ACM Transactions on Embedded Computing Systems (TECS), vol. 6, no. 1, pp. 8–es, 2007.
  * [4] A. Girard, C. Le Guernic, and O. Maler, “Efficient computation of reachable sets of linear time-invariant systems with inputs,” in Proc. of HSCC, pp. 257–271, Springer, 2006.
  * [5] A. D. Ames, S. Coogan, M. Egerstedt, G. Notomista, K. Sreenath, and P. Tabuada, “Control barrier functions: Theory and applications,” in Proc. of ECC, pp. 3420–3431, IEEE, 2019.
  * [6] A. D. Ames, J. W. Grizzle, and P. Tabuada, “Control barrier function based quadratic programs with application to adaptive cruise control,” in Proc. of CDC, pp. 6271–6278, IEEE, 2014.
  * [7] M. Rauscher, M. Kimmel, and S. Hirche, “Constrained robot control using control barrier functions,” in Proc. of IROS, pp. 279–285, IEEE, 2016\.
  * [8] A. D. Ames, X. Xu, J. W. Grizzle, and P. Tabuada, “Control barrier function based quadratic programs for safety critical systems,” IEEE Transactions on Automatic Control, vol. 62, no. 8, pp. 3861–3876, 2016.
  * [9] Y. Chen, H. Peng, and J. Grizzle, “Obstacle avoidance for low-speed autonomous vehicles with barrier function,” IEEE Transactions on Control Systems Technology, vol. 26, no. 1, pp. 194–206, 2017.
  * [10] S.-C. Hsu, X. Xu, and A. D. Ames, “Control barrier function based quadratic programs with application to bipedal robotic walking,” in Proc. of ACC, pp. 4542–4548, IEEE, 2015.
  * [11] Q. Nguyen and K. Sreenath, “Exponential control barrier functions for enforcing high relative-degree safety-critical constraints,” in Proc. of ACC, pp. 322–328, IEEE, 2016.
  * [12] W. Xiao and C. Belta, “Control barrier functions for systems with high relative degree,” in Proc. of CDC, pp. 474–479, IEEE, 2019.
  * [13] A. Taylor, A. Singletary, Y. Yue, and A. Ames, “Learning for safety-critical control with control barrier functions,” in Proc. of L4DC, pp. 708–717, PMLR, 2020.
  * [14] C. Wang, Y. Meng, Y. Li, S. L. Smith, and J. Liu, “Learning control barrier functions with high relative degree for safety-critical control,” in 2021 European Control Conference (ECC), pp. 1459–1464, IEEE, 2021.
  * [15] H. J. Kushner, Stochastic Stability and Control.  Academic Press, 1967.
  * [16] M. Sarkar, D. Ghose, and E. A. Theodorou, “High-relative degree stochastic control Lyapunov and barrier functions,” arXiv preprint arXiv:2004.03856, 2020.
  * [17] A. Clark, “Control barrier functions for complete and incomplete information stochastic systems,” in Proc. of ACC, pp. 2928–2935, IEEE, 2019.
  * [18] C. Santoyo, M. Dutreix, and S. Coogan, “A barrier function approach to finite-time stochastic system verification and control,” Automatica, vol. 125, p. 109439, 2021.
  * [19] C. Wang, Y. Meng, S. L. Smith, and J. Liu, “Safety-critical control of stochastic systems using stochastic control barrier functions,” in 2021 IEEE Conference on Decision and Control (CDC), pp. 5924–5931, 2021.
  * [20] R. Cheng, M. J. Khojasteh, A. D. Ames, and J. W. Burdick, “Safe multi-agent interaction through robust control barrier functions with learned uncertainties,” in 2020 59th IEEE Conference on Decision and Control (CDC), pp. 777–783, IEEE, 2020.
  * [21] P. Jagtap, G. J. Pappas, and M. Zamani, “Control barrier functions for unknown nonlinear systems using gaussian processes,” in 2020 59th IEEE Conference on Decision and Control (CDC), pp. 3699–3704, IEEE, 2020.
  * [22] V. Dhiman, M. J. Khojasteh, M. Franceschetti, and N. Atanasov, “Control barriers in bayesian learning of system dynamics,” IEEE Transactions on Automatic Control, 2021.
  * [23] L. Brunke, S. Zhou, and A. P. Schoellig, “Barrier bayesian linear regression: Online learning of control barrier conditions for safety-critical control of uncertain systems,” in Learning for Dynamics and Control Conference, pp. 881–892, PMLR, 2022.
  * [24] Y. Chen, C. Shang, X. Huang, and X. Yin, “Data-driven safe controller synthesis for deterministic systems: A posteriori method with validation tests,” arXiv preprint arXiv:2304.00729, 2023.
  * [25] A. Salamati, A. Lavaei, S. Soudjani, and M. Zamani, “Data-driven verification and synthesis of stochastic systems through barrier certificates,” arXiv preprint arXiv:2111.10330, 2021.
  * [26] R. Wajid, A. U. Awan, and M. Zamani, “Formal synthesis of safety controllers for unknown stochastic control systems using gaussian process learning,” in Learning for Dynamics and Control Conference, pp. 624–636, PMLR, 2022.
  * [27] C. Wang, Y. Meng, S. L. Smith, and J. Liu, “Data-driven learning of safety-critical control with stochastic control barrier functions,” in 2022 IEEE 61st Conference on Decision and Control (CDC), pp. 5309–5315, IEEE, 2022.
  * [28] C. Archambeau, D. Cornford, M. Opper, and J. Shawe-Taylor, “Gaussian process approximations of stochastic differential equations,” in Gaussian Processes in Practice, pp. 1–16, PMLR, 2007.
  * [29] S. Särkkä et al., Recursive Bayesian inference on stochastic differential equations.  Helsinki University of Technology, 2006.
  * [30] I. Castillo, J. Schmidt-Hieber, and A. Van der Vaart, “Bayesian linear regression with sparse priors,” 2015.
  * [31] C. Wang, Y. Li, Y. Meng, S. L. Smith, and J. Liu, “Learning control barrier functions with high relative degree for safety-critical control,” in Proc. of ECC, to appear, 2021.
  * [32] G. E. Box and G. C. Tiao, Bayesian inference in statistical analysis.  John Wiley & Sons, 2011.
  * [33] W. Xiao, C. A. Belta, and C. G. Cassandras, “Sufficient conditions for feasibility of optimal control problems using control barrier functions,” Automatica, vol. 135, p. 109960, 2022.