# Posterior Regularized Bayesian Neural Network Incorporating Soft and Hard
Knowledge Constraints

Jiayu Huang  Yutian Pang  Yongming Liu  Hao Yan<EMAIL_ADDRESS>
###### Abstract

Neural Networks (NNs) have been widely used in supervised learning due to
their ability to model complex nonlinear patterns, often presented in high-
dimensional data such as images and text. However, traditional NNs often lack
the ability for uncertainty quantification. Bayesian NNs (BNNS) could help
measure the uncertainty by considering the distributions of the NN model
parameters. Besides, domain knowledge is commonly available and could improve
the performance of BNNs if it can be appropriately incorporated. In this work,
we propose a novel Posterior-Regularized Bayesian Neural Network (PR-BNN)
model by incorporating different types of knowledge constraints, such as the
soft and hard constraints, as a posterior regularization term. Furthermore, we
propose to combine the augmented Lagrangian method and the existing BNN
solvers for efficient inference. The experiments in simulation and two case
studies about aviation landing prediction and solar energy output prediction
have shown the knowledge constraints and the performance improvement of the
proposed model over traditional BNNs without the constraints.

###### keywords:

Posterior regularization , Bayesian Neural Network , Knowledge constraint ,
Augmented Lagrangian method , Soft and hard constraint

We propose a PR-BNN model to incorporate various knowledge constraints into
BNNs.

Both deterministic and probabilistic views of soft and hard constraints are
defined.

Efficient algorithms are developed to solve models with soft or hard
constraints.

Augmented Lagrangian method is applied to solve PR-BNN with multiple hard
constraints.

In implementations, various knowledge applied in PR-BNN shows advantages over
BNNs.

##  1 Introduction

Neural Networks (NNs) have gained great popularity due to their better
accuracy in modeling datasets with inherent complex problem structures [8].
However, for critical tasks requiring accurate risk quantification in complex
engineered systems or applications, where system safety is important, the
deterministic prediction model provided by traditional NNs may not be
sufficient. Recently, combining the power of Bayesian inference and the
expressive power of the deep neural networks, Bayesian Neural Networks (BNNs)
have been proposed to combine high representation power with the quantifiable
uncertainty estimates [26]. More specifically, traditional NNs cannot quantify
epistemic uncertainty, which refers to the systematic uncertainty due to the
uncertainty of the models. BNNs put prior distributions on the Neural Network
weights. Furthermore, the prior distribution of the NN parameters also helps
reduce the bias of the model. Finally, the posterior distribution of the
predicted model help quantify the distribution of the output, which could lead
to better decision-making in the downstream tasks [18].

However, despite the reduction of the bias and good performance in many
prediction tasks, purely data-driven BNNs often have the following
limitations: 1) Due to the versatility of the Bayesian Neural Network models,
the results of the BNN models can still suffer from overfitting due to the
complexity of the BNN models, even though the prior regularization is applied.
This could be due to the lack of training data or lack of data in a particular
testing regime, which may result in poor generalization power. 2) Results from
BNN models trained purely based on data may not comply with physical laws or
engineering constraints, leading to a large bias or uninterpretable results.

In engineering systems, complex knowledge constraints may often exist. Common
knowledge constraints include the range constraint, smoothness constraint,
monotonic relationship between variables, the parametric relationship between
variables, etc. Incorporating such general knowledge constraints into the BNNs
may be crucial for complex engineering systems to reduce bias and improve
interpretability, especially in regions with limited training data.

In the line of knowledge-constrained learning research, there are two major
ways to incorporate constraints into machine learning models, parameter
constraints and output constraints. Parameter constraints put the constraints
on the model functions to regularize the model training. Parameter constraints
have been widely used in linear models, given that these model parameters are
often correlated to some physical meanings in engineering systems.
Furthermore, adding such parameter constraints for high-dimensional models is
especially important given that it can often avoid overfitting [20, 11]. These
constraints on the model parameters include the smoothness constraints, range
constraints [5], sparsity constraints [30], and low-rank constraints [27].
Much of the literature focuses on imposing informative prior distribution to
empower learning tasks with additional knowledge and constraint in Bayesian
inference [7]. Qualitative knowledge regarding the relationships among
variables can be characterized as signs or order constraints, which can be
applied via prior construction into Bayesian Network to exhibit more sensible
estimates [9].

However, in complex models such as deep learning models, model parameters
often do not have any physical meanings, and specifying a meaningful prior
distribution on the model parameters is often impossible. In the literature,
another way to integrate knowledge constraints into the machine learning model
is the output-based constraint. Instead of adding the constraint on the model
parameters, the constraints can be added directly to the output prediction of
the machine learning models as the regularization. In literature, Regularized
Expected Mean Log-likelihood Prior, Normal Wishart distribution has been
proposed as prior to combining some biological knowledge [9]. However,
modifying prior for general knowledge constraints is not straightforward for
many different applications due to the black-box models. This is especially
true when we have multiple knowledge constraints from different types or data-
dependent complex knowledge constraints, and the prior distributions may
become harder to define or even not mathematically exist [32]. Recently, OC-
BNN (Output-Constrained BNN) has been proposed in [29]. However, OC-BNN work
has some limitations on constraint type and prior distribution modification.
For more details about OC-BNN, please refer to Section 2.3.

Incorporating domain knowledge via Posterior Regularization (PR) has been
proposed into traditional machine learning models [13]. PR is defined as a
regularized maximum likelihood function, which has the ability to directly
constrain the posterior distribution. For more details about the PR framework,
please refer to Section 2.

In this paper, to incorporate more general knowledge constraints in BNNs, we
combine the PR framework with the regularized variational inference to create
the generic framework named Posterior Regularized BNN (PR-BNN). Our
contributions are listed below:

  * 1.

Based on the underlying properties of different knowledge constraints, we
define deterministic and probabilistic scenarios for both soft and hard
knowledge constraints.

  * 2.

The proposed PR-BNN applies posterior regularization to incorporate multiple
general constraints in the forms of score functions. Both soft and hard
constraints can be considered.

  * 3.

The algorithm, with the basis of different Bayesian inference methods, is
developed to solve the soft constraint problem under PR-BNN. The efficient
algorithm based on the Augmented Lagrangian form is proposed to solve the hard
constraint problem.

The rest of this paper is organized as follows. Section 2 is a review of
related methods, Section 3 describes the formulation of the Proposed PR-BNN
framework, and section 4 and 5 present two simulation experiments and two
implementations on real-world data, respectively. Section 6 concludes the
paper briefly.

##  2 Review of Related Methods

In this section, we would like to briefly introduce some methodologies related
to the proposed method.

###  2.1 Review of Bayesian Neural Network (BNNs)

In this section, the basic BNNs algorithm will be introduced. The notations
are listed in Table 1.

Table 1: Notations in BNNs Notations |  Meaning
---|---
$\mathbf{x}$ |  The training inputs
$\mathbf{y}$ |  The prediction outputs
$\mathcal{D}$ |  All training variables including input and output variables
$\mathbf{w}$ |  The BNNs parameters such as weights
$\theta$ |  The parameter vectors for distributions of $\mathbf{w}$
$x^{*}$ |  The new inputs
$y^{*}$ |  The predictions corresponding to new input $x^{*}$
$q$ |  The approximated variational distribution

Unlike the traditional NNs, model parameters in BNNs $\mathbf{w}$ are
distributions, which can be represented as $\mathbf{w}\sim t(\theta)$, where
$\theta$ is a deterministic parameter set for the distribution of
$\mathbf{w}$. Uncertainty will be invoked by the distribution of $\mathbf{w}$,
which is called epistemic uncertainty [16] and the over-fitting issue of NNs
can be relieved as well.

Given $n$ input data vectors $\mathbf{x}$ and outputs $\mathbf{y}$, we try to
find the posterior distributions of parameters $\mathbf{w}$, that is,
$p(\mathbf{w}|\mathbf{x},\mathbf{y})$. Then try to generate the true
distribution $p(\mathbf{y}|\mathbf{x})$ via estimated
$p(\mathbf{w}|\mathbf{x},\mathbf{y})$.

$\overbrace{p(\mathbf{w}|\mathbf{x},\mathbf{y})}^{\text{Posterior}}=\frac{\overbrace{p(\mathbf{y}|\mathbf{x},\mathbf{w})}^{\text{Likelihood}}\overbrace{p(\mathbf{w})}^{\text{Prior}}}{\underbrace{p(\mathbf{y}|\mathbf{x})}_{\text{Evidence}}}$
(1)

In practice, distribution space of $\mathbf{w}$ is often defined as a prior,
denoted as $p(\mathbf{w})$. $p(\mathbf{y}|\mathbf{x},\mathbf{w})$ is the
likelihood function based on training dataset. The normaliser in Equation 1,
$p(\mathbf{y}|\mathbf{x})$, is called evidence [12]. Thus, the posterior of
$\mathbf{w}$, $p(\mathbf{w}|\mathbf{x},\mathbf{y})$, is evaluated with
Equation 1.

Unlike the deterministic model, predictive results of BNNs depend on the
estimated posterior distribution of $\mathbf{w}$, instead of deterministic
$\mathbf{w}$. Predictive distribution will be

$p(y^{*}|x^{*})=\int
p(y^{*}|x^{*},\mathbf{w})p(\mathbf{w}|\mathbf{x},\mathbf{y})d(\mathbf{w})$ (2)

As shown in (2), the predicted output $y^{*}$ is the integration over
$\mathbf{w}$, which can be also regarded as an expectation formulated as
$\mathbb{E}_{p(\mathbf{w}|\mathbf{x},\mathbf{y})}[p(y^{*}|x^{*},\mathbf{w})]$.
Unfortunately, this expectation is intractable for any practical case in NNs
[4].

Variational inference [3, 31] has been widely used to solve the intractability
problem in posterior inference of BNNs. Generally, a simpler approximated
distribution $q(\mathbf{w}|\theta)$ is used to approximate
$p(\mathbf{w}|X,Y)$. Kullback-Leibler (KL) divergence [19] is adopted here to
measure the distance between distribution $q(\cdot)$ and $p(\cdot)$, which can
be denoted as $KL(q(\cdot)||p(\cdot))$. By minimizing KL divergence over
$\theta$, we are able to obtain optimal variational distribution $q(\cdot)$.

ELBO, a useful lower bound [17] on the negative log-likelihood function, has
been introduced to get the objective function in BNNs as shown in (3) in the
training process.

$\displaystyle\widehat{-ELBO}:=\underbrace{KL(q_{\theta}(\mathbf{w})||p(\mathbf{w}))}_{\text{KL
divergence}}\underbrace{-log(p(D|\mathbf{w}))}_{\text{Negative Log-Likelihood
(NLL)}}$ (3)

In the prediction, Monte Carlo (MC) sampling [15] can be applied to get the
mean prediction of output $y$, and the uncertainty can be estimated as in
Equation (4). $N$ is the number of MC samplings performed during prediction.

$\displaystyle p(y^{*}|x^{*})$
$\displaystyle=\mathbb{E}_{p(\mathbf{w}|\mathbf{x},\mathbf{y})}[p(y^{*}|x^{*},\mathbf{w})]$
(4)
$\displaystyle\approx\frac{1}{N}\sum^{N}_{i=1}p(y^{*}|x^{*},\hat{\mathbf{w}_{i}})$

Gaining from the NNs structure, BNNs can provide reliable prediction and
uncertainty quantification even for the high-dimensional (HD) data [6].
Besides, BNNs can quantify complex uncertainty, such as spatial/temporal
variations in HD data, via a rigorous framework [24].

###  2.2 Review of Posterior Regularization (PR)

In literature, the goal of posterior regularized constraint is to restrict the
space of mode posteriors as a way to guide the NN towards the behavior
governed by the engineering domain knowledge. Borrowing from the PR framework,
we define a general constraint function $f(\mathbf{x,y})$ on output
$\mathbf{y}$. It is worth noting that $f(\cdot)$ can be defined on both input
$\mathbf{x}$ and output $\mathbf{y}$, therefore, i.e.,
$f(\mathbf{x},\mathbf{y})$.

The methods can be represented as follows:

$\text{ PR:
}\max_{\Theta}\mathcal{L}(\Theta;\mathcal{D})+\mathbf{w}(p(\mathbf{y}\mid\mathcal{D},\Theta)),$
(5)

where $\mathcal{L}(\Theta;\mathcal{D})$ is the marginal likelihood of
$\mathcal{D}$, and $\mathbf{w}(\cdot)$ is a regularization function of the
posterior over latent or output variables $\mathbf{Y}$. The posterior here is
a posterior for a different learning model instead of the Bayesian posterior.
Here, the PR is used to constrain the output of non-Bayesian methods. To
generalize PR for the Bayesian method, RegBays [32] is proposed by using the
variational Bayesian inference form, as shown in Equation (6).

RegBayes:
$\displaystyle\inf_{q(\mathbf{w})}\mathrm{KL}(q(\mathbf{w})\|p(\mathbf{w}|\mathcal{D}))+U(\boldsymbol{\xi})$
(6) $\displaystyle=$
$\displaystyle\inf_{q(\mathbf{w}),\boldsymbol{\xi}}\mathrm{KL}\left(q(\mathbf{w})\|P(\mathbf{w})\right)-\int_{x}\log
P(D|\mathbf{w}))q(\mathbf{w})d\mathbf{w}+U(\boldsymbol{\xi})$
$\displaystyle\text{s.t.
}q(\mathbf{w})\in\mathcal{P}_{\text{post}}(\boldsymbol{\xi}),$

where $KL(q(\mathbf{w})\|p(\mathbf{w}\mid\mathcal{D}))$ is the KL-divergence
between a desired post-data posterior $q(\mathbf{w})$ over model weights
$\mathbf{w}$ and the true posterior distribution
$p(\mathbf{w}\mid\mathcal{D})$. The regularization $\mathbf{U}$ is a function
of slack variable $\xi$ (i.e., $\mathbf{U}(\xi)=\|\xi\|_{\beta}$ or
$\mathbf{U}(\xi)=\xi$). This function provides a flexible way to incorporate
additional information, such as domain knowledge. In this work, we will extend
RegBayes into a Bayesian deep learning framework to deal with complex data and
various types of knowledge constraints and solve the constrained BNNs.

###  2.3 Review of Output-constrained Bayesian Neural Network (OC-BNN)

Recently, OC-BNN has been proposed to integrate the domain knowledge to the
output variable $\mathbf{y}$ of a Bayesian Neural network. The
interpretability and generality of the output constraints are not hard to see
when we compare them with parameter constraints. More specifically, constraint
function $f(\mathbf{y})$ or $f(\mathbf{x},\mathbf{y})$ can be easily defined
in the knowledge about the output values or relationship between input and
output is given [29] .

Positive (negative) constraints $\mathcal{C}^{+(-)}$ are defined as satisfying
that iff $\forall x\in\mathcal{C}_{x}$, the output $y\in(\notin){C}_{y}(x)$.
For simplicity, $\mathcal{C}$ represents these two types of constraints in the
rest of this section.

Conditional Output-Constrained Prior $(\boldsymbol{COCP})$ is proposed by [29]
to incorporate output constraints into prior via Equation (7):

$p_{\mathcal{C}}(\mathbf{w})=\frac{1}{Z}p(\mathbf{w})\exp(f(\mathbf{y}))$ (7)

For multiple constraints, more exponential terms can be multiplied. The
objective function in this method is still variational inference for the BNNs
as described in Equation (3). Through modification on prior, the knowledge is
incorporated.

However, the prior distribution constant is often not easy to calculate or
even not well defined in many examples. It limits the types of constraints
that could be incorporated into the BNNs. For example, the constraint types
considered in OC-BNN are mainly about whether the output or the function of
output should be in or out of a pre-specified region. Finally, OC-BNN does not
have ways to decide the relative weights of different types of constraints,
and it doesn’t investigate how to guarantee hard knowledge constraints without
the use of rejective sampling.

In our work, we propose a novel knowledge-constrained BNNs framework named
Posterior Regularized BNNs (PR-BNNs). We define a more general constraint
function format and derive both soft and hard constraints and their
corresponding algorithms. PR has been applied to form the objective function
to simplify the problem with constraints in the BNN setting. For hard
constraints, an Augmented Lagrangian form of the posterior regularized
objective function has been derived, and the relative strength of the
constraint function can be automatically updated. We discuss the detailed
derivation in Section 3.

The main differences between OC-BNN work and our work include: First, in OC-
BNN work, the constraint type considered is mainly about conditional in or out
of region constraints, while in our work, we include more general functions to
represent constraints. Also, except for the soft constraints, we define
deterministic and probabilistic hard constraints to match the strictness of
knowledge. Second, different from OC-BNN, we show that incorporating domain
knowledge via posterior regularization, especially in the hard constraint, can
lead to a flexible framework that automatically updates the dual variable for
each piece of knowledge. This is especially important if we have multiple
constraints, given the proposed method can automatically decide the relative
importance of each constraint, thereby allowing for robust incorporation.
Please see [32] paper for more detailed discussions on the advantage of the
Posterior Regularization compared to the conditional output regularized prior
in the context of traditional machine learning methods.

##  3 Formulation

In this section, we will start with the introduction of the knowledge-
constraint function in this work, including the definition, types, and
examples for different types of knowledge constraints. Then we propose the
Posterior-Regularized Bayesian Neural Network (PR-BNN) framework for both soft
and hard constraints and efficient algorithms.

###  3.1 Definitions of Hard and Soft Knowledge Constraint via PR-BNN

In this section, we will define two types of constraints that can be used in
the PR-BNN, namely, the soft constraint and the hard constraint.

######  Definition 1.

Knowledge Function $f(\mathbf{x},\mathbf{y})$ is a score function of
$\mathbf{x}\text{ and }\mathbf{y}$, which is typically defined as the higher,
the better (i.e., satisfy the constraint).

######  Remark 1.

In this paper, we will use the idea of max-margin learning to define the
knowledge function $f(\mathbf{x,y})$ by
$f(\mathbf{x,y})=\min(0,s(\mathbf{x,y})+m)$, which $s(\mathbf{x,y})$ is a
score function. By maximizing the knowledge function $f(\mathbf{x,y})$, it
will push the knowledge function $f(\mathbf{x,y})$ to $0$ or the score
function $s(\mathbf{x,y})$ to be larger than a certain threshold $-m$ for a
particular $\mathbf{y}$ given input $\mathbf{x}$ , where m is a very small
positive value.

Furthermore, we assume that the predicted output $\hat{\mathbf{y}}$ are
generated from an NN model $\hat{\mathbf{y}}=g(\mathbf{x},\mathbf{w})$. In
BNN, it is commonly assumed that $\mathbf{w}$ is a random variable, which
follows a distribution $q(\mathbf{w})$. Therefore, we assume that the
randomness of $\mathbf{y}$ only comes from the uncertainty of $\mathbf{w}$. In
another word, the set of constraint $f$, $Q=\\{E_{q}[f(\hat{\mathbf{y}})]\geq
b\\}$ is equivalent on the constraint on $\mathbf{w}$ as
$Q_{\mathbf{w}}=\\{E_{q(\mathbf{w})}[f(g(\mathbf{x},\mathbf{w}))]\geq b\\}$.

######  Definition 2.

Hard Constraint: A posterior distribution $q(\mathbf{w})$ is said to satisfy
the hard constraint when $q(\mathbf{w})\in
Q_{\mathbf{w}}=\\{E_{q(\mathbf{w})}[f(\mathbf{x,y})]\geq b\\}$, where b is a
pre-specified constant. Both deterministic and probabilistic hard constraints
can be included by the above formula.

######  Proposition 1.

Deterministic Hard Constraint: If we assume $-c\leq f(\mathbf{x,y})\leq 0$
where $c\geq 0$, BNN $g(\cdot)$ is a bounded smooth function, and $b=0$,
$Q_{\mathbf{w}}=\\{E_{q(\mathbf{w})}[f(\mathbf{x,y})]\geq 0\\}$ since
$f(\mathbf{x,y})\leq 0$, we are able to get
$f(g(\mathbf{x},\mathbf{w}))=0,\quad\forall\mathbf{w}\in q(\mathbf{w})$.

This implies that for any posterior output $\mathbf{y}$ generated from the
BNN, the constraint should always be satisfied. In this case, a hard
constraint can be defined on the epistemic uncertainty, given all $w$ should
make $g(\mathbf{x},\mathbf{w})$ satisfy this constraint.

######  Definition 3.

Probabilistic Hard Constraint ($\epsilon-$satisfied constraint): We define
$f(\mathbf{x,y})=-I_{(\mathbf{y}\notin C)}$, which $\mathbf{y}\notin C$
implies that $y$ doesn’t satisfy the constraint. Therefore, the posterior set
with slack variable $\epsilon$ can be interpreted as the $\epsilon-$satisfied
constraint as $-P(\mathbf{y}\notin C)+\epsilon=E_{q}[f(\mathbf{x,y})]+\xi\geq
0$, where $\epsilon\in[0,1]$ is the strictness of the constraint. If
$\epsilon>0$, the constraint implies that $P(\mathbf{y}\notin C)\leq\epsilon.$

######  Definition 4.

Soft Constraint: A posterior distribution $q(\mathbf{w})$ is said to satisfy
the soft constraint defined by the knowledge function $f(\mathbf{x,y})$, when
$q(\mathbf{w})\in Q_{\mathbf{w}}=\\{E_{q(\mathbf{w})}[f(\mathbf{x,y})]+\xi\geq
0\\}$, where $\xi$ is a slack variable, instead of a known constant.

If we assume that $\xi$ is known, the method becomes the hard constraint
version where $b=-\xi$. However, the soft-constrained version typically
assumes that $\xi$ is not known, inspired by the maximum margin learning, we
would like to include the slack variable into the loss function.

Examples of $f(\mathbf{x,y})$ formulation for different types of knowledge
will be presented in Section 3.1.1. Furthermore, we will discuss the algorithm
for solving the soft and hard constraints in Section 3.2.

####  3.1.1 Examples of Knowledge Constraint Function

In this section, we would like to give some examples of the knowledge
constraint function on the output variable ($\mathbf{y}$). In order to relax
the constraint format, we assume a score function, which means that if the
predicted outputs accordance with the knowledge, then the value of
$f(\mathbf{x,y})$ would be high, and vice versa.

The types of constraints that can be included in this framework are rather
flexible and applicable to different fields. Here, some generic formats of
constraints will be listed, and how the design of the constraint functions
will be discussed in detail.

  * 1.

Conditional value constraint: Certain transformation functions $r(\cdot)$ of
the output values $\mathbf{y}$ should be close to a pre-specified value $C$ as
certain independent variables are within an input range
$[\mathbf{x}_{a},\mathbf{x}_{b}]$. Then the score function can be defined as
$f(\mathbf{x,y})=\min(-|r(\mathbf{x,y})-C|+m,0)$.

  * 2.

Boundary constraint: Certain transformation functions $r(\mathbf{x,y})$ should
satisfy the bound limit, including the upper and lower bounds (e.g.,
$[a(\mathbf{x}),b(\mathbf{x})]$). These bounds are not necessarily the same
value over a certain range $\mathbf{x}$. Instead, it can be a function of the
independent variables. If we have the constraint that $a(\mathbf{x})\leq
r(\mathbf{x,y})\leq b(\mathbf{x})$, the constraint function can be defined as
follows:
$f(\mathbf{x,y})=\min(0,-|r(\mathbf{x,y})-\frac{b(\mathbf{x})+a(\mathbf{x})}{2}|)$.

  * 3.

Monotonic constraint: Certain transformation learning $r(\mathbf{x,y})$ on the
output $\mathbf{y}$ will result in non-increasing or non-decreasing function.
In the time series problem, we can also define some special monotone
constraints which allow some small vibration but show some overall trend. For
non-decreasing constraint, $f(\mathbf{x,y})=\min(0,r(\mathbf{x,y})^{\prime})$.

  * 4.

Convex and Concave constraint: Certain transformation learning
$r(\mathbf{x,y})$ on the output $\mathbf{y}$ will result in convex and concave
functions. In some fields, such as physics, it is common to see some laws
including second derivative terms. In that case, we can include this law as a
knowledge constraint in BNNs by using the second derivative directly, as
follows: $f(\mathbf{x,y})=\min(0,r(\mathbf{x,y})^{\prime\prime})$.

######  Remark 2.

Multiple constraints: To incorporate multiple constraints, a linear
combination is a simple way to combine several available and useful
constraints with different formats and complexity, as follows:

$\mathbf{f}(\mathbf{x},\mathbf{y})=\Sigma_{i}(s_{i}f_{i}(\mathbf{x},\mathbf{y})),$

where $s_{i}$ is the weight of $i$-th knowledge constraint.

###  3.2 Posterior Regularized Bayesian Neural Networks (PR-BNNs)

In this section, we demonstrate core optimization problem solved in proposed
PR-BNNs and the corresponding algorithm to solve soft and hard version PR-
BNNs, respectively.

####  3.2.1 Proposed PR-BNNs

In this work, we focus on the output constraint given any valid input
$\mathbf{x}$, whose conditional expectation of the knowledge constraint
functions $f(\mathbf{x},\mathbf{y})$, is denoted as
$E[f(\mathbf{x},\mathbf{y})\mid\mathcal{D},\mathbf{w}]$, over the predictive
output $\mathbf{y}$ and input set $\mathbf{x}$ using dataset $\mathcal{D}$ and
trained weights $\mathbf{w}$. For the hard constraint version, the expectation
constraint over the estimated posterior distribution $q(\mathbf{w})$ should be
greater than $b$ to satisfy the knowledge constraint:
$E_{q(\mathbf{w})}E[f(\mathbf{x},\mathbf{y})|\mathcal{D},\mathbf{w}]\geq b$.
For multiple constraints, we can rewrite constraints into vector form as:
$E_{q(\mathbf{w})}E[\mathbf{f}(\mathbf{x},\mathbf{y})|\mathcal{D},\mathbf{w}]\geq\mathbf{b}$.

##### Hard constrained version

In other words, by applying this constraint into the variational form of the
Bayesian NN, the optimization problem in Equation (11) is derived. Bayesian
estimation solves with the constraints applied on the output and inputs of the
models as follows:

$\displaystyle\min_{q(\mathbf{w}),\mathbf{\xi}}KL[q(\mathbf{w})\|P(\mathbf{w})]-E_{q(\mathbf{w})}\log
P(\mathbf{y}|\mathcal{D},\mathbf{w})$ (8) $\displaystyle s.t.\quad
q(\mathbf{w})\in\mathcal{Q},\mathcal{Q}=\\{q(\mathbf{w}):E_{q(\mathbf{w})}E[\mathbf{f}(\mathbf{x},\mathbf{y})|\mathcal{D},\mathbf{w}]\geq\mathbf{b}\\},$

where $\mathcal{Q}$ denotes the regions, where the expected constraint
function involving the expected output is bounded by $b$.

##### Soft constrained version

In the soft-constrained version, it is often useful to allow some violations
of the constraints denoted by the slack variable $\xi$. Inspired by the
maximum margin learning, we would like to control the amount of violation
$\xi$ by a slack penalty $U(\xi)$.

To generalize the constraints, we would like to introduce the slack variable
$\xi$, where the constrained set $\mathcal{Q}_{\xi}$. To construct the
feasible non-empty constraints $\mathcal{Q}_{\xi}$ on the posterior
distribution estimation $q(\mathbf{w})$, we add a slack penalty $U(\xi)$ to
the regularized variational form of the Bayesian form:

$\displaystyle\min_{q(\mathbf{w}),\mathbf{\xi}}KL[q(\mathbf{w})\|P(\mathbf{w})]-E_{q(\mathbf{w})}\log
P(\mathbf{y}|\mathcal{D},\mathbf{w})+\lambda U(\xi)$ (9) $\displaystyle
s.t.\quad q(\mathbf{w})\in
Q_{\xi},Q_{\xi}=\\{q(\mathbf{w}):E_{q(\mathbf{w})}E[\mathbf{f}(\mathbf{x},\mathbf{y})|\mathcal{D},\mathbf{w}]\geq-\xi\\},$

where $\mathbf{w}$ is the weights of BNNs; $P(\mathbf{w}|D)$ is the true
posterior distribution of $\mathbf{w}$, $q(\mathbf{w})$ is used to estimate
the true posterior distribution $P(\cdot)$ due to the intractability of it;

For example, if $U(\xi)=\sum_{i}\lambda_{i}\xi_{i}$, Equation (9) becomes the
regularization form

$\min_{q(\mathbf{w})}KL[q(\mathbf{w})\|P(\mathbf{w})]-E_{q(\mathbf{w})}\log
P(\mathbf{y}|\mathcal{D},\mathbf{w})-\lambda
E_{q(\mathbf{w})}E[\mathbf{f}(\mathbf{x},\mathbf{y})|\mathcal{D},\mathbf{w}]$
(10)

In the following section, we will assume that $U(\xi)=\xi$ in the algorithm
development.

####  3.2.2 Optimization Algorithms to Solve PR-BNNs

In this section, we will derive the optimization algorithms for the PR-BNN. In
this paper, we focus on algorithms solving both the hard-constrained version
in Equation (8) and the soft-constrained version in Equation (9).

We assume the weights $\mathbf{w}$ in BNNs are distributions with parameter
$\theta$, so we need to search for optimal $\theta$ to solve the minimization
problem in Equation (10).

##### Soft constrained version

For the soft constrained version, we will use Equation (10) as an example, and
it is easy to derive that

$\displaystyle\arg\min_{\theta}KL[q(\mathbf{w|\theta})\|P(\mathbf{w})]-E_{q(\mathbf{w|\theta})}\log
P(\mathbf{y}|\mathcal{D},\mathbf{w})-\lambda
E_{q(\mathbf{w|\theta})}E[\mathbf{f}(\mathbf{x},\mathbf{y})|\mathcal{D},\mathbf{w}]$
(11) $\displaystyle=$
$\displaystyle\arg\min_{\theta}E_{q(\mathbf{w}|\theta)}\mathcal{F}(\mathbf{w},\theta),$

where $\mathcal{F}(\mathbf{w},\theta)=\log(q(\mathbf{w}|\theta))-\log
P(\mathbf{w})-\log P(\mathbf{y}|\mathcal{D},\mathbf{w})-\lambda
E[\mathbf{f}(\mathbf{x},\mathbf{y})|\mathcal{D},\mathbf{w}]$.

##### Hard constrained version

For the hard constraint version, we will first introduce the slack variable
$z_{i}$ of knowledge constraint $i$ in order to convert the original
inequality constraint into equality constraint. The optimization problem is
shown in Equation (12)

$\displaystyle\arg\min_{\theta}\quad
KL[q(\mathbf{w}|\theta)\|P(\mathbf{w}|D)],\quad$ (12) $\displaystyle s.t.\quad
E_{q(\mathbf{w}|\theta)}f_{i}(\mathbf{x,y})=z_{i},z_{i}\geq 0,i=1,\cdots,I$

The Augmented Lagrangian form can be derived as

$\displaystyle\min_{z,\theta}\quad KL[q(\mathbf{w}|\theta)\|P(\mathbf{w}|D)]-$
$\displaystyle\sum_{i}s_{i}(E_{q(\mathbf{w}|\theta)}f_{i}(\mathbf{x},\mathbf{y})-z_{i})+\sum_{i}\frac{1}{2}\rho_{i}\|E_{q(\mathbf{w}|\theta)}f_{i}(\mathbf{x},\mathbf{y})-z_{i}\|^{2},$
(13) $\displaystyle s.t.\quad z\geq 0.$

######  Proposition 2.

$z_{i}$ in Equation (13) can be solved as

$z_{i}=\max\\{E_{q(\mathbf{w}|\theta)}f_{i}(\mathbf{x},\mathbf{y})-\frac{s_{i}}{\rho_{i}},0\\}$
(14)

Finally, by plugging in Equation (14) to Equation (13), we can derive the
following proposition.

######  Proposition 3.

Solving $\theta$ in Equation (13) is equivalent to solve

$\min_{\theta}\quad
KL[q(\mathbf{w}|\theta)\|P(\mathbf{w}|D)]+\sum_{i}\phi(E_{q(\mathbf{w}|\theta)}f_{i}(\mathbf{w}),s_{i},\rho_{i}),$
(15)

where

$\phi(E_{q(\mathbf{w}|\theta)}f_{i}(\mathbf{w}),s_{i},\rho_{i})=\begin{cases}-s_{i}E_{q(\mathbf{w}|\theta)}f_{i}(\mathbf{w})+\frac{1}{2}\rho_{i}\left(E_{q(\mathbf{w}|\theta)}f_{i}(\mathbf{w})\right)^{2}&E_{q(\mathbf{w}|\theta)}f_{i}(\mathbf{w})\leq
s_{i}/\rho\\\
-\frac{1}{2}s_{i}^{2}/\rho_{i}&E_{q(\mathbf{w}|\theta)}f_{i}(\mathbf{w})>s_{i}/\rho\end{cases}.$

The proof is similar to [22]. Finally, the dual variable can be updated as
$s_{i}=\max\\{0,s_{i}-\rho_{i}E_{q(\mathbf{w}|\theta)}f_{i}(\mathbf{x},\mathbf{y})\\}$
and $\rho_{i}$ can be updated as $\rho_{i}=c\rho_{i}$, where $c$ is typically
close to $1$, such as $c=1.005$. Through the Equation (15), we are easy to
find that the violation of constraint will lead to higher $s_{i}$, and the
violation level also matters.

####  3.2.3 BNN Inference

There are multiple algorithms that can be used for BNN inference efficiently
by estimating optimize $\theta$. Overall, these algorithms based on Dropout
[12], Bayes by Back-propagation [4], Stein Variational Gradient Descent [21]
and so on.

##### Dropout

By using Dropout layers, nodes will be randomly removed, and the weights of
layer $i$, $\mathbf{w}_{i}$ can be represented as:

$w_{i}=M_{i}*\text{diag}([z_{i,j}]),$

where $M_{i}$ is the parameter matrix for layer $i$ and
$z_{i,j}\sim\text{Bernoulli}(p_{i})$ for each mode $j$. Thus the uncertainty
in $w$ comes from the distribution of $z_{i,j}$ for each variational parameter
[12].

Both Equation (11) and (15) contains the
$E_{q(\mathbf{w}|\theta)}h(\mathbf{w},\theta)$ for some function $h$. It is
easy to get sample $w_{n}$ from Bernoulli distribution $z^{n}_{i,j}$ and
estimate the objective function in Equation (11) and (15).

##### Bayes by Back-propagation

To optimize $\theta$ in $E_{q(\mathbf{w}|\theta)}h(\mathbf{w},\theta)$ by
stochastic gradient descent,
$\frac{\partial}{\partial\theta}E_{q(\mathbf{w}|\theta)}h(\mathbf{w},\theta)$
needs to be derived. Here, Bayes by Back-propagation [4] can be used by
utilizing the reparametrization trick. Here, we assume that $\mathbf{w}\sim
N(\mu,\sigma^{2})\text{, and }\theta=(\mu,\log(\sigma^{2})).$ By the
reparameterization trick, $\mathbf{w}$ can be written as
$\mathbf{w}=\mu+\sigma\epsilon,\epsilon\sim N(0,1)$. Then,
$q(\mathbf{w}|\theta)d\mathbf{w}=q(\epsilon)d\epsilon$ can be used to get the
gradient more easily. The gradient computation result is presented in Equation
(16).

$\frac{\partial}{\partial\theta}E_{q(\mathbf{w}|\theta)}h(\mathbf{w},\theta)=E_{q(\epsilon)}[\frac{\partial
h(\mathbf{w},\theta)}{\partial\mathbf{w}}\frac{\partial\mathbf{w}}{\partial\theta}+\frac{\partial
h(\mathbf{w},\theta)}{\partial\theta}]$ (16)

##### Stein Variational Gradient Descent

By defining a Stein operator and stein discrepancy, we can avoid calculating
the normalization constant in Bayesian rules and get the derivative of the
objective function more easily [21]. The gradient of the objective function
will only depend on

$\displaystyle\nabla(\log(P(\mathbf{w}|\theta))-\lambda\mathbf{f}(\mathbf{x},\mathbf{y}))$
$\displaystyle=\nabla(\log(P(D|\mathbf{w})P(\mathbf{w}))-\lambda\mathbf{f}(\mathbf{x},\mathbf{y}))$
(17) $\displaystyle=\nabla\log P(D|\mathbf{w})+\nabla\log
P(\mathbf{w})-\lambda\nabla\mathbf{f}(\mathbf{x},\mathbf{y}).$

The initial set of particles $w^{0}$ are sampled from the initial estimated
posterior distribution $q^{0}$. Through updating particles iteratively, we can
finally obtain a set of particles $w$ to match target distribution
$p(\mathbf{w})$.

Our proposed PR-BNNs can be solved via different bayesian inference method as
shown above. The complexity of the inference and computation depends on the
Bayesian inference method selected.

##  4 Implementation

In this section, we will evaluate the methods using some simulation examples.
The behavior of the algorithms in these simple examples would be helpful in
understanding the behavior of the proposed algorithm. More specifically, we
will evaluate the case of soft constraint and hard constraint in Section 4.1
and Section 4.2.

###  4.1 Simulation 1: Soft constraint

In this experiment, the same BNNs architecture and example in [29] is used.
The BNNs consist of a single 10-node RBF hidden layer and a Radial basis
function (RBF) as an activation function. For RBF activation function, given
an input $x$, the output after activation could be
$\Phi(x)=e^{-\frac{(x-\mu)^{2}}{\sigma^{2}}}$. For the prior, the paper
suggests an isotropic Gaussian
$\mathbf{w}\sim\mathcal{N}\left(\mathbf{0},\sigma_{\mathbf{w}}^{2}\mathbf{I}\right)$
that is first proposed by MacKay [23] due to the tractability.
$\sigma_{\mathbf{w}}^{2}=1$ is the shared variance for weights $\mathbf{w}$.
We have six data points in the training set. In this experiment, it is
hypothesized that the Knowledge-constrained model can learn those given side
information or knowledge. Bayes by Back-propagation (BBB) [4] is used in this
study. The constraint is in the same format as in [29]: $y\in[2.5,3]$ when
$\text{for }x\in[-0.3,0.3]$. After adding this conditional range constraint,
the posterior predictive results $y_{s}$ are within $[2.5,3]$ as the blue
curve in Figure 1. Without this constraint, the predictive results represented
in black in Figure 1 will not follow the conditional constraint.

Figure 1: Results generated from different Bayesian inference methods

###  4.2 Simulation 2: Hard constraint

##### Simulation Setup

As mentioned, incorporating knowledge constraints may help improve the model
performance when we have noisy or limited training data. In a simple
regression simulation, the training data is noisy and insufficient as shown in
Figure 2(a). In testing, we have data points of $x$ values between $0.08$ and
$1$, while data points available for training are of $x$ values only between
$0.1$ and $0.65$. The true function is $y=(\arctan(20x-10)-\arctan(-10))/3,$
for $x\in[0.08,1]$ as the blue curve shown in Figure 2(a) [1] .

For the constraints, we use: (i) boundary constraint: $a(x)\leq r(x)\leq
b(x)$. (ii) non-decreasing constraint.

  * 1.

Constraint 1: lower bound constraint: $y\geq a(x)=0$ is the lower bound, as
shown in red line of dashes shown in Figure 2. The constraint function can be
presented as $f_{1}(y)=\min\\{0,y\\}$;

  * 2.

Constraint 2: upper bound constraint: $y\leq b(x)=(\log(25x+1))/3+0.05$ is the
upper bound, as shown in the purple line of dashes. The constraint function is
$f_{2}(y)=\min\\{0,b(x)-y\\}$;

  * 3.

Constraint 3: monotone constraint: $\partial f/\partial x_{i}\geq 0$. The
constraint function is $f_{3}(y)=\min\\{0,y^{\prime}\\}$;

##### Methods for Evaluation

Here, we will compare the performance of proposed PR-BNN and two benchmark
models: OC-BNN and Baseline-BNN. More details of each method are presented
here:

  * 1.

Proposed PR-BNN: We apply a BNN with two fully-connected Bayesian linear
layers to fit the 1D function using data points in Figure 2(a). The first
layer has 100 nodes followed by a ReLU activation function and the second
layer is output layer with only 1 node. Besides, we add and update aleatoric
uncertainty to realize automatic parameter tuning for traditional error and
complexity loss.

  * 2.

Baseline-BNN: We use the exact Bayesian Neural Network architecture in the PR-
BNN but without adding the constraints.

  * 3.

OC-BNN: OC-BNN only supports limited types of constraints, such as positive
and negative constraints, where monotone constraints are not directly
supported. Furthermore, there is no guidance on how to set appropriate weights
for multiple constraint types. Therefore, we have made a modification to
convert the constraint on $y$ to the positive constraint on the first-order
gradient of $y$, denoted as $y^{\prime}\in[0,\infty)$. In this case, OC-BNN
can be extended and solve both boundary and monotone constraints. Finally, we
follow the definition of the paper to provide a simple way to incorporate
multiple constraints in the modified prior, including constraints as shown in
Equation (7).

For the case with multiple constraints, we can use rewrite Equation (7) as

$p_{\mathcal{C}}(\mathbf{w})=\frac{1}{Z}p(\mathbf{w})\exp(\Sigma_{i}(c_{i}*f_{i}(\mathbf{y})))=\frac{1}{Z}p(\mathbf{w})\prod\exp(c_{i}*f_{i}(\mathbf{y})),$

where $c_{i}$ is a hyperparameter for output constraint $f_{i}(\mathbf{y})$.

Since they haven’t proposed an efficient algorithm for hard constraints, we
treat the aforementioned three constraints as soft ones by introducing the
corresponding weight $c_{i}$. Parameter tuning is challenging, given that we
have many constraints with a lot of combinations of their weights. Let the
vector $c=[c_{1},c_{2},c_{3}]$ demotes weights for constraint 1, 2 and 3
respectively. To implement the OC-BNN work into this simulation case, we use
multiple parameter combinations of $c$ and record results.

##### Results

To evaluate the performance, the following evaluation metrics are used 1) Mean
Squared Error (MSE): MSE is used to evaluate the accuracy of mean prediction.
2) Epistemic Uncertainty Standard deviation (STD): Epistemic uncertainty is
computed by the standard deviation of the Monte Carlo sampling of all BNN
methods averaged over the input domain. 3) Violation value for constraint $i$
as $v_{i}$: $v_{i}$ is computed by the expected value of the constraint
function as $v_{i}=-E_{w}(f_{i}(y))$, where $f_{i}$ is the constrained
function. 4) Number of Violations for constraint $i$ as $n_{i}$: $n_{i}$ shows
the number of constraint $i$ violation predictions among 200 testing samples.
The results of these four groups of metrics of PR-BNN, OC-BNN, and BNN are
shown in Table 2.

For PR-BNN with hard constraints, all three constraints have been satisfied,
and the detailed results for each constraint are shown in the last row of
Table 2. The dual variables $s_{i}$ can be considered as the weights of
corresponding constraints, which are automatically decided by the augmented
Lagrangian algorithm. The initial setting is $s=[1,1,1]$, and after training,
final vector $s$ is $[4.36,14.41,6.41]$. Their rounded numbers are shown in
the last row of Table 2. We can observe that $s_{2}=14.41$ is far greater than
the other two dual variable values, indicating that the upper bound constraint
violation is more severe in the training stages and large $s_{2}$ accelerate
the process of learning this constraint.

Comparing PR-BNN and BNN, MSE is much smaller for the proposed PR-BNN method,
given that constraints provide useful information, especially for the range
without any data points. Without constraints, a severe violation of the upper
bound is shown in Figure 2(a), while most predictions and mean predictions
satisfy the bound constraint after constraints injection as shown in 2(d). The
MSE loss is reduced from 0.2320 in the BNN model to 0.0052 in the PR-BNN model
after adding constraints. The epistemic uncertainty is reduced from 0.059 in
the unconstrained model to 0.040 in the hard constrained model, benefiting
from the knowledge, and we can easy to observe this uncertainty reduction
through Figure 2(a) and (d).

For OC-BNN, we would like to evaluate different sets of weight coefficients to
combine the constraints. From Table 2, we can conclude that the upper bound
(Constraint 2) and Monotone constraint (Constraint 3) are competing during
training. Due to the space limit, we only include the visualization for
hyperparameter settings $c=[1,1,2]\text{ and }c=[1,1,4]$ in Figure 2(b) and
Figure 2(c) respectively but constraints violation in Table 2 shows details
for more experiments, where $v_{i}$ represents the violation value for
constraint $i$. From the visualization samples in Figure 2(b), we can observe
that even though we add double weights on the monotone constraint (Constraint
3), the upper bound constraint (Constraint 2) still pushes the predicted curve
too low to violate the monotone constraint. Thus, when we have multiple
constraints, it is hard to satisfy all of them via the OC-BNN without time-
consuming parameter tuning, and hard to find a balanced setting. As shown in
Figure 2(c), the prediction becomes better when we add more weights on
Constraint 3, but a very small violation value is still shown in the
corresponding row of Table 2. Besides, under the $[1,1,8]$ parameter setting,
the upper bound will be violated greatly since the weight of the monotone
constraint is too high. From the comparison in Table 2 we can see the proposed
PR-BNN can satisfy all the hard constraints without manually parameter tuning.

Table 2: Performance comparison on testing set  Framework | c/s | $v_{1}$ | $v_{3}$ | $v_{3}$ | $n_{1}$ | $n_{2}$ | $n_{3}$ | MSE | STD
---|---|---|---|---|---|---|---|---|---
BNN | - | 0 | 0.22 | 0 | 0 | 77 | 0 | 2.3e-1 | 0.059
OCBNN with different parameter vector (c) | [1,1,1] | 0 | 0 | 2.0e-3 | 0 | 0 | 84 | 3.6e-2 | 0.047
[1,1,2] | 0 | 0 | 1.4e-3 | 0 | 0 | 86 | 3.2e-2 | 0.047
[1,2,1] | 0 | 0 | 2.7e-3 | 0 | 0 | 84 | 5.9e-2 | 0.046
[2,1,1] | 0 | 0 | 1.7e-3 | 0 | 0 | 84 | 3.3e-2 | 0.052
[1,1,4] | 0 | 0 | 3.9e-4 | 0 | 0 | 83 | 8.3e-3 | 0.060
[1,1,8] | 0 | 0.11 | 0 | 0 | 67 | 0 | 8.6e-2 | 0.079
PR–BNN (s) | [4,14,6] | 0 | 0 | 0 | 0 | 0 | 0 | 5.2e-3 | 0.040
Figure 2: Experiment details and results comparison

##  5 Case studies

In this section, we will evaluate the proposed methods in two case studies.
First, we will investigate the effect of the knowledge constraints added in
the Aviation Landing Status Prediction in Section 5.1. Second, we will
evaluate how such constraints can be added to solar power prediction using the
knowledge from the inverter machine in Section 5.2.

###  5.1 Aviation Landing Trajectory Prediction

In air traffic management and safety control, aircraft landing is one of the
riskiest phases of flight with multiple possible adversities, such as
misalignment, hard landing, runway overrun, etc [28]. Therefore, in this work,
we focus on predicting the aircraft landing tracks with real-world aviation
radar recordings. We obtain the flight track data from the Sherlock Data
Warehouse [2]. Sherlock is a platform to collect, archive, process, and query
aviation datasets such that to support ATM research [25]. It contains aviation
flight data, meteorological data, and airline operations data from various
sources (e.g., FAA and NOAA). The flight tracks are the time-series aircraft
coordinates landed at runway 26R at KATL (Hartsfield - Jackson Atlanta
International Airport) from January 1st, 2019 to January 31st, 2019, for the
case studies. We perform data preprocessing to filter out the last several
timestamps since the tracks from Sherlock are the recordings for the entire
flight from take-off to landing. The preprocessing of the landing tracks can
be summarized as follows: a) perform altitude filtering to determine the
touchdown point (landing point on the runway) from the time series; b) conduct
runway checks to make sure the given track landed at the correct runway; c) do
completion checks to make sure the processed flight tracks are valid and only
including a fix time horizon prior to touchdown. Upon finishing data filtering
and preprocessing, the dataset used for model training, validation, and
testing consists of 7,542 landing trajectories with a fixed horizon of 50
timestamps before the touchdown.

Here, 5 different variables are collected to picture the flight status in the
final 50 seconds before touchdown. Altitude, latitude, and longitude can
describe the location of the airplane and we have 2-D landing trajectory
sequence examples including longitude and latitude for two flights shown in
Figure 3; airspeed and course can characterize the speed of the airplane
relative to the air mass [14] and heading direction, respectively. Data for
all variables is normalized using Min-Max scalar. Without loss of generality,
we use $x$ to represent any variable in this dataset and the normalization can
be shown in Equation (18).

$x_{min\\_max}=\frac{x-x_{min}}{x_{max}-x_{min}}$ (18)

However, the dataset generated by multiple sensors from airports and aircraft
is rather noisy, which troubles accurate prediction during the landing phase.
As a safety-critical problem, landing speed forecasting and uncertainty
quantification are of great importance for landing risk mitigation.
Integrating common knowledge in BNNs in the model of the landing phase helps
reduce prediction uncertainty and improve safety. Furthermore, accurate
prediction of the normal landing trajectory with the knowledge-constrained BNN
prediction can also help to detect anomalies in the landing phase that
seriously violate such knowledge constraints and send early warnings to
control tower.

Figure 3: Landing Trajectory Examples Figure 4: Bayesian LSTM Structure

Based on this dataset, we would use the first 35 seconds of data to predict
the flight status of the last 15 seconds before landing touchdown. A long
short-term memory model (LSTM) is used to model this time series data and via
Bayesian LSTM as described in [10], with modifications to incorporate
constraints.

The Bayesian LSTM model used in this case study consists of one Bayesian LSTM
layer and simply two fully connected layers.

We have 4 different constraints in this case:

  * 1.

Entire Trajectory: Altitude decreasing constrain. In practice, we allow
altitude fluctuation. The implementation of this constraint is described as:
In the last 15 seconds, the altitude at $t_{i}$ should greater than the
altitude at $t_{i+3}$.

  * 2.

Entire Trajectory: Range constraints on all five variables; For speed and
altitude, they should be non-negative. For course, longitude and latitude,
range constraints can be set based on the location and direction of the
runway.

  * 3.

Landing Point Constraint: Altitude range constraint in landing; Altitude range
will be very narrow in touchdown moment due to the narrow change in altitudes
of the runway.

  * 4.

Landing Point Constraint: Speed range constraint in landing; The touchdown
airspeed should be within a reasonable and safe range.

Here, we have added the constraints as the soft constraints, given that there
can actually be very rare violations of these constraints during some
uncontrollable rare cases. In this case, $\xi$ in Equation (9) is a slack
variable, instead of a known constant.

To show the performance of the PR-BNN model, we compare the 15-second sequence
prediction and last point prediction before and after adding constraints. The
Mean Squared Error (MSE), the uncertainty which is represented by Standard
Deviation (STD) and The continuous ranked probability score (CRPS) a
performance measurement for probabilistic forecasts is shown in Table 3, Table
4, and Table 5, respectively.

Table 3: Mean squared error results comparison Model | MSE without Constraints | MSE with Constraints
---|---|---
Variables | sequence | last point | sequence | last point
latitude | 1.806E-04 | 1.500E-04 | 1.813E-04 | 1.514E-04
longitude | 9.160E-05 | 1.838E-04 | 7.050E-05 | 9.710E-05
altitude | 3.079E-04 | 5.186E-04 | 2.450E-04 | 2.655E-04
speed | 1.773E-04 | 1.870E-04 | 1.709E-04 | 1.813E-04
course | 7.563E-04 | 5.462E-04 | 7.591E-04 | 5.477E-04
overall | 3.027E-04 | 3.171E-04 | 2.853E-04 | 2.486E-04
Table 4: Standard deviation results comparison Model | STD without Constraints | STD with Constraints
---|---|---
Variables | sequence | last point | sequence | last point
latitude | 2.907E-03 | 2.760E-03 | 2.411E-03 | 2.379E-03
longitude | 8.241E-03 | 1.066E-02 | 5.913E-03 | 4.854E-03
altitude | 8.149E-03 | 1.178E-02 | 5.182E-03 | 4.209E-03
speed | 3.983E-03 | 4.235E-03 | 3.255E-03 | 3.150E-03
course | 1.676E-03 | 1.554E-03 | 1.283E-03 | 1.297E-03
overall | 4.991E-03 | 6.198E-03 | 3.609E-03 | 3.178E-03
Table 5: CRPS results comparison Model | CRPS without Constraints | CRPS with Constraints
---|---|---
Variables | sequence | last point | sequence | last point
latitude | 1.091E-02 | 8.628E-03 | 9.305E-03 | 8.451E-03
longitude | 1.848E-02 | 1.641E-02 | 5.203E-03 | 5.990E-03
altitude | 1.841E-02 | 1.442E-02 | 1.008E-02 | 1.064E-02
speed | 2.180E-02 | 2.130E-02 | 8.682E-03 | 8.263E-03
course | 2.043E-02 | 1.670E-02 | 2.035E-02 | 1.686E-02
overall | 1.801E-02 | 1.549E-02 | 1.072E-02 | 1.004E-02
Figure 5: (a): Comparison for altitude prediction results based on 3 different
flights data; (b): Altitude constraint violation comparison between
constrained (454 violations) and unconstrained (675 violations) model.

From Table 3, the overall accuracy of the prediction of both the 15-second
sequence and the last landing point is higher after adding knowledge
constraints. Baseline BNN model obtains the lower error in latitude and course
predictions but those results are pretty close. Based on the knowledge of
these two variables, they are both related to the direction and one flight is
almost approaching along a straight line to the runway. Thus, the latitude and
course data near the landing points is noisy and no models can gain
significant advantage over the others.

From Table 4, an obvious reduction in the uncertainty of predictions is
exhibited when we use the knowledge-constrained model. Table 5 demonstrate the
probabilistic prediction performance improvement of our constrained framework
via CRPS scores. One of our constraints is boundary constraint on the altitude
of the last landing point, which is marked as the range in Figure 5 (a) by two
red lines of dashes. We can see the constraint violation of the green flight
without knowledge constraints, while the prediction satisfies the constraint
after applying the constraints. Figure 5 (b) shows the reduction of altitude
constraint violation with knowledge. Since in last 15 second, most fluctuation
of latitude is noise, we visualize the vertical 2D coordinate system including
altitude and longitude. Figure 6 shows the performance difference of BNN and
PR-BNN models and the corresponding ground true values for this flight as
approaching.

Figure 6: Sample prediction results comparison with uncertainty: The shaded
area shows the last 15-second predictions for two model and the ground true
power output values. The height and width of the markers are proportional to
the uncertainty of altitude and longitude predictions.

###  5.2 Photovoltaic Solar Energy Modelling

In a solar power plant, accurate prediction of the solar Photovoltaic (PV)
plant power is important for power modeling, forecasting, and anomaly
detection. In this case study, we use some ambient variables and operation
rules provided by professionals to explore the AC power generated in the
plant. Generally, a plant consists of multiple solar power inverters and each
inverter has several modules. In this dataset, we explore one inverter in a
solar plant using the proposed PR-BNN framework.

Our data is provided by Arizona Public Service in the red rock PV plant. The
core components in the solar energy system can be illustrated in Figure 7.
Various sensor data are available to measure the conditions at the system
level and inverter level, and the weather station information is also provided
on the site.

Figure 7: Solar Energy System

In this dataset with the 15-minute resolution, our target variable is the AC
power (W) of inverter A. The summer season is pretty long at the specific
plant location in Arizona, and the extreme heat in summer often causes issues
in this site, as experts described. Thus, we select only summer season data
from May to October in 2019 and 2020. We have 20698 samples used in training
and 6505 in testing. The independent variables are organized into two levels.
First, system-level variables contain some control-related variables such as
curtailment, combiner switch, feeder switch, and some weather-related factors
including wind speed, and irradiance. Second, the inverter-level data includes
two binary variables related to IGBT (Insulated-Gate Bipolar Transistor)
temperatures as shown in the second and third variables of Table 6. IGBT is a
core component of each module, and indicates the sometimes indicating the
working condition of the corresponding module. We have some knowledge about
the IGBT temperature. For example, if IGBT_derate value is 1, it means that at
least one module in inverter A has the derate condition due to the high IGBT
temperature. Daytime and date variables are included to represent the temporal
relationship in the simple BNN model. All variables used in this case study
are listed in Table 6.

Table 6: Variables and description Level | Viarable | Description | Unit
---|---|---|---
Inverter | Pac | AC Power | w
IGBT_derate |  $>$ IGBT temperature derate threshold | -
IGBT_off |  $>$ IGBT temperature off threshold | -
System | curtail_utility | Utility curtailment power set point | mw
curtail_on | Curtailment OFF/ON | -
P_watts | Main Xfmr3 Sec. power | mw
Feeder_B | Feeder 2 Breaker | -
WS | Wind Speed | m/s
IPOA | Irradiance - plane of array | w/m^2
day_label | The day in a year | -
time_label | The timestamp in a day | -

In the data prepossessing, min-max scalar has been applied for both input and
target variables for normalization as shown in Equation (18).

In the solar energy plant, it is well known that solar irradiance has a direct
impact on the power output, and it can explain most variability of the AC
power output, as shown in Figure 8. But for some variability, especially in
the dashed rectangular areas, the solar irradiance alone cannot explain such
variability. We would like to include other variables and even professional
knowledge into the model to learn complex patterns.

Figure 8: Sample variable visualization: x-axis is the time stamp; the y-axis
represents the normalized irradiance value and the ground true AC power output
values. The comparison of these two variables shows pattern mismatch during
the peak time or daytime on some days.

For the knowledge constraint, we follow the operational guide by APS that if
the Max IGBT temperature reaches a threshold given by the expert, as shown in
Figure 9, the corresponding module will be derated to protect the machine so
that IGBT temperature could be back to the safe value. As a result, the AC
power will go down given that other conditions are unchanged, such as
irradiance.

Figure 9: Derate condition based on Max IGBT temperature

A simple BNN with three Bayesian linear layers is applied as a benchmark BNN
model as shown in Figure 10. The three layers have 64, 32, and 1 node,
respectively, with Gaussian prior for Neural Network parameters, and the first
two layers are followed by the ReLU activation function. Based on the
practical goals, the daytime (9:30 am - 5:00 pm) and peak time (10:00 am -
2:00 pm) performances are of most interest. Thus, we also use Mean Square
Error (MSE), Root Mean Squared Error (RMSE) based on rescaled data to measure
the mean prediction accuracy, and Continuous Ranked Probability Score (CRPS)
based on rescaled data to evaluate the probabilistic forecasting quality.
Evaluation comparisons between the Benchmark BNN model and the proposed PR-BNN
model are shown in Table 7. After adding the derate constraint, the MSE is
reduced from 0.0018 to 0.0011. Performance improvements can be also observed
by the RMSE and CRPS.

Figure 10: Baseline Bayesian Linear Neural Network Structure Table 7: Performance comparison based on different time ranges Model | BNN | PR-BNN
---|---|---
Time | All | Daytime | Peaktime | All | Daytime | Peaktime
MSE | 0.0018 | 0.0045 | 0.0037 | 0.0011 | 0.0026 | 0.0017
RMSE | 0.0419 | 0.0672 | 0.0611 | 0.0327 | 0.0505 | 0.0415
CRPS | 0.0186 | 0.0448 | 0.0450 | 0.0114 | 0.0233 | 0.0212

Sample prediction results from both with and without constraint models for six
days are visualized in Figure 11. From the comparison, the PR-BNN model with
derate-constraint learns the derate effect in the solar power system more
accurately, which leads to improvement in overall accuracy in prediction.

The review of the corresponding solar irradiance as shown in Figure 8 can make
clear that the shape property of BNN baseline prediction is more similar to
the irradiance curve shape but cannot learn the derate pattern from the data
itself. Thus adding derate knowledge can help improve the accuracy of the
regression and further may help detect the anomaly in the energy plant.

Figure 11: Sample variable visualization: The x-axis is timestamp; The y-axis
represents the normalized ground true and predictions of two models. The
y-axis value of the red scatters no real meaning, but it indicates that there
should be a manual derate at a certain timestamp with a red dot.

##  6 Conclusion

The performances of deep learning models can be benefited from rich domain
knowledge or side information via BNNs, even for those cases with poor or
limited datasets in HD problems. The PR-BNN with a general form of constraint
functions based on posterior regularization is proposed. Due to the various
properties of domain knowledge, we include both types of soft and hard
constraints in this work and develop efficient optimization algorithms for
both cases.

After implementing both soft and hard knowledge constraints with proposed
algorithms, favorable results are demonstrated by our simulation results,
showing promising potentials for performance and accuracy improvement. Besides
simulation results, the proposed method is evaluated through two real-world
case studies from different fields, namely trajectory prediction and solar
energy modeling. Along with the accuracy metric, CRPS is adopted to evaluate
and compare the probabilistic forecast performance of BNNs and PR-BNNs. The
two case studies illustrate how to apply various professional or domain
knowledge to real-world applications and show desirable performance on both
deterministic and probabilistic prediction accuracy.

## Acknowledgment

The research reported in this paper was supported by funds from NASA
University Leadership Initiative program (Contract No. NNX17AJ86A, PI:
Yongming Liu, Technical Officer: Anupa Bajwa) and funds from Department of
Energy (Contract No. DE-EE0009354, PI: Hao Yan). The support is gratefully
acknowledged.

## References

  * Agrell [2019] Christian Agrell.  Gaussian processes with linear operator inequality constraints.  _arXiv:1901.03134 [cs, stat]_ , September 2019.  URL http://arxiv.org/abs/1901.03134.  arXiv: 1901.03134.
  * Arneson et al. [2019] Heather M Arneson, Pallavi Hegde, Michael E La Scola, Antony D Evans, Richard M Keller, and John E Schade.  Sherlock data warehouse.  2019\.
  * Blei et al. [2017] David M Blei, Alp Kucukelbir, and Jon D McAuliffe.  Variational inference: A review for statisticians.  _Journal of the American statistical Association_ , 112(518):859–877, 2017.
  * Blundell et al. [2015] Charles Blundell, Julien Cornebise, Koray Kavukcuoglu, and Daan Wierstra.  Weight uncertainty in neural network.  In _International Conference on Machine Learning_ , pages 1613–1622. PMLR, 2015.
  * Bouzerdoum and Pattison [1993] Abdesselam Bouzerdoum and Tim R Pattison.  Neural network for quadratic optimization with bound constraints.  _IEEE transactions on neural networks_ , 4(2):293–304, 1993.
  * Cabiscol [2019] Javier Antorán Cabiscol.  Understanding uncertainty in bayesian neural networks.  _Master of Philosophy (University of Cambridge)_ , 2019.
  * Dalton and Dougherty [2013] Lori A Dalton and Edward R Dougherty.  Optimal classifiers with minimum expected error within a bayesian framework—part i: Discrete and gaussian models.  _Pattern Recognition_ , 46(5):1301–1314, 2013\.
  * Dreyfus [2005] Gérard Dreyfus.  _Neural networks: methodology and applications_.  Springer Science & Business Media, 2005.
  * Esfahani and Dougherty [2014] Mohammad Shahrokh Esfahani and Edward R. Dougherty.  Incorporation of Biological Pathway Knowledge in the Construction of Priors for Optimal Bayesian Classification.  _IEEE/ACM Transactions on Computational Biology and Bioinformatics_ , 11(1), January 2014.  ISSN 1557-9964.  doi: 10.1109/TCBB.2013.143.
  * Esposito [2020] Piero Esposito.  Blitz - bayesian layers in torch zoo (a bayesian deep learing library for torch).  https://github.com/piEsposito/blitz-bayesian-deep-learning/, 2020\.
  * Feelders and Van Der Gaag [2012] Ad Feelders and Linda C Van Der Gaag.  Learning bayesian network parameters with prior knowledge about context-specific qualitative influences.  _arXiv preprint arXiv:1207.1387_ , 2012.
  * Gal and Ghahramani [2016] Yarin Gal and Zoubin Ghahramani.  Dropout as a bayesian approximation: Representing model uncertainty in deep learning.  In _international conference on machine learning_ , pages 1050–1059. PMLR, 2016.
  * Ganchev et al. [2010] Kuzman Ganchev, Joao Graça, Jennifer Gillenwater, and Ben Taskar.  Posterior regularization for structured latent variable models.  _The Journal of Machine Learning Research_ , 11:2001–2049, 2010.
  * Gracey [1980] William Gracey.  Measurement of aircraft speed and altitude.  Technical report, National Aeronautics and Space Administration Hampton VA Langley Reseach Center, 1980.
  * Hastings [1970] W Keith Hastings.  Monte carlo sampling methods using markov chains and their applications.  1970\.
  * Hüllermeier and Waegeman [2021] Eyke Hüllermeier and Willem Waegeman.  Aleatoric and epistemic uncertainty in machine learning: An introduction to concepts and methods.  _Machine Learning_ , 110(3):457–506, 2021.
  * Kingma and Welling [2013] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  _arXiv preprint arXiv:1312.6114_ , 2013.
  * Krzysztofowicz [1999] Roman Krzysztofowicz.  Bayesian theory of probabilistic forecasting via deterministic hydrologic model.  _Water Resour. Res_ , pages 2739–2750, 1999.
  * Kullback and Leibler [1951] Solomon Kullback and Richard A Leibler.  On information and sufficiency.  _The annals of mathematical statistics_ , 22(1):79–86, 1951.
  * Liao and Ji [2009] Wenhui Liao and Qiang Ji.  Learning bayesian network parameters under incomplete data with domain knowledge.  _Pattern Recognition_ , 42(11):3046–3056, 2009\.
  * Liu and Wang [2016] Qiang Liu and Dilin Wang.  Stein variational gradient descent: A general purpose bayesian inference algorithm.  _Advances in neural information processing systems_ , 29, 2016.
  * Liu et al. [2021] Xin-Wei Liu, Yu-Hong Dai, Ya-Kui Huang, and Jie Sun.  A novel augmented lagrangian method of multipliers for optimization with general inequality constraints.  _arXiv preprint arXiv:2106.15044_ , 2021.
  * MacKay [1995] David JC MacKay.  Probable networks and plausible predictions-a review of practical bayesian methods for supervised neural networks.  _Network: computation in neural systems_ , 6(3):469, 1995.
  * McDermott and Wikle [2019] Patrick L. McDermott and Christopher K. Wikle.  Bayesian Recurrent Neural Network Models for Forecasting and Quantifying Uncertainty in Spatial-Temporal Data.  _Entropy_ , 21(2):184, February 2019.  doi: 10.3390/e21020184.  URL https://www.mdpi.com/1099-4300/21/2/184.  Number: 2 Publisher: Multidisciplinary Digital Publishing Institute.
  * Pang et al. [2022] Yutian Pang, Xinyu Zhao, Jueming Hu, Hao Yan, and Yongming Liu.  Bayesian spatio-temporal graph transformer network (b-star) for multi-aircraft trajectory prediction.  _Knowledge-Based Systems_ , page 108998, 2022.
  * Specht [1990] Donald F Specht.  Probabilistic neural networks.  _Neural networks_ , 3(1):109–118, 1990.
  * Tai et al. [2015] Cheng Tai, Tong Xiao, Yi Zhang, Xiaogang Wang, et al.  Convolutional neural networks with low-rank regularization.  _arXiv preprint arXiv:1511.06067_ , 2015.
  * Wang et al. [2019] Yuhao Wang, Yutian Pang, Yongming Liu, Parikshit Dutta, and Bong-Jun Yang.  Aircraft trajectory prediction and risk assessment using bayesian updating.  In _Aiaa aviation 2019 forum_ , page 2936, 2019.
  * Yang et al. [2021] Wanqian Yang, Lars Lorch, Moritz A. Graule, Himabindu Lakkaraju, and Finale Doshi-Velez.  Incorporating Interpretable Output Constraints in Bayesian Neural Networks.  _arXiv:2010.10969 [cs, stat]_ , January 2021.  URL http://arxiv.org/abs/2010.10969.
  * Yang et al. [2020] Zebin Yang, Aijun Zhang, and Agus Sudjianto.  Enhancing explainability of neural networks through architecture constraints.  _IEEE Transactions on Neural Networks and Learning Systems_ , 2020\.
  * Zhang et al. [2018] Cheng Zhang, Judith Bütepage, Hedvig Kjellström, and Stephan Mandt.  Advances in variational inference.  _IEEE transactions on pattern analysis and machine intelligence_ , 41(8):2008–2026, 2018.
  * Zhu et al. [2014] Jun Zhu, Ning Chen, and Eric P Xing.  Bayesian inference with posterior regularization and applications to infinite latent svms.  _The Journal of Machine Learning Research_ , 15(1):1799–1847, 2014.