# First-order PDES for Graph Neural Networks: Advection And Burgers Equation
Models ††thanks: Note: This work was performed for an undergraduate research
project from May to August 2023.

Yifan Qu
Department of Applied Math
University of Waterloo
Ontario, Canada
<EMAIL_ADDRESS>
&Oliver Krzysik
Department of Applied Math
University of Waterloo
Ontario, Canada
<EMAIL_ADDRESS>
Hans De Sterck
Department of Applied Math
University of Waterloo
Ontario, Canada
<EMAIL_ADDRESS>
&Omer Ege Kara
Department of Applied Math
University of Waterloo
Ontario, Canada
<EMAIL_ADDRESS>

###### Abstract

Graph Neural Networks (GNNs) have established themselves as the preferred
methodology in a multitude of domains, ranging from computer vision to
computational biology, especially in contexts where data inherently conform to
graph structures. While many existing methods have endeavored to model GNNs
using various techniques, a prevalent challenge they grapple with is the issue
of over-smoothing. This paper presents new Graph Neural Network models that
incorporate two first-order Partial Differential Equations (PDEs). These
models do not increase complexity but effectively mitigate the over-smoothing
problem. Our experimental findings highlight the capacity of our new PDE model
to achieve comparable results with higher-order PDE models and fix the over-
smoothing problem up to 64 layers. These results underscore the adaptability
and versatility of GNNs, indicating that unconventional approaches can yield
outcomes on par with established techniques.

_K_ eywords Graph Neural Network $\cdot$ Partial Differential Equation $\cdot$
Advection Equation $\cdot$ Burgers Equation

##  1 Introduction

In recent years, there has been a great demand for the classification and
analysis of complex data such as images [1, 2], citation networks, molecular
graphs [3, 4, 5], and social networks [6]. Graphs with nodes representing
entities and edges representing relationships express these network
relationships most effectively. Nevertheless, the networks in real life are
intricate, hence how to extract valuable information from these graphs is a
challenge.

One of the significant limitations of traditional Graph Neural Networks (GNNs)
methods is the over-smoothing problem [7, 8, 9]. It is caused by multiple
layers of graph convolutions applied in the GNN model: after the information
propagates through several layers the representations tend to become
indistinguishable [10, 11]. It results in the blending of node and edge
features and lessens expressive strength. Some existing methods have explored
strategies to alleviate or overcome this challenge [9, 12, 13].

##  2 Related Work

Among the existing methods for GNN, Graph Convolution Network (GCN) models are
used most frequently. GCN was first proposed by Bruna et al. [14], which
considers spectral convolutions on graphs. The polynomial of the Laplacian
operator acting as the spectral convolution operator is a key consideration in
the GCN methods outlined in [15, 16]. The general GCN model structure includes
four main components: activation function $\sigma(\cdot)$; propagation
operator $\mathbf{P}$; feature matrix at $l^{th}$ layer
$\mathbf{F}^{(l)}\in\mathbb{R}^{N\times D_{l}}$, and layer-specific trainable
weight matrix $\mathbf{W}^{(l)}\in\mathbb{R}^{D_{l}\times D_{l+1}}$[16].

$\mathbf{F}^{(l+1)}=\sigma\left(\mathbf{P}\mathbf{F}^{(l)}\mathbf{W}^{(l)}\right).$
(1)

Here,
$\mathbf{P}=\mathbf{\tilde{D}}^{-\frac{1}{2}}\mathbf{\tilde{A}}\mathbf{\tilde{D}}^{-\frac{1}{2}}$
with degree matrix $\mathbf{D}\in\mathbb{R}^{N\times N}$ and adjacency matrix
$\mathbf{A}\in\mathbb{R}^{N\times N}$, and
$\mathbf{\tilde{A}}=\mathbf{A}+\mathbf{I}_{N},\mathbf{\tilde{D}}=\mathbf{D}+\mathbf{I}_{N}$
are the adjacency and degree matrix with self-loops on the original graph.
$\mathbf{I}_{N}\in\mathbb{R}^{N\times N}$ is the identity matrix. The formula
presented above is based on the non-negative smoothing operator $\mathbf{P}$,
which results in increasingly indistinguishable features as the layer count
increases. To address this issue, Eliasof et al. proposed an innovative
solution by introducing a trainable propagation parameter $\omega$ to weigh
each layer [12]. However, despite the improvement in GNN results for various
tasks like node classification, these improved GNNs still face limitations in
their applicability. For instance, a model that performs exceptionally well in
node classification may yield poor results in shape correspondence problems
[7, 12]. This problem of restricted applicability arises partly due to the
lack of robust theoretical support for GCN-based graph neural networks.

In recent years, attention has turned towards Partial Differential Equation
(PDE)-based convolutional neural networks for graph neural networks, offering
promising avenues for improvement [17]. PDE models, commonly used in computer
vision and constructing neural networks, have been adapted to devise notable
models for graph neural networks [17, 18, 19]. The underlying concept is to
view the graph neural network as a discretization of a dynamic system governed
by differential equations, where each layer represents a progressive time step
from the previous one, and the step size can be treated as a trainable
parameter. By employing forward Euler discretization in time, it becomes
feasible to discretize the parabolic forward propagation effectively. Existing
methods have successfully adapted hyperbolic and diffusive second-order PDEs
[13] for use in GNNs. These methods utilize a trainable parameter $\alpha$ to
control the weighting between the two PDEs. As a result, they have achieved
significantly improved results in both node classification tasks and shape
correspondence tasks. The combination of these two PDEs allows the model to
automatically select the most suitable PDE for different scenarios, enhancing
the overall performance and versatility of the GNN. The main goal of the
current paper is to explore the relevance of first-order PDEs. Notably, a
recent study by Eliasof et al. [20] introduced an advection-based GNN model
that incorporates an advection equation, as well as diffusion and reaction
components, in their investigation of spatiotemporal node forecasting
problems, coinciding with the time of our research.

##  3 Method

Existing methods have almost exclusively focused on second-order PDEs within
their models, with first-order PDEs remaining unexplored until recently. Our
objective is to incorporate first-order PDEs into our framework and assess
their effectiveness. First-order PDEs, such as the linear advection equation
and the nonlinear Burgers equation, offer several advantages in the context of
GNNs. A first advantage is that first-order PDEs are often simpler to
implement and computationally more efficient compared to higher-order PDEs.
They involve only first derivatives with respect to space or time, which can
lead to more straightforward numerical schemes. Second, first-order PDEs often
involve conservation laws, which can help ensure that important quantities
(e.g., mass, energy, momentum) are conserved during simulations. This is
crucial for maintaining physical realism in simulations and predictions. In
this paper, we will integrate two first-order PDE models, namely the advection
and Burgers equations, into GNNs for examination and evaluation.

###  3.1 Advection Model

The advection equation, as cited in [21], stands as a quintessential example
of a first-order hyperbolic PDE frequently harnessed in fluid dynamics
modeling. The general formulation for the linear advection equation in two
spatial dimensions is expressed as Equation 2:

$\displaystyle u_{t}+au_{x}+bu_{y}$ $\displaystyle=0.$ (2)

Here, $u$ represents the quantity being transported, and $t$ represents time,
$u_{t}$ is the partial derivative of $u$ with respect to time, while $a$ and
$b$ are constants that affect the rate of change of $u$ with respect to the
horizontal and vertical direction respectively. The Advection equation models
how $u$ changes as time progresses. Within the framework of GCNs, we let this
advection equation undergo a transformation. It adapts to the dynamic context
of graph structures illustrating how information flows across the nodes of a
graph. The integration of the advection equation into GCNs may offer some
advantages.

First and foremost, the advection equation may emerge as an adept guardian of
spatial information, surpassing its diffusion-based counterparts like the heat
equation. In the realm of information propagation, it excels at preserving
sharp, localized features [22]. The GCN may benefit from this characteristic,
ensuring that pertinent spatial nuances remain intact throughout the
dissemination process.

Second, as a first-order PDE, the advection equation may help with the over-
smoothing predicament that is often encountered in higher-order PDEs. It
produces a controlled flow of information while safeguarding the unique
identities of individual nodes. Consequently, even in the depths of intricate
GCN architectures, the risk of information loss is curtailed.

Last but not least, the advection equation’s remarkable simplicity makes the
model easier to implement. The advection equation involves a simple first-
order derivative. This simplicity significantly reduces the complexity of
implementing the model within a GNN framework, making it accessible to a broad
range of researchers and practitioners.

The adoption of the advection equation in our model is further informed by
discretizing equation 2 on a 2D regular grid and transforming this to a graph
with n nodes and m edges.

$\displaystyle\frac{u_{ij}^{l+1}-u_{ij}^{l}}{\Delta
t}+\nabla\cdot(a,b)u^{l}=0$ (3) $\displaystyle
u^{l+1}-u^{l}=-h\sigma(G^{T}A^{T}u^{l}K_{l})$ (4)

Equation 4 emerges as a result of the forward Euler discretization technique,
offering practical applicability to the evolving dynamics of graph-based
information diffusion. Here, $u^{l}\in\mathbb{R}^{n\times d}$ denotes the node
features residing within layer $l$, $h$ is the hyperparameter reflecting the
discretization step size, and $K_{l}\in\mathbb{R}^{d\times d}$ embodies the
$1\times 1$ convolution kernel positioned at layer $l$. Notably, ReLU is used
as the activation function $\sigma$ in our experimental setup. Furthermore,
the averaging operator $A^{T}$, with dimension $m\times n$, executes a pivotal
role in the conversion of node features from their intrinsic node-centric
space to a broader edge-centric perspective by averaging [13]. Finally,
$G\in\mathbb{R}^{n\times m}$ is the graph incidence matrix, and left
multiplication by $G^{T}$ encodes the divergence operator on the graph.

This amalgamation of mathematical principles and adaptable techniques
characterizes the innovation driving our model, extending the horizons of
graph-based machine learning.

###  3.2 Burgers Equation

Similar to the advection model, the introduction of the nonlinear Burgers
equation into the framework of GNN presents an approach to use the inherent
characteristics of this first-order partial differential equation to augment
the performance of graph-based models. The Burgers equation models the
behaviors of fluid flow and the formation of shockwaves. Equation 5, the
Burgers equation, encapsulates these dynamics:

$\displaystyle u_{t}+(u^{2}/2)_{x}$ $\displaystyle=0.$ (5)

Here, $u$ represents the velocity of the fluid or the quantity being
transported. In the context of fluid dynamics, it represents the fluid’s
velocity in the $x$-direction. To integrate this equation into the GNN
paradigm, we employ a forward Euler discretization technique applied to
Equation 5, yielding Equation 6. This discretization process allows us to
adapt Burgers equation to the dynamic landscape of evolving graphs, where node
quantities evolve over time:

$\displaystyle\frac{u^{l+1}-u^{l}}{h}+\frac{1}{2}~{}G^{T}A^{T}(u^{l}K_{l}\odot
u^{l}K_{l})$ $\displaystyle=0.$ (6)

Here, $u^{l}$ denotes the node feature within layer $l$, and the
hyperparameter $h$ corresponds to the step size applied during the
discretization process. Furthermore, $K_{l}$ represents the $1\times 1$
convolution kernel at layer $l$, while $A^{T}$ represents a $m\times n$
averaging operator that plays a pivotal role in transforming the node features
from the individual nodes to the edge space.

The use of the element-wise multiplication operator ($\odot$) signifies the
element-wise multiplication of the node feature matrices, reflecting the
interplay of these components within the context of our adapted Burgers
equation. On the other hand, this elementwise square function also makes the
model nonlinear. Hence a nonlinear activation function is no longer needed in
this model.

This adaptation uses a model from fluid dynamics in GNNs, endowing our model
with a different capability to capture and predict dynamic changes in complex
networks.

###  3.3 Mixing models

We have integrated elements from diffusion and wave equations into our
advection model. Drawing inspiration from Eliasof et al. [13], we have
introduced a trainable parameter, denoted as $\alpha$. This parameter holds
the key to dynamically selecting the most suitable model, whether it is based
on advection, diffusion, or wave dynamics, depending on the specific
characteristics of the problem at hand.

The nonlinear diffusion equation 7, and the wave equation 8 are two building
blocks for this adaptive framework:

$\displaystyle u^{l+1}=u^{l}-hG^{T}\sigma(Gu^{l}K_{l})K_{l}^{T},$ (7)
$\displaystyle u^{l+1}=2u^{l}-u^{l-1}-h^{2}G^{T}\sigma(Gu^{l}K_{l})K_{l}^{T}.$
(8)

The innovation lies in our ability to blend these equations with our linear
advection equation, leading to the two mixed dynamics expressed as Equations 9
and 10:

$\displaystyle
u^{l+1}-u^{l}=-(1-\alpha)h^{2}G^{T}\sigma(D_{D}Gu^{l}K_{l})K_{l}^{T}-h\alpha
G^{T}\sigma(D_{W}A^{T}u^{l}K_{l}),$ (9)
$\displaystyle(1-\alpha)(u^{l+1}-2u^{l}+u^{l-1})+\alpha(u^{l+1}-u^{l})=-(1-\alpha)h^{2}G^{T}\sigma(D_{D}Gu^{l}K_{l})K_{l}^{T}-h\alpha
G^{T}\sigma(D_{W}A^{T}u^{l}K_{l}).$ (10)

Here, $\alpha$ is a dynamic parameter, subject to training during the model’s
learning process, with a range between 0 and 1. This parameter orchestrates
the fusion of advection, diffusion, and wave dynamics, ensuring that the model
adaptively selects the most suitable mechanism for the specific problem under
consideration.

Moreover, we have also introduced trainable diagonal weighting parameters,
$D$, somewhat similar to the methodology established in [12]. This edge-wise
parameter plays a pivotal role in modulating the propagation and diffusion
operators, further enhancing the model’s adaptability and its ability to
effectively capture the underlying dynamics of the data. The incorporation of
the parameter $D$ in the second-order operators introduces anisotropic
diffusion, and in the first-order operator, it models the wave speed along the
edge. Anisotropic diffusion and varying edge-directed wave speed, in this
context, allow the model to capture and leverage directional dependencies
within the graph structure, enabling more accurate and contextually relevant
diffusion and propagation, ultimately enhancing performance in tasks where
understanding the diffusion and flow of information along specific edges or
connections is pivotal for achieving superior results in graph neural
networks.

In our experimentation and analysis, we anticipate that the diffusion mixing,
owing to its anisotropic diffusion properties, will excel in scenarios such as
node classification tasks, where preserving the local spatial information is
crucial. By adopting this approach, we aim to transcend the limitations of
traditional fixed models and pave the way for dynamic graph analysis that
tailors its behavior to the idiosyncrasies of each problem it encounters. We
denote the advection and diffusion model in equation 9 as AD and the advection
wave mixing model in equation 10 as AW.

##  4 Experiments

In this section, we will delve into the experimental phase, where the
capabilities of our proposed models across three distinct problems are
assessed: semi-supervised node classification, fully-supervised node
classification inspired by the work of London et al. [23], and the challenging
3D shape correspondence task, as motivated by Wu et al. [24], see also [13].

For semi-supervised node classification, the conventional approach has long
favored diffusion models as the primary mechanism, aiming to preserve
localized features and spatial information, which are often of paramount
importance in node-classification problems. Consequently, the fusion of
diffusion dynamics with advection is posited in our model of equation 9, since
it is expected to outperform other alternatives.

Conversely, in the context of the dense shape correspondence problem, the
prevailing wisdom dictates a departure from traditional smoothing processes
[13]. Unlike node classification, this problem benefits from the retention of
raw, unaltered data, since the intricate details are crucial. Therefore, the
wave equation is anticipated, when integrated with the advection equation, to
exhibit superior performance. This conjecture arises from the intrinsic nature
of the wave equation, which is adept at propagating, high-fidelity
information, harmonizing with the demands of the dense shape correspondence
task.

For all experiments, we followed the architecture presented in [13]. For semi
and fully-supervised node classification, we used the architecture in Table 4.
For dense shape correspondence, we used the architecture in Table 5.

Dataset | Cora | Cite. | Pubm. | Cham.
---|---|---|---|---
Classes | 7 | 6 | 3 | 5
Node | 2,708 | 3,327 | 19,717 | 2,277
Edge | 5,429 | 4,732 | 44,338 | 36,101
Features | 1,433 | 3,703 | 500 | 2,325
Table 1:  Statistics of datasets used in semi-and fully supervised node-
classification experiments.

###  4.1 Semi-supervised node-classification

For the semi-supervised node classification task, we evaluated the performance
of our model across three benchmark datasets: Cora, CiteSeer, and PubMed [23].
To maintain consistency with established practices [25], we adopted the
standard training/validation/testing split for all three datasets, featuring
20 nodes per class for training, 500 nodes for validation, and 1,000 nodes for
testing. The statistics are presented in table 1. The main purpose of the
validation set is to provide an unbiased evaluation of a model fit while
tuning hyperparameters and preventing overfitting. The best parameters are
presented in table 2. We executed our experiments following the training
protocol outlined in [7].

Our findings, as shown in Table 6, illuminate several key insights. Most
notably, we observed that our advection and combined models effectively
mitigated the over-smoothing problem, as evidenced by the sustained accuracy
with increasing layers. This underscores the efficacy of the approach in
preserving crucial local information.

Moreover, in line with our expectations, the advection and diffusion mixing
model outperformed its counterpart, the advection and wave mixing model.
Although our models fall short of achieving the same level of accuracy as
previous methods, we posit that this disparity can be attributed to the common
understanding that diffusion models are more optimal approaches for node
classification tasks and advection processes are not beneficial in this
context.

Dataset | Loss | LR | WD | Channels | Dropout | h
---|---|---|---|---|---|---
Cora | Relu | $4.6\times 10^{-5}$ | $1.2\times 10^{-4}$ | 64 | 0.5 | 0.6
Cite | Relu | $1.0\times 10^{-5}$ | $8.1\times 10^{-3}$ | 256 | 0.7 | 0.3
Pubm | Relu | $2.4\times 10^{-5}$ | $1.2\times 10^{-4}$ | 256 | 0.6 | 0.7
Table 2:  Semi-Supervised classification hyper-parameters

###  4.2 Fully-supervised node-classification

In accordance with the methodology established by Pei et al. in their work on
Geometric Graph Convolutional Networks ([26]), we conducted our experiments on
four distinct datasets, namely Cora, CiteSeer, PubMed, and Chameleon. To
maintain consistency, we adhered to the identical data partitioning scheme,
allocating 60% of the data for training, 20% for validation, and the remaining
20% for testing. Furthermore, to ensure the robustness of our findings, we
adopted the practice of averaging performance metrics over 10 random data
splits, as previously recommended by Pei et al. ([26]).

To optimize our model’s hyperparameters across varying depths, we conducted an
exhaustive grid search, the best parameters are presented in table 3. The
results of our models are summarized in Table 7. Our analysis indicates that
our approach yields comparable performance to existing methods, suggesting the
efficacy of our methodology in addressing the specific task at hand.

Dataset | Loss | LR | WD | Channels | Dropout | h
---|---|---|---|---|---|---
Cora | Relu | $2.3\times 10^{-5}$ | $1\times 10^{-4}$ | 64 | 0.5 | 0.2
Cite | Relu | $2.1\times 10^{-4}$ | $1.1\times 10^{-4}$ | 64 | 0.6 | 0.3
Pubm | Relu | $4.3\times 10^{-5}$ | $2.6\times 10^{-4}$ | 64 | 0.5 | 0.4
Cham | Relu | $8\times 10^{-4}$ | $9.2\times 10^{-5}$ | 64 | 0.6 | 0.5
Table 3:  Fully-Supervised classification hyper-parameters Input Size | Layer | Output Size
---|---|---
$n\times c_{in}$ |  $1\times 1$ droupout | $n\times c_{in}$
$n\times c_{in}$ |  $1\times 1$ convolution | $n\times c$
$n\times c$ | Relu | $n\times c$
$n\times c$ | Advection/Burgers Block | $n\times c$
$n\times c$ |  $1\times 1$ droupout | $n\times c$
$n\times c$ |  $1\times 1$ convolution | $n\times c_{out}$
Table 4: Architecture used for semi-and fully-supervised node classification

###  4.3 Dense Shape Correspondence

Within the realm of hyperbolic dynamics, the pursuit of identifying dense
correspondences between shapes represents a well-established avenue of
exploration. This pursuit is primarily driven by our interest in modeling
localized motion dynamics. In essence, the task of unearthing these
correspondences among shapes closely mirrors the endeavor to understand how
one shape can be systematically mapped onto another through a transformation
process.

To facilitate this investigation, we make use of the FAUST dataset, as
initially introduced by Bogo et al. [27]. This dataset encompasses ten
distinct human shapes, each scanned in ten different poses, culminating in a
total of 6,890 nodes per individual per pose In accordance with the data
partitioning methodology prescribed by Monti et al. [1], we designate the
initial 80 subjects for the purpose of training, reserving the remaining 20
subjects for the testing phase. In our pursuit, we have adopted the training
and testing methods outlined in the work of Eliasof et al. [13].

The results of our experiment are presented in Table 8. We can notice that the
wave mixing model consistently outperforms the diffusion mixing model,
aligning with our expectations. Although our findings align with those of
other researchers, indicating the efficacy of our approach for this particular
task, it is important to acknowledge the challenge posed by surpassing
existing methods, which have achieved an accuracy level of 99.9%.

Input Size | Layer | Output Size
---|---|---
$n\times 4$ |  $1\times 1$ convolution | $n\times c$
$n\times c$ | Relu | $n\times c$
$n\times c$ | Advection/Burgers Block | $n\times c$
$n\times c$ |  $1\times 1$ convolution | $n\times c$
$n\times c$ | Relu | $n\times c$
$n\times c$ |  $1\times 1$ convolution | $n\times 512$
$n\times 512$ | ELU | $n\times 512$
$n\times 512$ | Fully-connected | $n\times n$
Table 5: Architecture used for dense shape correspondence.

##  5 Conclusion

In this study, we have introduced a novel approach by formulating a first-
order PDE model for GNNs. Our motivation was to explore alternative avenues
for enhancing GNNs’ interpretability and performance and reducing their
complexity. We rigorously tested our proposed model against existing methods,
seeking to uncover any potential advantages.

Surprisingly, our results demonstrated an intriguing similarity in performance
with second-order PDEs. This suggests that our first-order PDE model can
indeed rival established methods in terms of predictive accuracy and
effectiveness.

This study underscores the versatility and adaptability of GNNs, suggesting
that even unconventional approaches, such as our first-order PDE models, can
achieve comparable results.

In summary, our research introduces new PDE-based models for Graph Neural
Networks. We recognize that there is room for further refinement and
optimization to enhance performance, an avenue we intend to explore in future
research endeavors. Our method shows that new models can match established
techniques. This study encourages further exploration of unconventional
approaches to advance GNNs and enhance our understanding of their
capabilities.

Table 6:  Semi-supervised node classification accuracy (%). The best results in each method are bolded. |  | Layers
---|---|---
Dataset | Method | 2 | 4 | 8 | 16 | 32 | 64
Cora | GCN | $\mathbf{81.1}$ | 80.4 | 69.5 | 64.9 | 60.3 | 28.7
GCN(Drop) | $\mathbf{82.8}$ | 82.0 | 75.8 | 75.7 | 62.5 | 49.5
GCNII | 82.2 | 82.6 | 84.2 | 84.6 | 85.4 | $\mathbf{85.5}$
GCNII* | 80.2 | 82.3 | 82.8 | 83.5 | 84.9 | $\mathbf{85.3}$
PDE-GCN(Diffusion) | 82.0 | 83.6 | 84.0 | 84.2 | 84.3 | $\mathbf{84.3}$
Advection | 70.10 | 71.1 | 72.1 | 71.3 | 72.3 | $\mathbf{73.1}$
Burgers | 71.10 | 70.1 | 72.3 | 70.0 | 72.4 | $\mathbf{73.1}$
Mix(AD) | 77.2 | 78.1 | 77.1 | 78.3 | 78.6 | $\mathbf{80.1}$
Mix(AW) | $\mathbf{76.3}$ | 75.3 | 76.2 | 74.0 | 74.3 | 72.9
Citeseer | GCN | $\mathbf{70.8}$ | 67.6 | 30.2 | 18.3 | 25.0 | 20.0
GCN (Drop) | $\mathbf{72.3}$ | 70.6 | 61.4 | 57.2 | 41.6 | 34.4
GCNII | 68.2 | 68.8 | 70.6 | 72.9 | 73.4 | $\mathbf{73.4}$
GCNII* | 66.1 | 66.7 | 70.6 | 72.0 | $\mathbf{73.2}$ | 73.1
PDE-GCN(Diffusion) | 74.6 | 75.0 | 75.2 | 75.5 | $\mathbf{75.6}$ | 75.5
Advection | 74.3 | 72.4 | 75.1 | 75.0 | 72.3 | $\mathbf{75.6}$
Burgers | 72.3 | 71.1 | 74.3 | 72.4 | 71.4 | $\mathbf{74.3}$
Mix(AD) | 73.10 | 74.2 | 75.2 | 75.0 | 75.3 | $\mathbf{75.5}$
Mix(AW) | 72.2 | 72.5 | 72.8 | 73.1 | 73.3 | $\mathbf{74.1}$
Pubmed | GCN | $\mathbf{79.0}$ | 76.5 | 61.2 | 40.9 | 22.4 | 35.3
GCN (Drop) | $\mathbf{79.6}$ | 79.4 | 78.1 | 78.5 | 77.0 | 61.5
GCNII | 78.2 | 78.8 | 79.3 | $\mathbf{80.2}$ | 79.8 | 79.7
GCNII* | 77.7 | 78.2 | 78.8 | $\mathbf{80.3}$ | 79.8 | 80.1
PDE-GCN(Diffusion) | 79.3 | $\mathbf{80.6}$ | 80.1 | 80.4 | 80.2 | 80.3
Advection | 70.10 | 71.1 | 72.1 | 71.3 | 72.3 | $\mathbf{73.1}$
Burgers | 71.10 | 70.1 | 72.3 | 70.0 | $\mathbf{73.1}$ | 72.4
Mix(AD) | 77.2 | 77.8 | 78.3 | $\mathbf{79.9}$ | 78.73 | 79.1
Mix(AW) | 76.5 | 77.1 | 76.6 | 77.3 | $\mathbf{78.3}$ | 76.1
Method | Cora | Cite. | Pubm. | Cham.
---|---|---|---|---
GCN | 85.77 | 73.68 | 88.13 | 28.18
GAT | 86.37 | 74.32 | 87.62 | 42.93
Geom-GCN | 85.19 | 77.99 | 90.05 | 60.31
APPNP | 87.87 | 76.53 | 89.40 | 54.30
Incep (Drop) | 86.86(8) | 76.83(8) | 89.18(4) | 61.71(8)
GCNII | 88.49(64) | 77.08(64) | 89.57(64) | 60.61(8)
GCNII II | 88.01(64) | 77.13(64) | $\mathbf{90.30(64)}$ | 62.48(8)
PDE-GCN (Diffusion) | 88.51(16) | 78.36(64) | 89.6(64) | 64.12(8)
PDE-GCN (Wave) | 87.71(32) | 78.13(16) | 89.16(16) | 61.57(64)
PDE-GCN (Mixing) | $\mathbf{88.60(16)}$ | $\mathbf{78.48(32)}$ | 89.93(16) | $\mathbf{66.01(16)}$
Advection | 77.47(32) | 76.45(64) | 87.23(64) | 60.42(32)
Burgers | 78.60(32) | 75.76(64) | 88.68(64) | 64.34(32)
Mix(AD) | 79.92(64) | 77.65(64) | 89.98(32) | 65.89(16)
Mix (AW) | 77.54(32) | 77.89(64) | 87.16(64) | 60.29(32)
Table 7:  Fully-supervised node classification accuracy (%). The best results in each dataset are bolded. Method | FAUST
---|---
ACNN | 63.8
MoNet | 89.1
FMNet | 98.2
SplineCNN | 99.2
PDE-GCN(Wave) | 99.9
PDE-GCN(Diffusion) | 64.2
Advection | 72.2
Burgers | 80.4
Mix(AD) | 78.2
Mix(AW) | 98.9
Table 8:  Dense Shape Correspondence (%). The best results in others’ methods
and our methods are bolded.

## Acknowledgments

I wish to express my gratitude to the Applied Mathematics department at the
University of Waterloo for affording me, as an undergraduate student, the
opportunity to successfully complete this project under the guidance of our
supervisor. I am also deeply appreciative of the invaluable technical
assistance provided by Yanming Kang, the insightful contributions from our
colleagues, and the unwavering support bestowed upon us by our families and
friends throughout the course of this research endeavor.

## References

  * [1] Federico Monti, Davide Boscaini, Jonathan Masci, Emanuele Rodola, Jan Svoboda, and Michael M Bronstein.  Geometric deep learning on graphs and manifolds using mixture model cnns.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 5115–5124, 2017.
  * [2] Guangtao Wang, Rex Ying, Jing Huang, and Jure Leskovec.  Improving graph attention networks with large margin-based constraints.  arXiv preprint arXiv:1910.11945, 2019.
  * [3] Will Hamilton, Zhitao Ying, and Jure Leskovec.  Inductive representation learning on large graphs.  Advances in neural information processing systems, NIPS 30, 2017\.
  * [4] John Jumper, Richard Evans, Alexander Pritzel, Tim Green, Michael Figurnov, Olaf Ronneberger, Kathryn Tunyasuvunakool, Russ Bates, Augustin Žídek, Anna Potapenko, et al.  Applying and improving alphafold at casp14.  Wiley Online Library. Proteins: Structure, Function, and Bioinformatics, 89(12):1711–1721, 2021.
  * [5] David Duvenaud, Dougal Maclaurin, Jorge Aguilera-iparraguirre, G Rafael, Timothy Hirzel, and Ryan P Adams.  Convolutional Networks on Graphs for Learning Molecular Fingerprints arXiv : 1509 . 09292v2 [ cs . LG ] 3 Nov 2015.  pages 1–9.
  * [6] Muhan Zhang and Yixin Chen.  Link Prediction Based on Graph Neural Networks.  NIPS, 2018.
  * [7] Ming Chen, Zhewei Wei, Zengfeng Huang, Bolin Ding, and Yaliang Li.  Simple and deep graph convolutional networks.  In International conference on machine learning, pages 1725–1735. PMLR, 2020.
  * [8] Lingxiao Zhao and Leman Akoglu.  Pairnorm: Tackling oversmoothing in gnns.  arXiv preprint arXiv:1909.12223, 2019.
  * [9] Deli Chen, Yankai Lin, Wei Li, Peng Li, Jie Zhou, and Xu Sun.  Measuring and relieving the over-smoothing problem for graph neural networks from the topological view.  In Proceedings of the AAAI conference on artificial intelligence, volume 34, pages 3438–3445, 2020.
  * [10] Felix Wu, Amauri Souza, Tianyi Zhang, Christopher Fifty, Tao Yu, and Kilian Weinberger.  Simplifying graph convolutional networks.  In International conference on machine learning, pages 6861–6871. PMLR, 2019.
  * [11] Qimai Li, Zhichao Han, and Xiao-Ming Wu.  Deeper insights into graph convolutional networks for semi-supervised learning.  In Proceedings of the AAAI conference on artificial intelligence, volume 32, 2018.
  * [12] Moshe Eliasof, Lars Ruthotto, and Eran Treister.  Improving graph neural networks with learnable propagation operators.  In International Conference on Machine Learning, pages 9224–9245. PMLR, 2023.
  * [13] Moshe Eliasof, Eldad Haber, and Eran Treister.  PDE-GCN: Novel Architectures for Graph Neural Networks Motivated by Partial Differential Equations.  Advances in Neural Information Processing Systems, 5(NeurIPS):3836–3849, 2021.
  * [14] Joan Bruna and Arthur Szlam.  Spectral Networks and Deep Locally Connected Networks on Graphs arXiv : 1312 . 6203v3 [ cs . LG ] 21 May 2014.  pages 1–14.
  * [15] Michaël Defferrard, Xavier Bresson, and Pierre Vandergheynst.  Convolutional Neural Networks on Graphs with Fast Localized Spectral Filtering.  Advances in Neural Information Processing Systems, pages 3844–3852, jun 2016.
  * [16] Thomas N. Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  5th International Conference on Learning Representations, ICLR 2017 - Conference Track Proceedings, pages 1–14, 2017.
  * [17] Ben Chamberlain, James Rowbottom, Maria I Gorinova, Michael Bronstein, Stefan Webb, and Emanuele Rossi.  Grand: Graph neural diffusion.  In International Conference on Machine Learning, pages 1407–1418. PMLR, 2021.
  * [18] Moshe Eliasof and Eran Treister.  Diffgcn: Graph convolutional networks via differential operators and algebraic multigrid pooling.  Advances in neural information processing systems, 33:18016–18027, 2020.
  * [19] Filipe De Avila Belbute-Peres, Thomas Economon, and Zico Kolter.  Combining differentiable pde solvers and graph neural networks for fluid flow prediction.  In international conference on machine learning, pages 2402–2411. PMLR, 2020.
  * [20] Moshe Eliasof, Eldad Haber, and Eran Treister.  Adr-gnn: Advection-diffusion-reaction graph neural networks.  arXiv preprint arXiv:2307.16092, 2023.
  * [21] Mario Lino, Stathi Fotiadis, Anil A Bharath, and Chris D Cantwell.  Multi-scale rotation-equivariant graph neural networks for unsteady eulerian fluid dynamics.  Physics of Fluids, 34(8), 2022.
  * [22] Manuel Miranda and Ernesto Estrada.  Degree-biased advection–diffusion on undirected graphs/networks.  Mathematical Modelling of Natural Phenomena, 17:30, 2022.
  * [23] Ben London and Lise Getoor.  Collective classification of network data.  Data Classification: Algorithms and Applications, pages 399–416, 2014.
  * [24] Zhirong Wu, Shuran Song, Aditya Khosla, Fisher Yu, Linguang Zhang, Xiaoou Tang, and Jianxiong Xiao.  3D ShapeNets: A deep representation for volumetric shapes.  Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 07-12-June:1912–1920, 2015.
  * [25] Zhilin Yang, William Cohen, and Ruslan Salakhudinov.  Revisiting semi-supervised learning with graph embeddings.  In International conference on machine learning, pages 40–48. PMLR, 2016.
  * [26] Hongbin Pei, Bingzhe Wei, Kevin Chen-Chuan Chang, Yu Lei, and Bo Yang.  Geom-gcn: Geometric graph convolutional networks.  arXiv preprint arXiv:2002.05287, 2020.
  * [27] Federica Bogo, Javier Romero, Matthew Loper, and Michael J. Black.  FAUST: Dataset and evaluation for 3D mesh registration.  In Proceedings IEEE Conf. on Computer Vision and Pattern Recognition (CVPR), pages 3794 –3801, Columbus, Ohio, USA, June 2014.


#  Safety-Aware Preference-Based Learning for Safety-Critical Control

Ryan K. Cosner1<EMAIL_ADDRESS>
Maegan Tucker1<EMAIL_ADDRESS>
Andrew J. Taylor1<EMAIL_ADDRESS>
Kejun Li1<EMAIL_ADDRESS>
Tamás G. Molnár1<EMAIL_ADDRESS>
Wyatt Ubellacker1<EMAIL_ADDRESS>
Anil Alan2<EMAIL_ADDRESS>
Gábor Orosz2<EMAIL_ADDRESS>
Yisong Yue1<EMAIL_ADDRESS>
Aaron D. Ames1<EMAIL_ADDRESS>
1 California Institute of Technology  Pasadena  CA  USA
2 University of Michigan  Ann Arbor  MI  USA

###### Abstract

Bringing dynamic robots into the wild requires a tenuous balance between
performance and safety. Yet controllers designed to provide robust safety
guarantees often result in conservative behavior, and tuning these controllers
to find the ideal trade-off between performance and safety typically requires
domain expertise or a carefully constructed reward function. This work
presents a design paradigm for systematically achieving behaviors that balance
performance and robust safety by integrating safety-aware Preference-Based
Learning (PBL) with Control Barrier Functions (CBFs). Fusing these
concepts—safety-aware learning and safety-critical control—gives a robust
means to achieve safe behaviors on complex robotic systems in practice. We
demonstrate the capability of this design paradigm to achieve safe and
performant perception-based autonomous operation of a quadrupedal robot both
in simulation and experimentally on hardware.

###### keywords:

Preference-Based Learning, Control Barrier Functions, Safety-Critical Control,
Robotics

##  1 Introduction

The increasing demands of modern engineering problems have required a
commensurate increase in the complexity of the underlying control systems
being used. The process of designing these complex control systems is often
accomplished by separating the design into individual subsystems such as
sensing, planning, and low-level control, which are later integrated. A
principal challenge in the integration of such complex systems is balancing
safety with performance at the system level. When each individual subsystem is
designed using over-approximations of worst-case scenarios, the integrated
system becomes extremely conservative and exhibits poor performance
(Singletary et al., 2021; Alan et al., 2022). The commonly employed
alternative is tuning the safety-performance trade-off of each component to
achieve the desired system-level behavior (Ma et al., 2017), which can be
challenging even for domain experts as the tuning is often done via
qualitative assessments.

For instance, for complex safety-critical systems, Control Barrier Functions
(CBFs) have become a popular tool for the constructive synthesis of model-
based controllers that endow nonlinear systems with rigorous guarantees of
safety (Ames et al., 2014, 2019; Hobbs et al., 2021). As these safety
guarantees are susceptible to inaccuracies in the models of a system’s
dynamics, actuators, and sensors, approaches have been proposed to deal with
model uncertainty (Wang et al., 2018; Taylor and Ames, 2020; Castañeda et al.,
2020; Taylor et al., 2020), disturbances (Jankovic, 2018; Kolathaya and Ames,
2018; Clark, 2019; Santoyo et al., 2019; Alan et al., 2022; Choi et al.,
2021), and measurement errors (Takano and Yamakita, 2018; Dean et al., 2020;
Cosner et al., 2021). These approaches can work well when deployed
independently, but can often lead to extremely conservative systems with poor
performance when used in conjunction. In practice, achieving performant system
behaviors with these methods is accomplished by conceding theoretical safety
guarantees for over-approximations of worst-case scenarios and tuning
controller robustness parameters.

Figure 1:  An overview of the design paradigm with Safety-Aware Preference-
Based Learning.

To reduce the burden on experts in controller tuning, we seek to incorporate
Preference-Based Learning (PBL) into the design of complex safety-critical
control systems. PBL has shown to be a powerful tool for converting subjective
user preferences of system behavior (e.g., behavior A has a superior
performance-safety trade-off than behavior B) into quantitative adjustments to
design parameters. The main advantage of online PBL is its ability to
interactively infer a user’s latent utility function using only subjective
feedback such as pairwise preferences and ordinal labels (Yue et al., 2012;
Shivaswamy and Joachims, 2012). This methodology has been demonstrated in
application for exoskeleton gait optimization (Tucker et al., 2020b), bipedal
locomotion (Tucker et al., 2021), spinal cord stimulation (Sui et al., 2018),
trajectory planning (Sadigh et al., 2017; Bıyık et al., 2020; Jain et al.,
2015), search engines (Raman et al., 2013), and recommender systems (De Gemmis
et al., 2009). For applications with actions that may be classified as safe or
unsafe, safety-critical PBL algorithms have been demonstrated to prevent
unsafe actions from being sampled (Sui et al., 2015, 2018; Berkenkamp et al.,
2016). However, these safety-critical algorithms require approximating worse-
case scenarios, and over-approximations may yield the characterization of
performant and safe actions as catastrophically unsafe. Thus, we seek to
formulate a safety-aware approach to PBL that balances avoiding sampling
unsafe actions with being overly conservative.

In this work we propose a design paradigm for synthesizing performant and
robust safety-critical controllers on real systems via safety-aware online PBL
(illustrated in Fig. 1). The contributions of this work are threefold. First,
we propose Safety-Aware LineCoSpar (SA-LineCoSpar), a modified version of
LineCoSpar (Tucker et al., 2020a) capable of high-dimensional preference-based
Bayesian optimization while also accounting for safety. Second, we combine the
robustness properties of Measurement-Robust CBFs (MR-CBFs) (Dean et al., 2020)
to measurement uncertainty and Input-to-State Safe CBFs (ISSf-CBFs) (Kolathaya
and Ames, 2018) to disturbances with reduced-order multi-layer safety-critical
control (Molnar et al., 2021) to achieve provable safety guarantees in a
parametric form amenable to SA-LineCoSpar. Finally, we deploy these two
methods together as a design paradigm for a safety-critical controller on a
quadrupedal robot in simulation and on hardware in laboratory and outdoor
settings. Additionally, this work is the first time that PBL has been used to
tune a CBF-based controller, and the first time these CBF methods have been
combined.

##  2 Safety-Aware Preference-Based Learning

Preference-Based Learning (PBL) provides an approach for searching complex
parameter spaces via subjective feedback, eliminating the need to explicitly
define an underlying reward function. This is particularly relevant for
safety-critical systems, as quantifying the user-preferred trade-off between
robustness and performance is generally difficult. Moreover, poorly defined
reward functions often result in “reward hacking” (Amodei et al., 2016), in
which undesirable actions achieve high rewards. In this work, we propose
Safety-Aware LineCoSpar (SA-LineCoSpar), outlined in Alg. 2. This algorithm is
a modification of the LineCoSpar algorithm (Tucker et al., 2020a), which
iteratively selects actions to query user for subjective feedback and updates
its belief of the user’s underlying utility function via Bayesian inference
accordingly.

Problem Setup: Let $\mathbf{a}$ denote an action, such as a collection of $l$
parameters used in a feedback controller, that takes values in a finite search
space $A\subset\mathbb{R}^{l}$. We assume that each action $\mathbf{a}\in A$
has an unknown utility to the user, defined by a function
$r(\mathbf{a}):A\to\mathbb{R}$. These utilities can be written in a vectorized
form
$\mathbf{r}_{A}=[r(\mathbf{a}_{1}),\dots,r(\mathbf{a}_{|A|})]^{\top}\in\mathbb{R}^{|A|}$.
In each iteration, $s\in\mathbb{N}$ actions are sampled from $A$. After these
actions are executed, the user is queried for two forms of feedback: pairwise
preferences and ordinal labels, in this work describing performance and
safety, respectively. This feedback is collected into the combined dataset
$D$.

Modeling the Utility Function: We model the unknown utilities $\mathbf{r}_{A}$
as a Gaussian process with prior
$\mathcal{N}(\bm{0},\boldsymbol{\Sigma}^{\text{pr}})$, where each element of
the covariance matrix
$\boldsymbol{\Sigma}^{\text{pr}}\in\mathbb{S}^{|A|\times|A|}_{\succ 0}$ is
computed as
$\boldsymbol{\Sigma}^{\text{pr}}_{ij}={k}(\mathbf{a}_{i},\mathbf{a}_{j})$ with
a kernel function $k:A\times A\to\mathbb{R}$ and $\mathbf{a}_{i}\in A$
denoting the $i^{\textup{th}}$ action in $A$. We select $k$ to be the squared
exponential kernel, yielding a prior given by the multivariate Gaussian:

$\displaystyle\mathcal{P}(\mathbf{r}_{A})=\frac{1}{(2\pi)^{\lvert
A\rvert/2}\lvert\boldsymbol{\Sigma}^{\text{pr}}\rvert^{1/2}}\text{exp}\left(-\frac{1}{2}\mathbf{r}_{A}^{\top}{\left(\boldsymbol{\Sigma}^{\text{pr}}\right)}^{-1}\mathbf{r}_{A}\right).$
(1)

Given a dataset $D$, the posterior is proportional to the likelihood and the
prior by Bayes’ theorem, i.e., $\mathcal{P}(\mathbf{r}_{A}\mid
D)\propto\mathcal{P}(D\mid\mathbf{r}_{A})\mathcal{P}(\mathbf{r}_{A})$. We
denote the maximum a posterior (MAP) estimate of the posterior by
$\hat{\mathbf{r}}_{A}\in\mathbb{R}^{|A|}$, which is defined as
$\hat{\mathbf{r}}_{A}\triangleq\operatornamewithlimits{argmax}_{\mathbf{r}_{A}\in\mathbb{R}^{|A|}}\mathcal{P}(\mathbf{r}_{A}\mid
D)$, noting that $\hat{\mathbf{r}}_{A}$ is equivalent to the minimizer of
$\mathcal{S}(\mathbf{r}_{A})=-\ln(\mathcal{P}(D\mid\mathbf{r}_{A}))+\frac{1}{2}\mathbf{r}_{A}^{T}\left(\Sigma^{\text{pr}}\right)^{-1}\mathbf{r}_{A}$.
We approximate the posterior as a multivariate Gaussian centered at
$\hat{\mathbf{r}}_{A}$ with the covariance
$\boldsymbol{\Sigma}_{A}\in\mathbb{S}^{|A|\times|A|}_{\succ 0}$ defined as
$\boldsymbol{\Sigma}_{A}=(\frac{\partial^{2}\mathcal{S}}{\partial\mathbf{r}_{A}^{2}}(\hat{\mathbf{r}}_{A}))^{-1}$.
This is known as the Laplace approximation of the distribution
$\mathcal{P}(\mathbf{r}_{A}\mid D)$, i.e., $\mathcal{P}(\mathbf{r}_{A}\mid
D)\approx\mathcal{N}(\hat{\mathbf{r}}_{A},\boldsymbol{\Sigma}_{A})$ (Chu and
Ghahramani, 2005). Lastly, since this process is intractable for high-
dimensional action spaces, we perform it over a subset $S\subset A$ of the
action space, forming a partial characterization of the utilities denoted by
$\mathcal{P}(\mathbf{r}_{S}\mid
D)\approx\mathcal{N}(\hat{\mathbf{r}}_{S},\boldsymbol{\Sigma}_{S})$, with
$\mathbf{r}_{S},\hat{\mathbf{r}}_{S}\in\mathbb{R}^{|S|}$.

Preference Likelihood Function: A pairwise preference is defined as a relation
between two actions $\mathbf{a}_{1},\mathbf{a}_{2}\in A$, where
$\mathbf{a}_{1}\succ\mathbf{a}_{2}$ if action $\mathbf{a}_{1}$ is preferred to
$\mathbf{a}_{2}$. Since user preferences are expected to be corrupted by
noise, we model individual pairwise preferences via a likelihood function:

$\mathcal{P}(\mathbf{a}_{1}\succ\mathbf{a}_{2}|r(\mathbf{a}_{1}),r(\mathbf{a}_{2}))=g_{p}\left(\frac{r(\mathbf{a}_{1})-r(\mathbf{a}_{2})}{c_{p}}\right)$
(2)

where $g_{p}:\mathbb{R}\to[0,1]$ is any monotonically-increasing link
function, and $c_{p}\in\mathbb{R}_{>0}$ accounts for preference noise. We
select $g_{p}$ to be the sigmoid function, i.e., $g_{p}(x)=1/(1+e^{-x})$.
Assuming conditional independence, the likelihood function for a collection of
$K\in\mathbb{N}$ preferences, $D_{p}$, can be modeled as the product of each
individual preference likelihood:

$\mathcal{P}(D_{p}|r(\mathbf{a}_{11}),r(\mathbf{a}_{12}),\cdots,r(\mathbf{a}_{K2}))=\prod_{k=1}^{K}\mathcal{P}(\mathbf{a}_{k1}\succ\mathbf{a}_{k2}|r(\mathbf{a}_{k1}),r(\mathbf{a}_{k2})),$
(3)

where $\mathbf{a}_{k1},\mathbf{a}_{k2}\in A$ are the preferred and non-
preferred actions, respectively, in the $k^{th}$ preference.

Ordinal Likelihood Function: We partition the action space into “unsafe” and
“safe” actions by leveraging the ordinal nature of these definitions (i.e.,
unsafe actions are always considered worse than safe actions). A user provides
this feedback as an ordinal label, which assigns an action to a discrete
ordered category such as “bad” and “good” (Chu et al., 2005). While ordinal
labels can be generalized to any number of ordinal categories (c.f. Li et al.
(2021)), we utilize just two categories to represent “unsafe” and “safe”. In
this case, the action space is decomposed into two disjoint sets, $A=O_{1}\cup
O_{2}$, with $\mathbf{a}\in O_{1}$ if $r(\mathbf{a})<\beta$ and $\mathbf{a}\in
O_{2}$ if $r(\mathbf{a})\geq\beta$, with the ordinal threshold
$\beta\in\mathbb{R}$. As with preferences, we assume that ordinal label
feedback is corrupted by noise and is modeled as:

$\displaystyle\mathcal{P}(\mathbf{a}\in O_{1}\mid
r(\mathbf{a}))=g_{o}\left(\frac{\beta-r(\mathbf{a})}{c_{o}}\right),\qquad\mathcal{P}(\mathbf{a}\in
O_{2}\mid
r(\mathbf{a}))=1-g_{o}\left(\frac{\beta-r(\mathbf{a})}{c_{o}}\right),$ (4)

where $g_{o}:\mathbb{R}\to[0,1]$ is any monotonically-increasing link function
and $c_{o}$ quantifies the noise in the ordinal label feedback. Again, we
select $g_{o}$ to be the sigmoid function $g_{o}(x)=1/(1+e^{-x})$. Assuming
conditional independence of ordinal label queries, the likelihood function for
a collection of $M\in\mathbb{N}$ ordinal labels, $D_{o}$, can be modeled as
the product of each individual ordinal likelihood:

$\displaystyle\mathcal{P}(D_{o}\mid
r(\mathbf{a}_{1}),\cdots,r(\mathbf{a}_{k}))=$
$\displaystyle\prod_{k=1}^{M}\mathcal{P}\left(\mathbf{a}_{k}\in O_{o(k)}\mid
r(\mathbf{a}_{k})\right),$ (5)

where $\mathbf{a}_{k}\in A$ refers to the action corresponding to the
$k^{\textup{th}}$ ordinal label, $o(k)\in\\{1,2\\}$. For our simulation and
experiments, the hyperparameters $c_{p}$, $c_{o}$, $\beta$ are determined in
advance. Lastly, assuming conditional independence of the feedback mechanisms,
the combined likelihood function is calculated as the product of the
individual likelihoods, $\mathcal{P}(D\mid r)=\mathcal{P}(D_{p}\mid
r)\mathcal{P}(D_{o}\mid r)$.

Sampling New Actions: In the first iteration ($i=1$), $s\in\mathbb{N}$ actions
are sampled randomly from $A$, recorded as the set of visited actions
$V_{1}=\\{\mathbf{a}_{1}^{(1)},\dots,\mathbf{a}_{1}^{(s)}\\}$, executed on the
system, and the preferences and ordinal labels are collected into a dataset
$D_{1}$. In each subsequent iteration ($i>1$), $s$ new actions are sampled
using Thompson sampling, which is shown to have desirable regret minimization
properties (Chapelle and Li, 2011). Ideally, Thompson sampling draws $s$
samples from the posterior $\mathcal{P}(\mathbf{r}_{A}\mid D_{i-1})$, i.e
$\mathbf{r}^{(j)}\sim\mathcal{P}(\mathbf{r}_{A}\mid D_{i-1})$ for
$j\in\\{1,\dots,s\\}$, and the action $\mathbf{a}_{i}^{(j)}\in A$ maximizing
each $\mathbf{r}^{(j)}$ is selected to execute on the system. These sampled
actions $\\{\mathbf{a}_{i}^{(1)},\dots,\mathbf{a}_{i}^{(s)}\\}$ are
concatenated with $V_{i-1}$ to produce $V_{i}$, executed on the system, and
the resulting preferences and ordinal labels are concatenated with $D_{i-1}$
to produce $D_{i}$. However, since it is intractable to approximate
$\mathcal{P}(\mathbf{r}_{A}\mid D)$ for high-dimensional action spaces, we
utilize a dimensionality-reduction technique introduced in Tucker et al.
(2020a) that instead updates the posterior over a subset $S_{i}\subset A$.
Motivated by Kirschner et al. (2019), we construct the subset as
$S_{i}=L_{i}\cup V_{i-1}$, where $L_{i}\subset A$ is the collection of
$e\in\mathbb{N}$ actions in $A$ closest to a randomly drawn line
$\ell_{i}\subset\mathbb{R}^{l}$. This line is drawn to intersect with the
believed best action, computed as
$\hat{\mathbf{a}}_{i-1}^{*}=\operatornamewithlimits{argmax}_{\mathbf{a}\in
V_{i-1}}\hat{\mathbf{r}}_{V_{i-1}}(\mathbf{a})$ where
$\hat{\mathbf{r}}_{V_{i-1}}$ is the MAP estimate of the posterior
$\mathcal{P}(\mathbf{r}_{V_{i-1}}\mid D_{i})$. See Tucker et al. (2020a) for
more details.

Safety-Aware Sampling: It is important to avoid unsafe actions during
sequential decision making in certain applications, such as learning robotic
controllers on hardware, where low-reward actions might lead to physical
damage of the platform. Safe exploration algorithms (Sui et al., 2015, 2018;
Berkenkamp et al., 2016) considered the setting where actions below a
prespecified safety threshold are catastrophic and must be avoided at all
cost. In our work, since we construct controllers that account for safety, we
adopt a more optimistic learning approach called safety-aware. In this case,
actions labeled by a human as “unsafe” are not catastrophic but undesirable.
Thus, the algorithm avoids these actions; whereas the safe exploration
algorithms guarantee that no such actions are sampled which can be sometimes
exceedingly conservative in settings like ours.

To achieve this safety-awareness, we leverage the approach introduced in Li et
al. (2021), which uses ordinal labels to identify a region of interest (ROI)
in $A$. In this work, the ROI is defined to be the actions labeled as “safe”.
In each iteration $i$ we estimate an ROI within the set $S_{i}$ as:

$\displaystyle S^{\text{ROI}}_{i}=\\{\mathbf{a}\in
S_{i}\mid\hat{\mathbf{r}}_{S_{i}}(\mathbf{a})+\lambda\boldsymbol{\sigma}_{S_{i}}(\mathbf{a})>\beta\\},$
(6)

where $\hat{\mathbf{r}}_{S_{i}}(\mathbf{a})$ and
$\boldsymbol{\sigma}_{S_{i}}(\mathbf{a})$ are the posterior mean and standard
deviation, respectively, evaluated at the action $\mathbf{a}\in S_{i}$. The
variable $\lambda\in\mathbb{R}$ determines how conservative the algorithm
would be in estimating the safety region, as illustrated in Figure 2. We see
that lower values of $\lambda$ result in fewer unsafe actions being sampled,
with only a slight effect on sample-efficiency. The restriction to
$S_{i}^{\text{ROI}}$ is added to LineCoSpar by only considering actions in
$S_{i}^{\text{ROI}}$ during Thompson sampling. We refer to this as Safety-
Aware LineCoSpar (SA-LineCoSpar), with the full algorithm outlined in Alg. 2.

[H] Safety-Aware LineCoSpar Inputinput $s$ uniform random actions
($V_{1}\subset A$), corresponding feedback ($D_{1}$),  for $i=2,\dots,N$ do
Update posterior over $V_{i-1}$
$\hat{\mathbf{a}}^{*}_{i-1}\leftarrow\operatornamewithlimits{argmax}_{\mathbf{a}\in
V_{i-1}}\hat{\mathbf{r}}_{V_{i-1}}(\mathbf{a})$ $L_{i}\leftarrow$ New linear
subspace intersecting $\hat{\mathbf{a}}^{*}_{i-1}$ Construct subspace
$S_{i}=L_{i}\cup V_{i-1}$ Update the model posterior over $S_{i}$ Determine
region of interest $S^{\text{ROI}}_{i}$ for $j=1,\dots,s$ do
$r^{(j)}\sim\mathcal{N}(\hat{\mathbf{r}}_{S_{i}},\boldsymbol{\Sigma}_{S_{i}})$
$\mathbf{a}_{i}^{(j)}\leftarrow\operatornamewithlimits{argmax}_{\mathbf{a}\in
S^{\text{ROI}}_{i}}r^{(j)}$ Deploy
$\\{\mathbf{a}_{i}^{(1)},\dots\mathbf{a}_{i}^{(s)}\\}$ on system
$V_{i}\leftarrow
V_{i-1}\cup\\{\mathbf{a}_{i}^{(1)},\dots\mathbf{a}_{i}^{(s)}\\}$
$D_{i}\leftarrow D_{i-1}\cup\text{new prefs.}\cup\text{new ord. labels}$
Figure 2: A comparison of SA-LineCoSpar and standard LineCoSpar on a synthetic
utility function (drawn from the Gaussian prior) averaged over 50 runs with
standard error shown by the shaded region. The safety-aware criteria reduces
the number of sampled unsafe actions with a minimal effect on the prediction
error, defined as $\lvert\hat{\mathbf{a}}_{i}^{*}-\mathbf{a}^{*}\rvert$ with
$\hat{\mathbf{a}}_{i}^{*}\triangleq\operatornamewithlimits{argmax}_{\mathbf{a}}\hat{\mathbf{r}}_{S_{i}}$
and
$\mathbf{a}^{*}\triangleq\operatornamewithlimits{argmax}_{\mathbf{a}}r(\mathbf{a})$.

##  3 Robust Safety-Critical Control

In this section, we formalize robust safety and discuss safe controller
synthesis through the use of Control Barrier Functions (CBFs), that ultimately
yield controllers whose parameters are to be updated with SA-LineCoSpar.
Consider the following nonlinear control-affine system:

$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x})+\mathbf{g}(\mathbf{x})(\mathbf{v}+\mathbf{d}(t)),$
(7)

with state ${\mathbf{x}\in\mathbb{R}^{n}}$, input
${\mathbf{v}\in\mathbb{R}^{m}}$, functions
${\mathbf{f}:\mathbb{R}^{n}\to\mathbb{R}^{n}}$ and
${\mathbf{g}:\mathbb{R}^{n}\to\mathbb{R}^{n\times m}}$ assumed to be locally
Lipschitz continuous on their domains, and piecewise continuous disturbance
signal $\mathbf{d}:\mathbb{R}_{\geq 0}\to\mathbb{R}^{m}$ for which we define
$\|\mathbf{d}\|_{\infty}\triangleq\sup_{t\geq 0}\|\mathbf{d}(t)\|$. Specifying
the input via a controller $\mathbf{k}:\mathbb{R}^{n}\to\mathbb{R}^{m}$ that
is locally Lipschitz continuous on its domain yields the closed-loop system:

$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x})+\mathbf{g}(\mathbf{x})(\mathbf{k}(\mathbf{x})+\mathbf{d}(t)).$
(8)

We assume for any initial condition
${\mathbf{x}(0)=\mathbf{x}_{0}\in\mathbb{R}^{n}}$ and disturbance
$\mathbf{d}$, this system has a unique solution $\mathbf{x}(t)$ for all
$t\in\mathbb{R}_{\geq 0}$. We consider this system safe if its state
$\mathbf{x}(t)$ remains in a safe set $\mathcal{C}\subset\mathbb{R}^{n}$,
defined as the 0-superlevel set of a continuously differentiable function
$h:\mathbb{R}^{n}\times\mathbb{R}^{p}\to\mathbb{R}$:

$\mathcal{C}=\\{\mathbf{x}\in\mathbb{R}^{n}:h(\mathbf{x},\boldsymbol{\rho})\geq
0\\},$ (9)

where $\boldsymbol{\rho}\in\mathbb{R}^{p}$ are constant application-specific
parameters. We say the set $\mathcal{C}\subset\mathbb{R}^{n}$ is forward
invariant if for every $\mathbf{x}_{0}\in\mathcal{C}$ the solution
$\mathbf{x}(t)$ to (8) satisfies $\mathbf{x}(t)\in\mathcal{C}$ for all $t\geq
0$. The system (8) is safe with respect to $\mathcal{C}$ if $\mathcal{C}$ is
forward invariant. Ensuring the safety of the set $\mathcal{C}$ in the absence
of disturbances can be achieved through Control Barrier Functions (CBFs):

######  Definition 1 (Control Barrier Functions (CBF) (Ames et al., 2014))

The function $h$ is a Control Barrier Function (CBF) for (7) on $\mathcal{C}$
if there exists $\alpha\in\mathcal{K}_{\infty}^{\rm e}$111We say that a
continuous function ${\alpha:\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}}$ is
class $\mathcal{K}_{\infty}$ (${\alpha\in\mathcal{K}_{\infty}}$) if
${\alpha(0)=0}$, $\alpha$ is strictly monotonically increasing, and
$\lim_{r\to\infty}\alpha(r)=\infty$. We say that a continuous function
$\alpha:\mathbb{R}\to\mathbb{R}$ is class $\mathcal{K}_{\infty}^{\rm e}$
(${\alpha\in\mathcal{K}_{\infty}^{\rm e}}$) if ${\alpha(0)=0}$, $\alpha$ is
strictly monotonically increasing, $\lim_{r\to\infty}\alpha(r)=\infty$, and
$\lim_{r\to-\infty}\alpha(r)=-\infty$. such that for all
$\mathbf{x}\in\mathbb{R}^{n}$:

$\sup_{\mathbf{v}\in\mathbb{R}^{m}}\underbrace{\frac{\partial
h}{\partial\mathbf{x}}(\mathbf{x},\boldsymbol{\rho})\mathbf{f}(\mathbf{x})}_{L_{\mathbf{f}}h(\mathbf{x},\boldsymbol{\rho})}+\underbrace{\frac{\partial
h}{\partial\mathbf{x}}{(\mathbf{x},\boldsymbol{\rho})}\mathbf{g}(\mathbf{x})}_{L_{\mathbf{g}}h(\mathbf{x},\boldsymbol{\rho})}\mathbf{v}>-\alpha(h(\mathbf{x},\boldsymbol{\rho})).$
(10)

While it may be possible to synthesize controllers that render a given set
$\mathcal{C}$ safe in the presence of disturbances (Jankovic, 2018), this may
result in overly-conservative behavior. Instead, we consider how safety
properties degrade with disturbances via the following definition.

######  Definition 2 (Input-to-State Safety (Kolathaya and Ames, 2018))

The system (8) is Input-to-

State Safe (ISSf) with respect to $\mathcal{C}$ if there exists
$\gamma\in\mathcal{K}_{\infty}$ such that for all $\delta\in\mathbb{R}_{\geq
0}$ and disturbances $\mathbf{d}:\mathbb{R}_{\geq 0}\to\mathbb{R}^{m}$
satisfying $\|\mathbf{d}\|_{\infty}\leq\delta$, the set
$\mathcal{C}_{\delta}\subset\mathbb{R}^{n}$ defined as:

$\displaystyle\mathcal{C}_{\delta}=\\{\mathbf{x}\in\mathbb{R}^{n}:h(\mathbf{x},\boldsymbol{\rho})\geq-\gamma(\delta)\\}$
(11)

is forward invariant. The function $h$ is an Input-to-State Safe Control
Barrier Function (ISSf-CBF) for (7) on $\mathcal{C}$ with parameter
$\varphi\in\mathbb{R}_{\geq 0}$ if there exists
$\alpha\in\mathcal{K}_{\infty}^{\rm e}$ such that for all
$\mathbf{x}\in\mathbb{R}^{n}$:

$\sup_{\mathbf{v}\in\mathbb{R}^{m}}L_{\mathbf{f}}h(\mathbf{x},\boldsymbol{\rho})+L_{\mathbf{g}}h(\mathbf{x},\boldsymbol{\rho})\mathbf{v}-\varphi\|L_{\mathbf{g}}h(\mathbf{x},\boldsymbol{\rho})\|^{2}>-\alpha(h(\mathbf{x},\boldsymbol{\rho})).$
(12)

The parameter $\boldsymbol{\rho}\in\mathbb{R}^{p}$ contains information about
the system’s environment that designates a given state as safe or unsafe, such
as the location and size of obstacles. In highly unstructured environments the
system may need to estimate $\boldsymbol{\rho}$, denoted by
$\widehat{\boldsymbol{\rho}}\in\mathbb{R}^{p}$ from measurements of the
environment, such as camera data. The process of converting complex
measurements into the representation $\widehat{\boldsymbol{\rho}}$ is often
imperfect, leading to error between the estimator
$\widehat{\boldsymbol{\rho}}$ and the true value of $\boldsymbol{\rho}$. In
this setting, safety can be achieved via Measurement-Robust Control Barrier
Functions (MR-CBFs):

######  Definition 3 (Measurement-Robust Control Barrier Functions (Dean et
al., 2020))

The

function $h$ is a Measurement-Robust Control Barrier Function (MR-CBF) for (7)
on $\mathcal{C}$ with parameters $a,b\in\mathbb{R}_{\geq 0}$ if there exists
$\alpha\in\mathcal{K}_{\infty}^{\rm e}$ such that for all
$\widehat{\boldsymbol{\rho}}\in\mathbb{R}^{p}$ and
$\mathbf{x}\in\mathbb{R}^{n}$:

$\sup_{\mathbf{v}\in\mathbb{R}^{m}}L_{\mathbf{f}}h(\mathbf{x},\widehat{\boldsymbol{\rho}})+L_{\mathbf{g}}h(\mathbf{x},\widehat{\boldsymbol{\rho}})\mathbf{v}-a-b\|\mathbf{v}\|>-\alpha(h(\mathbf{x},\widehat{\boldsymbol{\rho}})).$
(13)

The following theorem summarizes the safety results achieved with these
various types of CBFs:

######  Theorem 3.1.

Consider the set $\mathcal{C}$ defined in (9).

  1. 1.

If $h$ is a CBF for (7) on $\mathcal{C}$, $\mathbf{d}(t)=\mathbf{0}$ for
$t\in\mathbb{R}_{\geq 0}$ and $\widehat{\boldsymbol{\rho}}=\boldsymbol{\rho}$,
then there exists a controller $\mathbf{k}$ such that (8) is safe with respect
to $\mathcal{C}$.

  2. 2.

If $h$ is an ISSf-CBF for (7) on $\mathcal{C}$ with parameter $\varphi$ and
$\widehat{\boldsymbol{\rho}}=\boldsymbol{\rho}$, then there exists a
controller $\mathbf{k}$ such that (8) is ISSf with respect to $\mathcal{C}$
with $\gamma(\delta)=-\alpha^{-1}(-\delta^{2}/4\varphi)$ where
$\alpha^{-1}\in\mathcal{K}_{\infty}^{\rm e}$.

  3. 3.

Assume $L_{\mathbf{f}}h$, $L_{\mathbf{g}}h$, and $\alpha\circ h$ are Lipschitz
continuous on their domains, and assume that
$\|\widehat{\boldsymbol{\rho}}-\boldsymbol{\rho}\|\leq\epsilon$ for some
$\epsilon\in\mathbb{R}_{\geq 0}$. Then there exists
$\underline{a},\underline{b}\in\mathbb{R}_{\geq 0}$ such that if $h$ is an MR-
CBF for (7) on $\mathcal{C}$ with parameters $a,b\in\mathbb{R}_{\geq 0}$
satisfying $a\geq\underline{a}$ and $b\geq\underline{b}$, and
$\mathbf{d}(t)=\mathbf{0}$ for $t\in\mathbb{R}_{\geq 0}$, then there exists a
controller $\mathbf{k}$ such that (8) is safe with respect to $\mathcal{C}$.

##  4 Integrating Safety-Aware Preference-Based Learning with Safety-Critical
Control

In this section we propose a design paradigm that leverages SA-LineCoSpar to
select parameters for a CBF-based controller that achieves performance and
safety for a multi-layered control system.

Multi-Layered System Dynamics: Many real-life engineering systems have high-
dimensional state spaces and complex dynamics. Hence control systems are often
designed as a set of interconnected subsystems, such as a low-dimensional
subsystem that provides reference signals capturing safe behavior and a high-
dimensional subsystem that tracks these reference signals. In particular,
consider the following cascaded nonlinear control-affine system resulting as a
modification of (7):

$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x})+\mathbf{g}(\mathbf{x})\boldsymbol{\kappa}(\boldsymbol{\xi}),\qquad\dot{\boldsymbol{\xi}}=\mathbf{f}_{\boldsymbol{\xi}}(\mathbf{x},\boldsymbol{\xi})+\mathbf{g}_{\boldsymbol{\xi}}(\mathbf{x},\boldsymbol{\xi})\mathbf{u},$
(14)

with additional states $\boldsymbol{\xi}\in\mathbb{R}^{n_{\boldsymbol{\xi}}}$,
control input $\mathbf{u}\in\mathbb{R}^{m_{\boldsymbol{\xi}}}$ and functions
$\boldsymbol{\kappa}:\mathbb{R}^{n_{\boldsymbol{\xi}}}\to\mathbb{R}^{m}$,
$\mathbf{f}_{\boldsymbol{\xi}}:\mathbb{R}^{n}\times\mathbb{R}^{n_{\boldsymbol{\xi}}}\to\mathbb{R}^{n_{\boldsymbol{\xi}}}$,
and
$\mathbf{g}_{\boldsymbol{\xi}}:\mathbb{R}^{n}\times\mathbb{R}^{n_{\boldsymbol{\xi}}}\to\mathbb{R}^{n_{\boldsymbol{\xi}}\times
m_{\boldsymbol{\xi}}}$ assumed to be locally Lipschitz continuous on their
domains. We note that the input $\mathbf{v}$ from (7) was replaced by
$\boldsymbol{\kappa}(\boldsymbol{\xi})$. These dynamics may represent Euler-
Lagrange systems such as robots, where $\mathbf{x}$ reflects base position,
$\boldsymbol{\xi}$ captures base velocities and joint positions and
velocities, and the input $\mathbf{u}$ reflects the torques applied to the
joints.

Given this cascaded system, we utilize the low-dimensional subsystem to ensure
that $\mathcal{C}$ is ISSf by making two assumptions. First, we assume the
safe set $\mathcal{C}$ can be described as in (9), such that it only depends
on the states $\mathbf{x}$ and parameters $\boldsymbol{\rho}$, and not the
states $\boldsymbol{\xi}$. For example, in the context of a robotic system,
this assumption is justified if safety is described as keeping the base
position of the robot away from obstacles. Second, we assume there exists a
controller
$\boldsymbol{\pi}:\mathbb{R}^{n}\times\mathbb{R}^{n_{\boldsymbol{\xi}}}\times\mathbb{R}^{m}\to\mathbb{R}^{m_{\boldsymbol{\xi}}}$
and $\mu\in\mathbb{R}_{\geq 0}$ such that for any continuous, bounded signal
$\mathbf{s}:\mathbb{R}_{\geq 0}\to\mathbb{R}^{m}$, the closed-loop system:

$\dot{\boldsymbol{\xi}}=\mathbf{f}_{\boldsymbol{\xi}}(\mathbf{x},\boldsymbol{\xi})+\mathbf{g}_{\boldsymbol{\xi}}(\mathbf{x},\boldsymbol{\xi})\boldsymbol{\pi}(\mathbf{x},\boldsymbol{\xi},\mathbf{s}(t))$
(15)

satisfies the following implication:

$\|\boldsymbol{\kappa}(\boldsymbol{\xi}(0))-\mathbf{s}(0)\|\leq\mu\implies\|\boldsymbol{\kappa}(\boldsymbol{\xi}(t))-\mathbf{s}(t)\|\leq\mu,\quad
t\in\mathbb{R}_{\geq 0}.$ (16)

This assumption reflects that a separate controller may be designed for the
high-dimensional dynamics to track well-behaved reference signals synthesized
via the low-dimensional model. In particular, if a continuous controller
$\mathbf{k}:\mathbb{R}^{n}\to\mathbb{R}^{m}$ is designed for the low-
dimensional system (7) and
$\|\boldsymbol{\kappa}(\boldsymbol{\xi}(0))-\mathbf{k}(\mathbf{x}(0))\|\leq\mu$,
then we have that the controller $\boldsymbol{\pi}$ ensures
$\|\boldsymbol{\kappa}(\boldsymbol{\xi}(t))-\mathbf{k}(\mathbf{x}(t))\|\leq\mu$
for $t\in\mathbb{R}_{\geq 0}$. With this assumption in mind, we may study the
ISSf behavior of the closed-loop system:

$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x})+\mathbf{g}(\mathbf{x})(\mathbf{k}(\mathbf{x})+\mathbf{d}(t)),\qquad\dot{\boldsymbol{\xi}}=\mathbf{f}_{\boldsymbol{\xi}}(\mathbf{x},\boldsymbol{\xi})+\mathbf{g}_{\boldsymbol{\xi}}(\mathbf{x},\boldsymbol{\xi})\boldsymbol{\pi}(\mathbf{x},\boldsymbol{\xi},\mathbf{k}(\mathbf{x})),$
(17)

with the disturbance defined as
$\mathbf{d}(t)=\boldsymbol{\kappa}(\boldsymbol{\xi}(t))-\mathbf{k}(\mathbf{x}(t))$
satisfying $\|\mathbf{d}\|_{\infty}\leq\mu$.

Combined Robust CBFs for PBL: We now combine the robustness properties of MR-
CBFs and ISSf-CBFs to account for measurement uncertainty and the disturbance,
$\mathbf{d}$, allowing us to make robust safety guarantees for the full system
(17). This is formalized in the following theorem:

######  Theorem 4.1.

Given the set $\mathcal{C}$ defined in (9), suppose the functions
$L_{\mathbf{f}}h$, $L_{\mathbf{g}}h$, $\|L_{\mathbf{g}}h\|^{2}$, and
$\alpha\circ h$ are Lipschitz continuous on their domains, and assume that
$\|\widehat{\boldsymbol{\rho}}-\boldsymbol{\rho}\|\leq\epsilon$ for some
$\epsilon\in\mathbb{R}_{\geq 0}$. Then there exists
$\underline{a},\underline{b}\in\mathbb{R}_{\geq 0}$ such that if $h$
satisfies:

$\sup_{\mathbf{v}\in\mathbb{R}^{m}}L_{\mathbf{f}}h(\mathbf{x},\widehat{\boldsymbol{\rho}})+L_{\mathbf{g}}h(\mathbf{x},\widehat{\boldsymbol{\rho}})\mathbf{v}-\varphi\|L_{\mathbf{g}}h(\mathbf{x},\widehat{\boldsymbol{\rho}})\|^{2}-a-b\|\mathbf{v}\|>-\alpha(h(\mathbf{x},\widehat{\boldsymbol{\rho}})),$
(18)

for all $\mathbf{x}\in\mathbb{R}^{n}$ and some $a,b\in\mathbb{R}_{\geq 0}$
satisfying $a\geq\underline{a}$ and $b\geq\underline{b}$, then there exists a
controller $\mathbf{k}:\mathbb{R}^{n}\to\mathbb{R}^{m}$ such that (17) is ISSf
with respect to $\mathcal{C}$ with
$\gamma(\delta)=-\alpha^{-1}(-\delta^{2}/(4\varphi))$.

The proof of this theorem can be found in Appendix A. As in Gurriet et al.
(2018), (18) can be incorporated as a constraint into a safety filter on a
locally Lipschitz continuous nominal controller
$\mathbf{k}_{\textrm{nom}}:\mathbb{R}^{n}\to\mathbb{R}^{m}$. We call this
filter the Tunable Robustified Optimization Program (TR-OP) with tunable
parameters $\alpha,\varphi,a,$ and $b$.

$\displaystyle\mathbf{k}(\mathbf{x})={\rm
arg}\\!\\!\min_{\mathbf{v}\in\mathbb{R}^{m}}$
$\displaystyle\quad\|\mathbf{v}-\mathbf{k}_{\rm nom}(\mathbf{x})\|^{2}$ (TR-
OP) s.t. $\displaystyle
L_{\mathbf{f}}h(\mathbf{x},\widehat{\boldsymbol{\rho}}_{i})+L_{\mathbf{g}}h(\mathbf{x},\widehat{\boldsymbol{\rho}}_{i})\mathbf{v}-\varphi\|L_{\mathbf{g}}h(\mathbf{x},\widehat{\boldsymbol{\rho}}_{i})\|^{2}-a-b\|\mathbf{v}\|\geq-\alpha
h(\mathbf{x},\widehat{\boldsymbol{\rho}}_{i}),$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\forall
i\in\\{1,\dots,N_{o}\\}.$

Here we use a linear class $\mathcal{K}_{\infty}^{\rm e}$ function with
coefficient $\alpha\in\mathbb{R}_{>0}$. If we wish to enforce multiple safety
constraints, such as in obstacle avoidance with several obstacles,
$\widehat{\boldsymbol{\rho}}_{i}$ can be used to indicate the measured
parameters of the $i^{th}$ obstacle, with $N_{o}\in\mathbb{N}$ being the total
number of obstacles. Enforcing this constraint for $N_{o}>1$ can be viewed as
Boolean composition of safe sets (Glotfelter et al., 2018). Additionally, this
safety filter is a Second-Order Cone Program (SOCP) (Boyd and Vandenberghe,
2004) for which an array of solvers exist including ECOS (Domahidi et al.,
2013).

Integrating Learning to Tune the Control Barrier Function: The parameter
selection process of TR-OP is particularly important, since the parameters
$\underline{a}$ and $\underline{b}$ guaranteed to exist by Theorem 4.1 are
worst-case approximations of the uncertainty generated using Lipschitz
constants. Such approximations often lead to undesired conservatism and may
render the system incapable of performing its goal (as seen in Figure 3).
Thus, as illustrated in Figure 1, we propose utilizing SA-LineCoSpar to
identify user-preferred parameters of TR-OP. This relaxes the worst-case over-
approximation to experimentally realize performant and safe behavior. This
design paradigm relies on the tunable construction of TR-OP, allowing us to
define the actions for SA-LineCoSpar to $\mathbf{a}=(\alpha,\varphi,a,b)$. We
note the construction of TR-OP assures that unsafe actions are not necessarily
catastrophic, as any $\alpha,\varphi,a,b>0$ endows the system with a degree of
robustness to disturbances and measurement error. This assurance allows us to
utilize a safety-aware approach as opposed to more conservative safety-
critical approach to learning.

##  5 Experimental Results

Figure 3: (Left) Actions sampled during simulation in 30 iterations with 3 new
actions in each iteration. The preferred action,
$\hat{\mathbf{a}}_{30}=(3,0.6,0.5,0.015)$, is shown in black and white. A
conservative action, $\mathbf{a}=(2,0.5,0.0651,0.485)$, is indicated by the
black circle, where $a$ and $b$ were determined by estimating the Lipschitz
coefficients present in the proof of Theorem 4.1. The conservative action
fails to progress whereas SA-LineCoSpar provides an action which successfully
navigates between obstacles. (Center) The minimum value of $h$ that occurred
in each iteration. Triangles, diamonds, and squares represent actions that are
sampled randomly, by PBL in simulation and on hardware in an indoor setting,
respectively. Colors correlate to iteration number. The lower bound
$-\gamma(\delta)$ for the expanded set $\mathcal{C}_{\delta}$ with $\delta=1$
is plotted. The preferred actions for simulation and hardware experiments are
circled. (Right) Seven additional iterations of 3 actions executed indoors.
The preferred action, $\hat{\mathbf{a}}_{37}^{*}=(4,0.6,0.4,0)$, successfully
traverses between the obstacles.

We applied the proposed design paradigm to a perception-based obstacle
avoidance task with a Unitree A1 quadrupedal robot (Figure 1) in simulation
and on hardware for both indoor and outdoor environments (see video: vid ).
The action space $A$ and hyperparameters of PBL are defined in Table 1. We
used the unicycle model as our simplified model (7) with the nominal
controller $\mathbf{k}_{\textrm{nom}}$:

$\displaystyle\underbrace{\begin{bmatrix}\dot{x}\\\ \dot{y}\\\
\dot{\psi}\end{bmatrix}}_{\dot{\mathbf{x}}}=\underbrace{\begin{bmatrix}0\\\
0\\\
0\end{bmatrix}}_{\mathbf{f}(\mathbf{x})}+\underbrace{\begin{bmatrix}\cos\psi&0\\\
\sin\psi&0\\\
0&1\end{bmatrix}}_{\mathbf{g}(\mathbf{x})}\left(\underbrace{\begin{bmatrix}v\\\
\omega\end{bmatrix}}_{\mathbf{v}}+\mathbf{d}(t)\right),$
$\displaystyle\mathbf{k}_{\textrm{nom}}(\mathbf{x})=\begin{bmatrix}K_{v}d_{\textrm{g}}+C\\\
-K_{\omega}(\sin\psi-(y_{\textrm{g}}-y)/d_{g})\end{bmatrix},$ (19)

where $(x,y)$ is the planar position of the robot, $\psi$ is the yaw angle,
$(x_{\textrm{g}},y_{\textrm{g}})$ is the goal position of the robot,
$d_{\textrm{g}}=\|(x_{\textrm{g}}-x,y_{\textrm{g}}-y)\|$ is the distance to
the goal, and $K_{v},K_{\omega}$, and $C$ are positive constants. Obstacle
avoidance is encoded via the 0-superlevel set of the function:

$h(\mathbf{x},\boldsymbol{\rho}_{i})=d_{\textrm{obs},i}-r_{\textrm{obs}}-\zeta\cos(\psi-\theta_{i}),$
(20)

where $\boldsymbol{\rho}_{i}=[x_{\textrm{obs},i},y_{\textrm{obs},i}]$ is the
location of the $i^{th}$ obstacle,
$d_{\textrm{obs},i}=\|(x_{\textrm{obs},i}-x,y_{\textrm{obs},i}-y)\|$ and
$\theta_{i}=\arctan((y_{\textrm{obs},i}-y)/(x_{\textrm{obs},i}-x))$ are the
distance and angle from the $i^{th}$ obstacle, $r_{\textrm{obs}}$ is the sum
of the radii of the obstacle and robot, and $\zeta>0$ determines the effect of
the heading angle on safety. The controller used to drive the system is the
TR-OP with the nominal controller $\mathbf{k}_{\textup{nom}}$ from (19). In
practice, infeasibilities of this safety filter were considered unsafe and the
inputs were saturated such that $v\in[-0.2,0.3]\,\textrm{m/s}$ and
$\omega\in[-0.4,0.4]\,\textrm{rad/s}$. The velocity command $\mathbf{v}$ is
computed at 20 Hz and error introduced by this sampling scheme is captured by
the tracking error $\mathbf{d}(t)$. Tracking of $\mathbf{v}$ is performed by
an inverse dynamics quadratic program (ID-QP) walking controller designed
using the concepts in Buchli et al. (2009), which realizes a stable walking
gait for (17) at 1 kHz.

Figure 4: The preferred action, $\hat{\mathbf{a}}_{40}^{*}=(5,0.1,0.4,0.02)$,
after simulation, indoor experiments, and 3 additional iterations of 3 actions
in an outdoor environment is shown alongside views from the onboard camera.

hyperparameter | value
---|---
$\lambda$ | $-0.5$
$\beta$ | 0

name | min. | max. | $\Delta$
---|---|---|---
$\alpha$ | 0.5 | 5 | 0.5
$\varphi$ | 0 | 1 | 0.1
$a$ | 0 | 1 | 0.1
$b$ | 0 | 0.05 | 0.005

Table 1: The safety-aware hyperparameters, and action space bounds (min. and
max.) with discretizations $\Delta$.

Simulation results: We simulated the quadruped executing the proposed
controller with parameters provided by SA-LineCoSpar. The resulting
trajectories and the position of the obstacles are shown in Figure 3. We ran
30 iterations, with 3 new actions sampled in each iteration ($s=3$), and
obtained user preferences and ordinal labels in between each set of actions.
To simulate perception error, the measurements of the obstacles were shifted
by $-0.1$ m in the $y$-direction. The parameters found with SA-LineCoSpar
allow the robot to navigate between obstacles. For comparison, a conservative
action is also shown, which is safe but fails to progress towards the goal.
SA-LineCoSpar eliminates this conservatism with only minor safety violations
and determines a parameter set which is both safe and performant.

Hardware results: After simulation, we continued learning on hardware
experiments in a laboratory setting for 7 additional iterations until the user
was satisfied with the experimental behavior. The robot and obstacle positions
were estimated using Intel RealSense T265 and D415 cameras to perform SLAM and
segmentation. Centroids of segmented clusters in the occupancy map were used
as the measured obstacle positions $\widehat{\boldsymbol{\rho}}_{i}$. The true
robot and obstacle positions were obtained for comparison using an OptiTrack
motion capture system. The results of these experiments can be seen in Figure
3. Afterwards, three additional iterations were conducted outdoors on grass
until again the user was satisfied with the experimental behavior. The
resulting best trajectory can be seen in Figure 4. The preferred action was
also tested on a variety of other obstacle arrangements to confirm its
generalizability. The performance of the final preferred action for these
obstacle configurations can be seen in the supplementary video (vid, ).

##  6 Conclusion

In this work we proposed a design paradigm for control systems in which the
robust safety requirements of a provably safe, but conservative controller are
relaxed, and controller parameters are instead chosen using a Preference-Based
Learning algorithm called SA-LineCoSpar. Using our algorithm, we were able to
learn a set of parameters that leads to user-preferred balance between safety
and robustness on a quadrupedal robot platform. Future work includes applying
this framework to other platforms such as bipedal robots, autonomous vehicles,
and assistive devices, and to more complicated environments like obstacles
with time-varying parameters.

## References

  * (1) Supplementary video, https://vimeo.com/649746749.
  * Alan et al. (2022) A. Alan, A. J. Taylor, C. R. He, G. Orosz, and A. D. Ames.  Safe controller synthesis with tunable input-to-state safe control barrier functions.  _Control Systems Letters_ , 6:908–913, 2022.
  * Ames et al. (2014) A. Ames, J. Grizzle, and P. Tabuada.  Control barrier function based quadratic programs with application to adaptive cruise control.  In _Conference on Decision & Control (CDC)_, pages 6271–6278. IEEE, 2014.
  * Ames et al. (2019) A. D. Ames, S. Coogan, M. Egerstedt, G. Notomista, K. Sreenath, and P. Tabuada.  Control barrier functions: Theory and applications.  In _European Control Conference (ECC)_ , pages 3420–3431. IEEE, 2019\.
  * Amodei et al. (2016) D. Amodei, C. Olah, J. Steinhardt, P. Christiano, J. Schulman, and D. Man.  Concrete problems in AI safety.  _arXiv preprint arXiv:1606.06565_ , 2016.
  * Berkenkamp et al. (2016) F. Berkenkamp, A. P. Schoellig, and A. Krause.  Safe controller optimization for quadrotors with Gaussian processes.  In _International Conference on Robotics and Automation (ICRA)_ , pages 491–496. IEEE, 2016.
  * Bıyık et al. (2020) E. Bıyık, N. Huynh, M. J. Kochenderfer, and D. Sadigh.  Active preference-based Gaussian process regression for reward learning.  _arXiv preprint arXiv:2005.02575_ , 2020.
  * Boyd and Vandenberghe (2004) S. Boyd and L. Vandenberghe.  _Convex optimization_.  Cambridge University Press, 2004.
  * Buchli et al. (2009) J. Buchli, M. Kalakrishnan, M. Mistry, P. Pastor, and S. Schaal.  Compliant quadruped locomotion over rough terrain.  In _2009 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , pages 814–820. IEEE, 2009.
  * Castañeda et al. (2020) F. Castañeda, J. J. Choi, B. Zhang, C. J. Tomlin, and K. Sreenath.  Gaussian process-based min-norm stabilizing controller for control-affine systems with uncertain input effects.  _arXiv preprint arXiv:2011.07183_ , 2020.
  * Chapelle and Li (2011) O. Chapelle and L. Li.  An empirical evaluation of Thompson sampling.  _Advances in Neural Information Processing Systems_ , 24:2249–2257, 2011.
  * Choi et al. (2021) J. J. Choi, D. Lee, K. Sreenath, C. J. Tomlin, and S. L. Herbert.  Robust control barrier-value functions for safety-critical control.  _arXiv preprint arXiv:2104.02808_ , 2021.
  * Chu and Ghahramani (2005) W. Chu and Z. Ghahramani.  Preference learning with Gaussian processes.  In _International Conference on Machine Learning (ICML)_ , pages 137–144, 2005.
  * Chu et al. (2005) W. Chu, Z. Ghahramani, and C. K. I. Williams.  Gaussian processes for ordinal regression.  _Journal of Machine Learning Research_ , 6(7), 2005.
  * Clark (2019) A. Clark.  Control barrier functions for complete and incomplete information stochastic systems.  In _American Control Conference (ACC)_ , pages 2928–2935. IEEE, 2019\.
  * Cosner et al. (2021) R. K. Cosner, A. W. Singletary, A. J. Taylor, T. G. Molnar, K. L. Bouman, and A. D. Ames.  Measurement-robust control barrier functions: Certainty in safety with uncertainty in state.  _arXiv preprint arXiv:2104.14030_ , 2021.
  * De Gemmis et al. (2009) M. De Gemmis, L. Iaquinta, P. Lops, C. Musto, F. Narducci, and G. Semeraro.  Preference learning in recommender systems.  _Preference Learning_ , 41:41–55, 2009.
  * Dean et al. (2020) S. Dean, A. J. Taylor, R. K. Cosner, B. Recht, and A. D. Ames.  Guaranteeing safety of learned perception modules via measurement-robust control barrier functions.  In _Conference on Robotics Learning (CoRL)_ , 2020.
  * Domahidi et al. (2013) A. Domahidi, E. Chu, and S. Boyd.  ECOS: An SOCP solver for embedded systems.  In _European Control Conference (ECC)_ , pages 3071–3076. IEEE, 2013\.
  * Glotfelter et al. (2018) P. Glotfelter, J. Cortés, and M. Egerstedt.  Boolean composability of constraints and control synthesis for multi-robot systems via nonsmooth control barrier functions.  In _Conference on Control Technology and Applications (CCTA)_ , pages 897–902. IEEE, 2018.
  * Gurriet et al. (2018) T. Gurriet, A. Singletary, J. Reher, L. Ciarletta, E. Feron, and A. Ames.  Towards a framework for realizable safety critical control through active set invariance.  In _International Conference on Cyber-Physical Systems (ICCPS)_ , pages 98–106. IEEE Press, 2018.
  * Hobbs et al. (2021) K. Hobbs, M. Mote, M. Abate, S. Coogan, and E. Feron.  Run time assurance for safety-critical systems: An introduction to safety filtering approaches for complex control systems.  _arXiv preprint arXiv:2110.03506_ , 2021.
  * Jain et al. (2015) A. Jain, S. Sharma, T. Joachims, and A. Saxena.  Learning preferences for manipulation tasks from online coactive feedback.  _The International Journal of Robotics Research_ , 34(10):1296–1313, 2015.
  * Jankovic (2018) M. Jankovic.  Robust control barrier functions for constrained stabilization of nonlinear systems.  _Automatica_ , 96:359–367, 2018.
  * Kirschner et al. (2019) J. Kirschner, M. Mutny, N. Hiller, R. Ischebeck, and A. Krause.  Adaptive and safe Bayesian optimization in high dimensions via one-dimensional subspaces.  In _International Conference on Machine Learning_ , pages 3429–3438. PMLR, 2019.
  * Kolathaya and Ames (2018) S. Kolathaya and A. D. Ames.  Input-to-state safety with control barrier functions.  _Control Systems Letters_ , 3(1):108–113, 2018\.
  * Li et al. (2021) K. Li, M. Tucker, E. Bıyık, E. Novoseller, J. W. Burdick, Y. Sui, D. Sadigh, Y. Yue, and A. D. Ames.  Roial: Region of interest active learning for characterizing exoskeleton gait preference landscapes.  In _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 3212–3218. IEEE, 2021.
  * Ma et al. (2017) W. Ma, S. Kolathaya, E. R. Ambrose, C. M. Hubicki, and A. D. Ames.  Bipedal robotic running with DURUS-2D: Bridging the gap between theory and experiment.  In _International Conference on Hybrid Systems: Computation and Control (HSCC)_ , pages 265–274. ACM, 2017.
  * Molnar et al. (2021) T. G. Molnar, R. K. Cosner, A. W. Singletary, W. Ubellacker, and A. D. Ames.  Model-free safety-critical control for robotic systems.  _arXiv preprint arXiv:2109.09047_ , 2021.
  * Raman et al. (2013) K. Raman, T. Joachims, P. Shivaswamy, and T. Schnabel.  Stable coactive learning via perturbation.  In _International conference on machine learning_ , pages 837–845. PMLR, 2013.
  * Sadigh et al. (2017) D. Sadigh, A. D. Dragan, S. Sastry, and S. A. Seshia.  Active preference-based learning of reward functions.  In _Robotics: Science and Systems_ , 2017.
  * Santoyo et al. (2019) C. Santoyo, M. Dutreix, and S. Coogan.  Verification and control for finite-time safety of stochastic systems via barrier functions.  In _2019 IEEE Conference on Control Technology and Applications (CCTA)_ , pages 712–717. IEEE, 2019.
  * Shivaswamy and Joachims (2012) P. Shivaswamy and T. Joachims.  Online structured prediction via coactive learning.  _arXiv preprint arXiv:1205.4213_ , 2012.
  * Singletary et al. (2021) A. Singletary, K. Klingebiel, J. R. Bourne, N. A. Browning, P. Tokumaru, and A. D. Ames.  Comparative analysis of control barrier functions and artificial potential fields for obstacle avoidance.  In _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2021.
  * Sui et al. (2015) Y. Sui, A. Gotovos, J. Burdick, and A. Krause.  Safe exploration for optimization with Gaussian processes.  In _International Conference on Machine Learning (ICML)_ , pages 997–1005, 2015.
  * Sui et al. (2018) Y. Sui, V. Zhuang, J. W. Burdick, and Y. Yue.  Stagewise safe Bayesian optimization with Gaussian processes.  In _International Conference on Machine Learning (ICML)_ , 2018.
  * Takano and Yamakita (2018) R. Takano and M. Yamakita.  Robust constrained stabilization control using control Lyapunov and control barrier function in the presence of measurement noises.  In _Conference on Control Technology and Applications (CCTA)_ , pages 300–305. IEEE, 2018.
  * Taylor and Ames (2020) A. J. Taylor and A. D. Ames.  Adaptive safety with control barrier functions.  In _American Control Conference (ACC)_ , pages 1399–1405. IEEE, 2020\.
  * Taylor et al. (2020) A. J. Taylor, A. Singletary, Y. Yue, and A. D. Ames.  Learning for safety-critical control with control barrier functions.  _Proceedings of Machine Learning Research (PMLR)_ , 120:708–717, 2020.
  * Tucker et al. (2020a) M. Tucker, M. Cheng, E. Novoseller, R. Cheng, Y. Yue, J. W. Burdick, and A. D. Ames.  Human preference-based learning for high-dimensional optimization of exoskeleton walking gaits.  In _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pages 3423–3430. IEEE, 2020a.
  * Tucker et al. (2020b) M. Tucker, E. Novoseller, C. Kann, Y. Sui, Y. Yue, J. W. Burdick, and A. D. Ames.  Preference-based learning for exoskeleton gait optimization.  In _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 2351–2357. IEEE, 2020b.
  * Tucker et al. (2021) M. Tucker, N. Csomay-Shanklin, W.-L. Ma, and A. D. Ames.  Preference-based learning for user-guided HZD gait generation on bipedal walking robots.  In _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 2804–2810. IEEE, 2021.
  * Wang et al. (2018) L. Wang, E. A. Theodorou, and M. Egerstedt.  Safe learning of quadrotor dynamics using barrier certificates.  In _International Conference on Robotics and Automation (ICRA)_ , pages 2460–2465. IEEE, 2018.
  * Yue et al. (2012) Y. Yue, J. Broder, R. Kleinberg, and T. Joachims.  The k-armed dueling bandits problem.  _Journal of Computer and System Sciences_ , 78(5):1538–1556, 2012.

##  Appendix A Proof of Theorem 4.1

In this appendix we present the proof of Theorem 4.1 which establishes the
input-to-state safety of $\mathcal{C}$ in the context of disturbances and
measurement uncertainties. Recall that Theorem 4.1 is given as:

###### Theorem

4.1 Given the set $\mathcal{C}$ defined in (9), suppose the functions
$L_{\mathbf{f}}h$, $L_{\mathbf{g}}h$, $\|L_{\mathbf{g}}h\|^{2}$, and
$\alpha\circ h$ are Lipschitz continuous on their domains, and assume that
$\|\widehat{\boldsymbol{\rho}}-\boldsymbol{\rho}\|\leq\epsilon$ for some
$\epsilon\in\mathbb{R}_{\geq 0}$. Then there exists
$\underline{a},\underline{b}\in\mathbb{R}_{\geq 0}$ such that if $h$
satisfies:

$\sup_{\mathbf{v}\in\mathbb{R}^{m}}L_{\mathbf{f}}h(\mathbf{x},\widehat{\boldsymbol{\rho}})+L_{\mathbf{g}}h(\mathbf{x},\widehat{\boldsymbol{\rho}})\mathbf{v}-\varphi\|L_{\mathbf{g}}h(\mathbf{x},\widehat{\boldsymbol{\rho}})\|^{2}-a-b\|\mathbf{v}\|>-\alpha(h(\mathbf{x},\widehat{\boldsymbol{\rho}})),$
(21)

for all $\mathbf{x}\in\mathbb{R}^{n}$ and some $a,b\in\mathbb{R}_{\geq 0}$
satisfying $a\geq\underline{a}$ and $b\geq\underline{b}$, then there exists a
controller $\mathbf{k}:\mathbb{R}^{n}\to\mathbb{R}^{m}$ such that (17) is ISSf
with respect to $\mathcal{C}$ with
$\gamma(\delta)=-\alpha^{-1}(-\delta^{2}/(4\varphi))$ where
$\alpha^{-1}\in\mathcal{K}_{\infty}^{\rm e}$.

######  Proof A.1.

First, we show that satisfying (21) for a particular set of
$\underline{a},\underline{b}$ implies satisfaction of (12). For this we
choose:

$\displaystyle\underline{a}=\epsilon(\mathcal{L}_{L_{\mathbf{f}}h}+\mathcal{L}_{\alpha\circ
h}+\mathcal{L}_{\varphi\|L_{\mathbf{g}}h\|^{2}}),$
$\displaystyle\underline{b}=\epsilon\mathcal{L}_{L_{\mathbf{g}}h},$ (22)

where $\mathcal{L}$ indicates the Lipschitz coefficient of the subscripted
function with respect to argument $\boldsymbol{\rho}$. Let us define the
function
$c:\mathbb{R}^{n}\times\mathbb{R}^{p}\times\mathbb{R}^{m}\to\mathbb{R}$ such
that:

$c(\mathbf{x},\boldsymbol{\rho},\mathbf{v})\triangleq
L_{\mathbf{f}}h(\mathbf{x},\boldsymbol{\rho})+L_{\mathbf{g}}h(\mathbf{x},\boldsymbol{\rho})\mathbf{v}-\varphi\|L_{\mathbf{g}}h(\mathbf{x},\boldsymbol{\rho})\|^{2}+\alpha(h(\mathbf{x},\boldsymbol{\rho})).$
(23)

Using this definition we have for all $\mathbf{v}\in\mathbb{R}^{m}$ that:

$\displaystyle c(\mathbf{x},\boldsymbol{\rho},\mathbf{v})$
$\displaystyle=c(\mathbf{x},\widehat{\boldsymbol{\rho}},\mathbf{v})+c(\mathbf{x},\boldsymbol{\rho},\mathbf{v})-c(\mathbf{x},\widehat{\boldsymbol{\rho}},\mathbf{v}),$
(24) $\displaystyle\geq
c(\mathbf{x},\widehat{\boldsymbol{\rho}},\mathbf{v})-\underbrace{\epsilon(\mathcal{L}_{L_{\mathbf{f}}h}+\mathcal{L}_{\alpha\circ
h}+\mathcal{L}_{\varphi\|L_{\mathbf{g}}h\|^{2}})}_{\underline{a}}-\underbrace{\epsilon\mathcal{L}_{L_{\mathbf{g}}h}}_{\underline{b}}\|\mathbf{v}\|$
(25) $\displaystyle\geq
c(\mathbf{x},\widehat{\boldsymbol{\rho}},\mathbf{v})-a-b\|\mathbf{v}\|.$ (26)

Above we added zero in (24) and used the Lipschitz coefficients and the worst
case uncertainty $\epsilon$ to achieve the bound in (25). Since
$\sup_{\mathbf{v}\in\mathbb{R}^{m}}c(\mathbf{x},\widehat{\boldsymbol{\rho}},\mathbf{v})-a-b\|\mathbf{v}\|>0$
holds based on (21), inequality (26) implies that (12) holds for the true
parameters, $\boldsymbol{\rho}$. Since (12) holds, the conditions of Theorem
3.1 point 2 are satisfied and thus $\mathcal{C}$ is ISSf with
$\gamma(\delta)=-\alpha^{-1}(-\delta^{2}/(4\varphi))$ where
$\alpha^{-1}\in\mathcal{K}_{\infty}^{\rm e}$.