# Can physical information aid the generalization ability of Neural Networks
for hydraulic modeling?

Gianmarco Guglielmo  Andrea Montessori  Jean-Michel Tucny  Michele La Rocca
Pietro Prestininzi

###### Abstract

Application of Neural Networks to river hydraulics is fledgling, despite the
field suffering from data scarcity, a challenge for machine learning
techniques. Consequently, many purely data-driven Neural Networks proved to
lack predictive capabilities. In this work, we propose to mitigate such
problem by introducing physical information into the training phase. The idea
is borrowed from Physics-Informed Neural Networks which have been recently
proposed in other contexts. Physics-Informed Neural Networks embed physical
information in the form of the residual of the Partial Differential Equations
(PDEs) governing the phenomenon and, as such, are conceived as neural solvers,
i.e. an alternative to traditional numerical solvers. Such approach is seldom
suitable for environmental hydraulics, where epistemic uncertainties are
large, and computing residuals of PDEs exhibits difficulties similar to those
faced by classical numerical methods. Instead, we envisaged the employment of
Neural Networks as neural operators, featuring physical constraints formulated
without resorting to PDEs. The proposed novel methodology shares similarities
with data augmentation and regularization. We show that incorporating such
soft physical information can improve predictive capabilities.

###### keywords:

Physics-Informed Neural Networks , river hydraulics modeling , free surface ,
soft physical information , neural operator , data augmentation

††journal:

[first]organization=Roma Tre University,addressline=Department of Civil,
Computer Science, and Aeronautical Technologies Engineering, Via Vito Volterra
62, 00146, city=Rome, country=Italy

[second]organization=Center for Life Nano- & Neuro-Science, Fondazione
Istituto Italiano di Tecnologia (IIT),addressline=Viale Regina Elena 295,
00161, city=Rome, country=Italy

##  1 Introduction

As common to environmental hydraulics, river flood mapping requires models
able to explore scenarios not entirely encompassed by the set of available
observations, i.e. featuring predictive/generalization capabilities. An
illustrative example is the simulation of high hydraulic hazard scenarios by
means of models calibrated only on observations of medium to low flow regimes
[16].

Recently, the need for reliable flood maps of ungauged basins has become
urgent [7, 13, 20], the latter ranging from small catchments to scarcely
populated large regions in developing countries. In addition to the
difficulties related to the lack of high quality and/or scarce measurements,
mapping vast flood-prone areas by means of physically based models requires a
considerable computational burden, even assuming drastic conceptual
simplifications [36]. Classical modeling approaches [28] may thus become
unfeasible in such cases, and resorting to the exploitation of similarities
with other basins has been envisaged. Indeed, Machine Learning (ML) has been
proposed [15] as a fruitful way to overcome the above issues, with Neural
Networks (NNs) specifically employed in the context of environmental
hydraulics [27]. A comprehensive review has been developed by [5]. An example
of prediction of urban pluvial flood by means of Convolutional Neural Networks
(CNNs) is given by [30], while a CNN for fluvial flood inundation is developed
by [24].

However, the above referenced first attempts seem to point out a substantial
inability of the trained NNs to generalize, that is to provide a reasonable
prediction for scenarios even slightly different from the ones over which they
were trained on [35]. The main reason of such flaw is the lack of sufficiently
large training datasets [18], resulting in overfitted models. An effect
closely related to overfitting and, as such, ascribable to the scarcity of
calibration data, is the frequent violation of conservation laws [6],
resulting in marked non-physical results.

A recent new paradigm is represented by Physics‐Informed Machine Learning
(PIML) [25]. PIML consists of augmenting the informational content of the
training datasets by introducing physically‐based constraints into existing ML
models. By doing so, it is theoretically possible to reduce the machine’s
learning time, increase its generalization capabilities [22], and obtain
physically‐based results. The physical constraints typically require the
machine output to comply with conservation equations [21].

The widespread approach of PIML generated a new class of NNs, namely the
Physics-Informed Neural Networks (PINNs), introduced by [38] and intended as
neural solvers [19], i.e. a NN specifically designed to find the behavior of a
unique physical phenomenon described by a known differential problem. PINNs
allow for solving both forward and inverse problems governed by a known
partial differential equation (PDE), initial and boundary conditions [31].
They typically modify the loss function by adding a term containing the
residual of the PDE. Most of the highly specialized algorithms proposed for
computing the derivatives of PDEs rely on Automatic Differentiation [4],
although Numerical Differentiation or a combination of both approaches [11]
have also been considered.

Recently, PINNs have been applied to a wide range of benchmark problems of
fluid mechanics [8], in particular as neural solvers for Navier-Stokes
equations [23], and have shown the potential to become more practical than
classical techniques in many circumstances [14]. The advantages of such
approach, compared to numerical models, are certainly the possibility of
seamlessly integrating empirical data into the model and the ability to
address problems that are mathematically ill-posed. Additionally, it enables
tackling issues whose dimensionality is so high that the use of traditional
solvers would result in dramatically elevated computational burden.

However, we believe employing PINNs as neural solvers is hardly applicable to
river hydraulics. Indeed, PINNs are designed to approximate the solution to a
PDE within a specific domain, and every change requires training a new PINN.
Furthermore, what is being sought is not a surrogate for the numerical model
solving for a single instance of a physical system, which would be subject to
the same, if not additional, limitations and difficulties previously
described. Instead, the aim is to design a neural operator [19] that can
effectively approximate an entire class of physical problems.

In the field of flood simulations, the challenge for ML applications is to
employ scarce observations for the training phase while being able to
generalize reliably over new scenarios. Therefore, in this work we assess
whether the performance of a deep learning model can be enhanced by
incorporating soft physical information, i.e. without resorting to the
calculation of the residual of a differential equation (hard physical
information). This feature can be a significant advantage in applications
characterized by a large epistemic uncertainty, i.e. the conceptual model of
the phenomenon is partially or entirely unknown. In the context of flood
simulations, where hydraulically relevant distributed parameters such as
roughness, lithology, topography etc. pose significant uncertainties,
introducing hard physical information like either the residual of the 2D
Shallow-Water equations [37] or the 1D Saint-Venant equations [32], may not be
advisable.

In this work, we extend the concept of PINNs with respect to two fundamental
points:

\- the NN is not used to obtain a tailored solution of a single problem
(neural solver), instead the aim is to construct a neural operator capable of
capturing the relationships between input and output functions across an
entire class of physical phenomena (neural operator);

\- the insertion of physical information does not necessarily have to be
linked to the residual of the differential equation (hard physical
information), but can involve other physical properties of the system and the
data (soft physical information). This method can be interpreted as a special
form of physical data augmentation, involving both inputs and outputs.

We then assessed the gain in generalization capability of such physical aided
neural operators compared to the purely Data-Driven (DD) ones in the context
of environmental hydraulics.

The problem solved by our NNs consists of a simple yet non-trivial physical
problem, namely the reconstruction of a steady state, one-dimensional, water
surface profile in a rectangular channel. This problem was investigated by [9]
for steady flow using PINNs as neural solvers. Its simplicity stems from both
the low dimensionality of the feature space, which allows for the adoption of
small NNs, and the possibility to easily generate dataset of exact solutions
at will. Moreover, the choice of an idealized case over a more realistic one
mitigates the risk of producing biased results as a consequence of possible
prevailing aspects of the specific case, thus preventing the effects of the
proposed modifications to be properly assessed. The highly informative content
of the analyzed problem comes mainly from the hidden complexity of the
underlying physics, associated to the possible occurrence of a hydraulic jump,
a phenomenon which does not belong to the set of solutions of the underlying
differential equations.

The paper is structured as follows: section Methodology contains the overall
plan followed in the work; the Numerical experiments section details the
structures of the NNs and describes how the soft physical information is
employed; then the Discussion of results follows, defining the assessment
procedures. Conclusions are then drawn and perspectives are advanced.

##  2 Methodology

###  2.1 Data-Driven and Physics-Informed Neural Networks

Neural Networks [29] aim to capture (i.e. learn) complex mappings between
inputs $\mathbf{x}$ and outputs $\mathbf{\hat{y}}$, described by an unknown
non-linear function $G^{*}$:

$G^{*}(\mathbf{x})=\mathbf{\hat{y}}$ (1)

By tuning the values of the parameter vector $\bm{\theta}$, the NN learns to
mimic $G^{*}$ through the approximate function $G$:

$G(\mathbf{x};\bm{\theta})=\mathbf{y}$ (2)

where $\mathbf{y}$ are the predicted values.

The NN parameters are learned from data by leveraging information derived from
the training set, namely the couples $(\mathbf{x},\mathbf{\hat{y}})$. In a
purely DD training process, the network calibrates its parameters by
minimizing a loss function $\mathcal{L}$ between $\mathbf{y}$ and
$\mathbf{\hat{y}}$:

$arg\ min_{\bm{\theta}}\mathcal{L}(\mathbf{y},\mathbf{\hat{y}})$ (3)

A common choice for the loss function in a DD model for a regression problem
(the output variables are real unbounded numbers) is the Mean Squared Error
(MSE):

$MSE=\frac{\sum_{i=1}^{N}(y_{i}-\hat{y_{i}})^{2}}{N}$ (4)

but other metrics can be used depending on the nature of the data and the
goals of the model. For further insights into Feed-Forward Neural Networks
(FFNNs) and their training process, additional information can be found in A.

PINNs are a new paradigm of ML models that combine principles from physics-
based modeling with the flexibility of NNs. The key concept is to embed well-
known physical laws into the training phase of a NN, integrating data with
mathematical models.

PINNs have been used as neural solvers for physical systems governed by a
known set of Partial Differential Equations (PDEs), addressing both forward
problems (solving the equation) and inverse problems (determining unknown
parameters). PINNs can even handle situations where the problem is
mathematically ill-posed and therefore classical numerical methods are
unfeasible.

The strategy pursued by PINNs lies in building the loss function as follows:

$\mathcal{L}=\lambda\cdot\mathcal{L}_{DD}(\mathbf{y},\mathbf{\hat{y}})+(1-\lambda)\cdot\mathcal{L}_{P}(\mathbf{x},\mathbf{y})$
(5)

where $\mathcal{L}_{DD}$ is a metric which measures the distance between
predicted $\mathbf{y}$ and observed data $\mathbf{\hat{y}}$ (e.g. MSE,
equation (4)), $\mathcal{L}_{P}$ measures the residual of the PDE sampled at a
set of points within the PDE domain, $\lambda$ balances the data and physical
contributions, ranging between 0 and 1.

###  2.2 Soft physical information

As for the PINNs, the proposed inclusion of soft physical information in the
training phase of a neural operator is achieved through an additional physical
term $\mathcal{L}_{P}$ in the loss function, as shown in equation (5). Various
physical principles can be encoded into this term, enabling the model to
capture a broader spectrum of physical behaviors and relationships hidden in
the data.

While the term $\mathcal{L}_{DD}$ still depends solely on the predicted
outputs $\mathbf{y}$ and the true outputs $\mathbf{\hat{y}}$, the term
$\mathcal{L}_{P}$ is a metric that can also involve the inputs $\mathbf{x}$.
Therefore, the loss function $\mathcal{L}$ incorporating soft physical
information reads:

$\mathcal{L}=\lambda\cdot\mathcal{L}_{DD}(\mathbf{y},\mathbf{\hat{y}})+(1-\lambda)\cdot\mathcal{L}_{P}(\mathbf{x},\mathbf{y},\mathbf{\hat{y}})$
(6)

This approach shares similarities with two techniques, namely the introduction
of a regularization term and data augmentation.

The former limits the growth of weights during the training phase, thus
mitigating possible overfitting effects [34]. In this sense, the physical term
$\mathcal{L}_{P}$ acts as a regularization term.

Data augmentation [40] is a technique used to artificially increase the size
of the training dataset by applying various transformations to the input data.
The goal is to enhance the model’s performance by exposing it to a more
diverse set of examples, which helps improving generalization and robustness.
Our approach allows for an enrichment of the informational content of the
dataset by computing new physically-based quantities from input data
$\mathbf{x}$ and the predicted outputs $\mathbf{y}$. The $\mathcal{L}_{P}$
term compares such quantities with their true values, which depend on input
data $\mathbf{x}$ and the true outputs $\mathbf{\hat{y}}$. From this
perspective, this approach can be considered as a kind of physically-based
data augmentation and it allows for enriching the informational capacity of
the data.

As this approach involves only the loss function, it is a priori applicable to
a wide range of physical systems and compatible with all NN architectures.

###  2.3 Synthetic case definition

We construct a specific synthetic case study to test and show the proposed
methodology. We assess the performance of several approaches in the
reconstruction of the water surface profile along a 1D channel induced by the
presence of a weir placed at the outlet cross section. Such problem mimics a
commonly occurring scenario of determining the area affected by the presence
of an inline structure in a river. Additionally, if a supercritical flow
regimes develops in the upstream part of the channel, a hydraulic jump occurs
(Figure 1). The water profile encompassing a hydraulic jump does not belong to
the set of solutions of the differential equation generating the gradually-
varied water profile. Indeed, it is an internal boundary whose location (and
strength of the local discontinuity) needs to be solved for through additional
information (see B). Therefore, approaching the problem within a classical
PINN framework (namely the neural solver) would be non-trivial, since the jump
condition would require a tailored mathematical description.

The solution of the physical problem, assuming a cylindrical rectangular
channel and steady flow conditions, is a function $F^{*}$ which maps from
$\mathbb{R}^{6}$ to $\mathbb{R}$:

$\hat{h}=F^{*}(x,s,b,n,z_{d},Q)$ (7)

with $\hat{h}$ representing the true flow depth; $x$, the distance from the
dam; $s$, the channel slope; $b$, the channel width; $n$, the Manning
coefficient (related to the channel’s roughness); $z_{d}$, the height of the
dam; and $Q$, the water discharge.

The gradual varied flow assumption rules out the presence of discontinuities
in the depth profile. However, the transition from supercritical to
subcritical flow occurs through the appearance of a hydraulic jump which needs
to be accounted for as an internal boundary whose location and strength need
to be solved for through additional information.

Figure 1: Sketch of the analyzed physical problem, i.e. the steady state water
profile in a rectangular prismatic channel. Upper panel describes the
subcritical case; lower panel depicts the mixed regime case. For the meaning
of the symbols, the reader is referred to the text.

##  3 Numerical experiments

###  3.1 Generation of the training dataset

A set of 10390 profiles were generated by uniformly varying $s$, $b$, $n$,
$z_{d}$, and $Q$, solving equation 22 using a Finite-Difference scheme, with a
fixed spatial discretization of $\Delta x$ set to 10 meters, covering a total
length of 5000 meters. This results in profiles composed of 501 data points.

Therefore, we can interpret the problem as approximating $F^{*}$ in equation
(7) through a NN, having at our disposal a uniform sampling of the six-
dimensional function’s domain.

The obtained profiles were then randomly divided into training profiles (70 %,
totaling 7274 profiles), while the remaining ones constitute the validation
set and the test set (both consisting of 1558 profiles).

The input data to the various models are always normalized using the Standard
Scaler, which is a commonly used data pre-processing technique to scale the
features to have a mean of zero and a standard deviation of one. This is
necessary to prevent troublesome issues during training due to differences in
the numerical values of the features [2].

###  3.2 The three architectures

Three FFNN architectures were examined, namely the Single-Point (SP), the
INTegrator (INT), and the Vector-To-Sequence (VTS). We would like to point out
that the aim of this work is not to select the best architecture for solving
the problem at hand, but to evaluate the effects of a physically informed
training. Moreover, while it was also possible to employ Recurrent Neural
Networks, interpreting the temporal series output as the spatial series of the
water depth (namely the water profiles), we chose to limit the analysis to
FFNNs, aligning with the FFNN architecture of classical PINNs.

A sketch of the operating principles of SP, INT, and VTS architectures is
illustrated in Figure 2 and the corresponding network topologies are depicted
in Figure 3. Further details are in order.

Figure 2: Profile reconstruction using the three different approaches; input
and output are depicted as blue and orange circles, respectively; Single-Point
(upper panel) outputs the flow depth at a specific stationing; Integrator
(middle panel) requires the neighboring downstream value, regardless of the
stationing; Vector-To-Sequence outputs the entire vector of flow depth (i.e.
the whole water profile) at once.

Figure 3: The employed FFNN architectures. Left to right: Single-Point,
Integrator, Vector-To-Sequence.

####  3.2.1 Single-Point

A verbatim translation of the problem formulated in Eq. (7) is to employ the
NN to approximate the function $F^{*}$ with $F$:

$h=F(x,s,b,n,z_{d},Q;\bm{\theta})$ (8)

In this approach, the FFNN takes the six quantities governing the phenomenon
as inputs and predicts a single value for the water depth $h$ at the
stationing $x$ as an output.

To reconstruct the whole profile $\mathbf{h}$, a SP model needs to be run with
all the desired stationing values.

####  3.2.2 Integrator

In this method, a NN is utilized in the guise of a numerical integrator, that
is aimed at determining the local water depth based on its value at the
adjacent stationing. An eigenanalysis of the differential problem underpinning
a steady-state 1D free surface flow like the one chosen in this study, would
require to adopt either a downstream or upstream oriented solution direction
based on the local flow regime, namely super or subcritical, respectively.
However, in the context of surrogate models like the one based on ML, due to
the lack of any physical support, such requirement can be overlooked. In the
following, an upstream moving algorithm has been chosen. This methodology
requires a rearrangement of the dataset into the pairs
$\left([\hat{h}_{i},s,b,n,z_{d},Q],\hat{h}_{i+1}\right)$, where $\hat{h}_{i}$
and $\hat{h}_{i+1}$ must be uniformly spaced by $\Delta x$, where $i$ spans
all stationing points of all profiles.

Due to the above structure, the INT approach misses any information regarding
the distance from the dam. As a result, unlike the SP, it cannot predict the
height at any distance from the dam but only at multiples of the fixed $\Delta
x$ set by the dataset.

$h_{i+1}=INT(h_{i},s,b,n,z_{d},Q;\bm{\theta})$ (9)

Starting from a known downstream boundary condition ($h_{1}$), the integrator
$INT$ can be applied recursively. At each stationing, the output from the
downstream one serves as input, eventually leading to the reconstruction of
the whole profile.

$h_{i}=INT(h_{i-1},...)\circ INT(h_{i-2},...)\circ...\circ
INT(h_{1},s,b,n,z_{d},Q;\bm{\theta})$ (10)

It is essential to note that any classical numerical integrator would require
a repeated check for the occurrence of a hydraulic jump, as well as an ad hoc
procedure for its solution. Instead, this model can be applied flawlessly
across such discontinuity.

Two further specific features of the Integrator approach are in order.
Firstly, at each location, due to its recursive application, INT outputs a
depth value whose error depends on the accuracy of the previous applications.
As a consequence, the model’s accuracy generally decreases in the marching
direction of the algorithm. A side advantage of such upstream marching
algorithm lies in the possibility to incorporate physical conservation
balances between the current and previous location. Secondly, the recursive
application of INT introduces a subtle advantage over the other purely DD
approaches: indeed its initialization (occurring at the most downstream
stationing in our case) represents an implicit imposition of a physical
constraint.

####  3.2.3 Vector-To-Sequence

The Vector-To-Sequence employs a FFNN (namely $V$) that receives as input the
five parameters determining the profile solution and predicts as output the
entire vector $\mathbf{h}$.

$\mathbf{h}=V(s,b,n,z_{d},Q;\bm{\theta})$ (11)

The dataset must be rearranged in pairs
$\left([s,b,n,z_{d},Q],\hat{\mathbf{h}}\right)$. Just like for the INT
approach, the input does not include stationing data, thus implying that the
spacing of the output matches the one of the training dataset.

An advantage of this architecture is that it allows for the implementation of
physical loss terms whose formulation requires the knowledge of the whole
profile, e.g. the volume of water.

###  3.3 Physical training strategies

Each of the abovementioned three NN architectures has undergone both a purely
DD and a physically informed training, the latter consisting in exploiting the
local values of the:

  * 1.

Specific energy, as in equation (23), EN strategy in the following;

  * 2.

Froude number, as in equation (26), FR strategy in the following.

Since the VTS approach outputs the entire profile at once, it was also
possible to test three additional training strategies using:

  * 1.

The volume of water flowing in the river (namely, the area under the water
profile), VOL strategy in the following;

  * 2.

The downstream boundary condition, BC strategy in the following;

  * 3.

The energy differential equation, PDE strategy in the following.

The loss function for the purely DD training strategy is formulated in terms
of MSE between real and predicted water depth. For the physics-informed
training strategies, this function is augmented with an additional loss term,
as shown in Equation (6).

The physical loss terms adopted for each of the above list of training
strategies, read:

$\mathcal{L}_{P}^{EN}=\frac{\sum_{i=1}^{N}(E(h_{i})-E(\hat{h}_{i}))^{2}}{N}$
(12)
$\mathcal{L}_{P}^{FR}=\frac{\sum_{i=1}^{N}(Fr(h_{i})-Fr(\hat{h}_{i}))^{2}}{N}$
(13) $\mathcal{L}_{P}^{VOL}=\sum_{i=1}^{N}h_{i}-\sum_{i=1}^{N}\hat{h}_{i}$
(14) $\mathcal{L}_{P}^{BC}=|h_{1}-\hat{h}_{1}|$ (15)
$\mathcal{L}_{P}^{PDE}=\sum_{i=2}^{N-1}\left(\frac{E_{i+1}(h_{i})-E_{i-1}(h_{i})}{2\Delta
x}+i-J_{i}(h_{i})\right)$ (16)

The first four constraints above are what we refer to as soft physical
information. The last constraint, which is a hard physical information, is
only formally borrowed from the classical PINN paradigm: indeed it is here
employed in the context of a neural operator, i.e. aimed at improving its
training phase. In equations (12) - (13), $N$ is to be considered as the batch
size, that is the number of training examples utilised at each training
iteration, for the SP and INT architectures (which have point data), while for
VTS, $N$ in equations (12) - (16) represents the length of the data series.
Formally, to obtain the loss function at each iteration, for the VTS, the
average of the values obtained over the batch size must be considered (though
this is not shown here to avoid burdening the notation).

As introduced in Section 2.2, it is worth recalling that the $\mathcal{L}_{P}$
term acts as a regularization term, and the physical training strategies can
be interpreted as a form of physical data augmentation. Hereinafter, the term
”model” refers to the combination of an architecture and a training strategy.

###  3.4 Hyperparameters

Hyperparameters are values to be set before the training process and that are
not updated during the training phase. They encompass crucial features such as
the number of hidden layers and neurons, the optimization algorithm along with
its learning rate, and parameters related to early stopping. For the SP, we
initially assume 3 layers with 30 neurons each. Similarly, for the INT, we
also use 3 layers with 30 neurons. As for the VTS, we opt for 3 layers with 40
neurons each, considering the higher number of outputs the network needs to
produce compared to the other architectures. The ReLU activation function [33]
is consistently employed for each neuron. The number of layers is fixed for
all NNs as the analysis of its influence is beyond the scope of this work.

We implemented each neural network using TensorFlow [1] and Keras [12]. We
employ a learning rate reduction technique, specifically ReduceLROnPlateau,
within the Adam [26] optimization algorithm. The initial learning rate is set
to $0.001$ and is progressively reduced when approaching a minimum of the loss
function. We also employ an early stopping criterion during the training
phase, based on the MSE calculated on the validation set. All these
hyperparameters are fixed within each model to achieve consistent results.

However, it is important to recall here that the focus of this work is not to
evaluate the best model to solve the issue at hand, but rather the assessment
of the effects of the proposed methodology within each architecture. Thus any
influence of the topology that might favor one architecture over the other can
be overlooked. When adopting the physical training strategy, the value for the
hyperparameter $\lambda$ in equation (5) has been chosen as the one yielding
the best performance in the interval $[0,1]$.

##  4 Discussion of results

In order to acquaint the reader with the physical problem at hand, Figure 4
depicts a typical water level profile encompassing a hydraulic jump, as
predicted by the SP, the latter having been trained with and without the
inclusion of physical information. The reference solution, namely the profile
resulting from the Finite Difference integration, is depicted as well.
Incidentally, it is evident how the two physics-informed models outperform the
purely DD one. A quantitative analysis is in order.

Figure 4: Typical outcome of the comparison between the employed models,
namely DD, EN, FR, and the reference solution, FD: only the case of the SP
architecture is shown for the sake of clarity. The depicted profile,
determined by the quantities in the box, contains a hydraulic jump.

In this section we show the results obtained from the application of the three
different FFNN architectures, either in the purely DD or in the physics-
informed training strategies.

In the following, results yielded by the models in the initial configuration
are presented and analyzed. The training has been carried out on the entire
dataset. Then, the subsections contain results from several stress tests
carried out by:

i) reducing the size of the training dataset, while keeping the NN complexity
fixed, thus mimicking a ”small data regime”;

ii) varying the NN complexity in a both data-scarce and data-rich scenario,
therefore exploring both overfitting and underfitting conditions.

iii) applying models for extrapolation, i.e. seeking predictions beyond the
range of values covered during training, a scenario common to technical
applications.

In all the above tests, the effects of the physics-informed training strategy
are evaluated. It is indeed widely reported in the literature [17, 25, 41]
that the inclusion of physical information into ML models is beneficial
especially in conditions of small data regimes.

We employed two key metrics to evaluate the performance in reconstructing
water profiles, namely the Normalized Mean Absolute Error (NMAE) and the
Normalized Nash-Sutcliffe Efficiency (NNSE). The NMAE, assuming the dam’s
height as a representative length scale for the flow depth, is defined as:

$NMAE=\frac{\sum_{i=1}^{N}|h_{i}-\hat{h_{i}}|}{Nz_{d}}$ (17)

The NNSE is formulated as:

$NNSE=\frac{1}{(2-NSE)}$ (18)

where

$NSE=1-\frac{\sum_{i=1}^{N}(h_{i}-\hat{h}_{i})^{2}}{\sum_{i=1}^{N}(h_{i}-\bar{h})^{2}}$
(19)

and $\bar{h}$ is the average depth of the profile. The rationale behind the
scaling in Eq. (18) is to bound the NSE value between $0$ and $1$, thus
avoiding asymptotic tendency towards negative infinite values.

The NNSE is employed in addition to the NMAE due to its ability to amplify
errors made close to the location of the hydraulic jump. Indeed, NNSE
penalizes errors in areas where the flow depth is close to the mean value, and
the hydraulic jump is bounded by the two conjugate depths which are the
closest to the mean value within each profile.

A statistical analysis of the above two metrics shows a symmetrical shape of
the cumulative frequencies distributions attained over the test set, thus
allowing to employ their mean values as representative of the performance of
each approach.

###  4.1 Complete training dataset

The cumulative frequency distribution of NMAE and NNSE is employed to provide
a comprehensive picture of the model performance over the whole test set. In
Figure 5 we depict the distributions for the SP architecture and the tree
training strategies. Due to their definitions in Eqs. (17) and (18), a better
performing model features a higher cumulative frequency curve for NMAE and and
a lower one for NNSE metrics. The values of the error metrics associated to
the profiles shown in Figure 4 are reported as markers in Figure 5 for the
reader’s convenience.

Figure 5:  Cumulative frequency distributions of NMAE and NNSE values over the
test set. Markers represent the positioning of the results for the example
case shown in Figure 4

Figure 6: Box plot of the NMAE and NNSE distributions for SP (upper panels)
and INT (lower panels) architectures; boxes extend from the 25th to 75th
percentile; whiskers are placed at 10th and 90th percentiles; mean values and
median values are respectively shown as triangles and orange segments. An
overall improvement of predictive capabilities is detectable for both
architecture employing physical training strategies. Results for VTS are shown
in Figure 7.

Figure 7: Box plot of the NMAE and NNSE distributions for the VTS
architecture. Box plot features are the same as in Figure 6. VTS also yields
better predictions when trained with all types of the proposed physical
constraints, but the volume (VOL) one.

The overall picture suggests that the SP approach stands out with the best
performance, while the integrator displays the highest variance in its
results. The introduction of the physical information for all architectures
enhances their performance, with minimal differences between EN and FR
training strategies. In the context of VTS, all but the VOL training
strategies induce a slight improvement in the model’s performance (Figure 7).
Constraining the volume of the water profile seems to consistently provide a
deviant information content which detrimentally impacts the model’s
performance. Indeed, the same volume is shared by a large number of possible
output profiles, so that the inclusion of the physical loss term results in
the generation of additional local minima within the purely DD loss function.
This example unveils that, despite apparently providing additional and
physically sound information (i.e. the volume of water in the channel), the
training phase may be diverted towards non physical solutions. This peculiar
behavior could be shared with other physical constraints.

###  4.2 Stress tests: unbalance between training dataset size and model
complexity

We aim at training the above models so to induce either overfitting or
underfitting. The former occurs when the model achieves a high fit to the
training data while loses its ability to generalize to new unseen data. It can
be fostered by either reducing the size of the training set or increasing the
complexity of the model (i.e. the number of neurons within each layer). The
reduction in the number of parameters is carried out while maintaining the
three-layer structure and the same number of neurons among layers. The
underfitting, conversely, manifests as a poor performance of the model caused
by its lack of complexity with respect to the informational content of the
training data.

The data-scarce regime is obtained by excluding entire water profiles from the
training dataset, since two of the three models (INT and VTS) require the
stationing information to be retained in the reduced dataset.

In the following, in order to readily compare the performance of the models,
only the average values of the cumulative frequency distributions of the two
metrics are shown. Figure 8 shows that, while the performances progressively
worsen by reducing the training set as expected, the physical training
strategy improves the performance of all models and mitigates the overfitting,
thus effectively enhancing generalization capabilities. The data-scarce
scenario mimics a frequently occurring picture in the environmental hydraulics
framework, when the response of scarcely gauged basins is sought.

Figure 8: Stress test: mean values of the distributions of NMAE and NNSE
against the size of the training dataset. The physically based training
strategies consistently mitigate the loss of accuracy of predictions when
dataset size decreases.

Hereinafter, we vary the number of parameters in the models in both a data-
scarce and data-rich regimes. The latter employs the whole available training
dataset while the former only retains 5%.

In Figure 9, an expected trend emerges as the number of parameters decreases:
performance consistently deteriorates across all three architectures. This
discernible pattern clearly unveils the occurrence of underfitting. The lower
complexity architectures now struggle to capture the variability spanned by
the complete dataset. As a consequence, the models are insufficiently
expressive, resulting in suboptimal performance metrics. In such underfitting
regime, it is noteworthy that the introduction of a physical training strategy
does not yield a consistent improvement in generalization results. This
behavior is the consequence of the requirement for the low-complexity model,
which is already struggling to discern patterns in the data, to fulfill
additional constraints, such as the physical ones.

Instead, with reference to Figure 9, it is evident that when increasing model
complexity up to the occurrence of overfitting (on the right of the elbow of
the curve), the physics-informed models consistently yield better results than
the purely DD ones.

Figure 9: Stress test: mean values of the distributions of NMAE and NNSE
against models’ complexity (full training dataset). In the yellow-shaded
region (underfitting regime), physical constraints do not consistently lead to
improvements.

In the context of a small data regime, however, different behaviors emerge.
Figure 10 illustrates that, in situations with limited available data, model
performance shows loose dependence on the number of parameters and, in some
cases, it even deteriorates as the complexity of the models increases. Also in
this clear scenario of overfitting, models enhanced with the physical training
strategy consistently ensure an improvement in performance. The inclusion of
the physical loss term acts as a regularization mechanism, improving
generalization. Moreover, the physical term involves both input and output
data, akin to a form of data augmentation, providing the model with further
insights into the nature of the system to be interpreted.

Figure 10: Stress test: mean values of the distributions of NMAE and NNSE
against models’ complexity (5% of training data). In data scarce scenario, the
beneficial effect of adding the physical information emerges.

###  4.3 Extrapolation

Contrary to conventional practices in machine learning frameworks, we also
explore the extrapolation abilities of the models. New profiles were generated
by considering at least one of the values of the five inputs $s,b,n,z_{d},Q$
extending by 10% above or below the range covered by the training set. Out of
these profiles, 1558 were randomly selected (same size of test set), obtaining
a modified test set (referred to as EXT). The performance curves of the models
on the EXT set when decreasing the size of the training set (as previously
seen for the small data regime stress test), are illustrated in Figure 11.
These curves are compared with the results obtained using the standard test
set (interpolation). Beside the expected drop in overall performance compared
to the interpolation cases, the results confirm that models trained with
physical information exhibit greater generalization capabilities. These
results are of crucial importance for the perspective application to flood
mapping. Indeed, in such field predictions are often sought not only for
scenarios falling within the quantitative range of available observations,
though previously unexplored, but also for cases featuring values of the
observed quantities falling out of the range of the recorded series.

Figure 11: Stress test: assessment of models’ accuracy in predicting scenarios
out of the range of the training dataset (extrapolation, EX series, depicted
with light traces and markers); each panels show the same stress test results
for the models on scenarios encompassed by the trained dataset (same data
shown in 8, depicted with bold traces and markers). Though performance is
lower for extrapolated scenarios, the beneficial effect of physical
information is still evident.

##  5 Conclusions

In this work, we propose the incorporation of soft physical information in the
training strategy of neural operators to improve their generalization
capabilities in the context of environmental hydraulics.

Neural operators are designed to mimic the relationships between input and
output functions across an entire class of physical phenomena, as opposed to
neural solvers, which are aimed at finding the solution to a specific
realization. Indeed, the latter approach may not be useful in environmental
hydraulics, where a large number of scenarios need to be explored.

Moreover, the proposed soft physical information does not need to be
formulated in terms of the governing PDEs, providing a considerable advantage
in presence of large epistemic uncertainty, as in river hydraulics.

We tested the effect of soft physical information on a controllable and yet
highly informative test case commonly used in hydraulics, namely a 1D steady
state, mixed-regime, water profile in a rectangular channel. This physical
system often develops a discontinuity, represented by the hydraulic jump,
whose solution challenges the employed models.

Three NN architectures were trained on synthetic data alone and in conjunction
with soft physical information. The latter provide better predictive skills,
particularly in the presence of overfitting. Underfitting does not seem to be
mitigated by the added physical information. Moreover, in order to replicate a
common practice in the environmental modeling, the NNs predictions are
evaluated also for scenarios falling beyond the range of the training data
(extrapolation). Even in this case, the physical information allows for a
measurable gain in performance. These results are of great relevance to the
possible applications of NNs to flood mapping, where cases featuring values of
the observed quantities falling out of the range of the recorded series need
to be predicted. Finally, our analysis unveils the possible occurrence of a
detrimental effect of an apparently informative physical constraint. In this
sense, our work calls for the development of a systematic framework to measure
the informational content of physical constraints.

## Acknowledgments

We sincerely thank Dr. Andrea Gemma for his valuable support in setting up the
initial version of the code.

## Funding

The G. G. doctoral scholarship is granted under the National Recovery and
Resilience Plan (PNRR, NextGenerationEU).

J.-M. T. thanks the FRQNT “Fonds de recherche du Québec – Nature et
technologies (FRQNT)” for financial support (Research Scholarship No. 314328).

## Declaration of interests

None

##  Appendix A Neural Networks: Architectures and Training

Neural Networks were inspired by the architecture of nerve cells. A single
element of the network, called neuron, receives a series of $I$ inputs
$a_{i}$, calculates a weighted sum of these inputs (which may include a bias
value $c$), and then applies an activation function $\sigma$ to produce the
neuron’s output $z$.

$z=\sigma\left(\sum_{i=1}^{I}w_{i}a_{i}+c\right)$ (20)

Among the most commonly used activation functions are the Sigmoid, the
Hyperbolic Tangent, the Rectified Linear Unit (ReLU), and the Leaky ReLU [3]
function. The choice of the activation function affects the model outcome.

NNs can be categorised into Feed-Forward Neural Networks (FFNN), where the
signal propagates only from input to output, and Recurrent Neural Networks
(RNN), where the output of a neuron is fed back as its input. Both are
composed of an input layer, a series of $H$ hidden layers, and an output
layer. All networks employed in this work are FFNNs.

The approximate function $F$ represented by a FFNN can then be seen as the
composition of several functions $f$, one for each hidden layer:

$F=f_{H}\circ f_{H-1}\circ...\circ f_{1}$ (21)

A FFNN with more than one layer can theoretically represent any continuous, or
even discontinuous function, to arbitrary precision, given a sufficient number
of neurons.

The weights optimization phase (training) of NNs exploits a gradient descent
optimizer, such as the Adam algorithm [26], to iteratively find the minimum of
the loss function. The learning rate determines the step size taken during
each iteration of the optimization process and plays a crucial role to ensure
convergence. The gradient descent algorithm needs to evaluate the gradients
$\partial{\mathcal{L}}/\partial{w}$ of the loss function with respect to every
weight $w$ in the network. Such computation is carried out by means of the
Backpropagation algorithm [39], which efficiently calculates these gradients
by propagating the error backward through the network.

##  Appendix B Dataset generation

The dataset was generated by solving the specific energy equation, expressed
as [10]:

$\frac{dE}{dx}=i-J$ (22)

Here, $E$ represents the specific energy relative to the channel bottom,
which, for a rectangular channel, can be calculated as:

$E=h+\frac{Q^{2}}{2gb^{2}h^{2}}$ (23)

$g$ is the gravity, $b$ is the channel width, $h$ is the flow depth; $J$ is
the energy grade slope, calculated using the Chezy relation:

$J=\frac{n^{2}Q^{2}}{\Omega^{2}R^{4/3}}$ (24)

$n$ is the Manning coefficient, $\Omega=bh$ is the cross-sectional flow area
and $R=\Omega/(b+2h)$ is the hydraulic radius. The downstream boundary
condition is represented by the broad-crested weir equation:

$h=z_{d}+\left(\frac{3\sqrt{3}Q}{2\sqrt{2g}b}\right)^{2/3}$ (25)

where $z_{d}$ is the weir height.

The possible transition from supercritical ($Fr>1$) to subcritical ($Fr<1$)
flow occurs through a hydraulic jump, where:

$Fr=\frac{Q}{\Omega\sqrt{gh}}$ (26)

Examining the conservation of momentum within the fluid volume encompassing
the hydraulic jump, two cross-sections denoted as 1 (upstream, supercritical
flow) and 2 (downstream, subcritical flow) can be defined [10]. The
conservation of momentum can be expressed as:

$\Pi_{p1}+M_{1}=\Pi_{p2}+M_{2}$ (27)

Here, $\Pi_{p1}$ and $\Pi_{p2}$ represent the upstream and downstream values
of the hydrostatic force, while $M_{1}$ and $M_{2}$ denote the upstream and
downstream values of the momentum flux, respectively. For a rectangular
prismatic channel, equation (27) becomes:

$\frac{y_{1}^{2}}{2}+\frac{Q^{2}}{gb^{2}y_{1}}=\frac{y_{2}^{2}}{2}+\frac{Q^{2}}{gb^{2}y_{2}}$
(28)

which can be arranged in the following expression, establishing a mathematical
relation between $y_{1}$ and $y_{2}$, that are called conjugate depths:

$y_{1}=\frac{y_{2}}{2}(-1+\sqrt{1+8F_{2}^{2}})$ (29)

$F_{2}$ is the Froude number (equation (26)) at the cross-section 2.

The water profile is obtained by solving the differential equation (22) using
a first order FD scheme, with a constant spatial discretization $\Delta x$. An
upstream marching algorithm is employed, encompassing the solution of the
hydraulic jump if necessary.

## References

  * Abadi et al. [2015] Abadi, M., Agarwal, A., Barham, P., Brevdo, E., Chen, Z., Citro, C., Corrado, G.S., Davis, A., Dean, J., Devin, M., Ghemawat, S., Goodfellow, I., Harp, A., Irving, G., Isard, M., Jia, Y., Jozefowicz, R., Kaiser, L., Kudlur, M., Levenberg, J., Mané, D., Monga, R., Moore, S., Murray, D., Olah, C., Schuster, M., Shlens, J., Steiner, B., Sutskever, I., Talwar, K., Tucker, P., Vanhoucke, V., Vasudevan, V., Viégas, F., Vinyals, O., Warden, P., Wattenberg, M., Wicke, M., Yu, Y., Zheng, X., 2015.  TensorFlow: Large-scale machine learning on heterogeneous systems.  URL: https://www.tensorflow.org/. software available from tensorflow.org.
  * Ali et al. [2014] Ali, P.J.M., Faraj, R.H., Koya, E., Ali, P.J.M., Faraj, R.H., 2014.  Data normalization and standardization: a technical report.  Mach Learn Tech Rep 1, 1–6.
  * Apicella et al. [2021] Apicella, A., Donnarumma, F., Isgrò, F., Prevete, R., 2021.  A survey on modern trainable activation functions.  Neural Networks 138, 14–32.
  * Baydin et al. [2018] Baydin, A.G., Pearlmutter, B.A., Radul, A.A., Siskind, J.M., 2018.  Automatic differentiation in machine learning: a survey.  Journal of Marchine Learning Research 18, 1–43.
  * Bentivoglio et al. [2022] Bentivoglio, R., Isufi, E., Jonkman, S.N., Taormina, R., 2022.  Deep learning methods for flood mapping: a review of existing applications and future research directions.  Hydrology and Earth System Sciences 26, 4345–4378.
  * Berkhahn et al. [2019] Berkhahn, S., Fuchs, L., Neuweiler, I., 2019.  An ensemble neural network model for real-time prediction of urban floods.  Journal of Hydrology 575, 743–754.
  * Blöschl [2013] Blöschl, G., 2013.  Runoff prediction in ungauged basins: synthesis across processes, places and scales.  Cambridge University Press.
  * Cai et al. [2021] Cai, S., Mao, Z., Wang, Z., Yin, M., Karniadakis, G.E., 2021.  Physics-informed neural networks (pinns) for fluid mechanics: A review.  Acta Mechanica Sinica 37, 1727–1738.
  * Cedillo et al. [2022] Cedillo, S., Núñez, A.G., Sánchez-Cordero, E., Timbe, L., Samaniego, E., Alvarado, A., 2022.  Physics-informed neural network water surface predictability for 1d steady-state open channel cases with different flow types and complex bed profile shapes.  Advanced Modeling and Simulation in Engineering Sciences 9, 1–23.
  * Cengel and Cimbala [2013] Cengel, Y., Cimbala, J., 2013.  Ebook: Fluid mechanics fundamentals and applications (si units).  McGraw Hill.
  * [11] Chiu, P.H., Wong, J.C., Ooi, C.C., Dao, M.H., Ong, Y.S., .  Marrying the benefits of automatic and numerical differentiation in physics-informed neural network .
  * Chollet et al. [2015] Chollet, F., et al., 2015.  Keras.  https://keras.io.
  * Cole et al. [2006] Cole, S., Moore, R., Bell, V., Jones, D., 2006.  Issues in flood forecasting: ungauged basins, extreme floods and uncertainty, in: Frontiers in Flood Forecasting, 8th Kovacs Colloquium, UNESCO, Paris, pp. 103–122.
  * Cuomo et al. [2022] Cuomo, S., Di Cola, V.S., Giampaolo, F., Rozza, G., Raissi, M., Piccialli, F., 2022.  Scientific machine learning through physics–informed neural networks: Where we are and what’s next.  Journal of Scientific Computing 92, 88.
  * Dasgupta et al. [2023] Dasgupta, R., Das, S., Banerjee, G., Mazumdar, A., 2023.  Revisit hydrological modeling in ungauged catchments comparing regionalization, satellite observations, and machine learning approaches.  HydroResearch .
  * Di Baldassarre and Montanari [2009] Di Baldassarre, G., Montanari, A., 2009.  Uncertainty in river discharge observations: a quantitative analysis.  Hydrology and Earth System Sciences 13, 913–921.
  * Eichelsdörfer et al. [2021] Eichelsdörfer, J., Kaltenbach, S., Koutsourelakis, P.S., 2021.  Physics-enhanced neural networks in the small data regime.  arXiv preprint arXiv:2111.10329 .
  * Guo et al. [2021] Guo, Z., Leitao, J.P., Simões, N.E., Moosavi, V., 2021.  Data-driven flood emulation: Speeding up urban flood predictions by deep convolutional neural networks.  Journal of Flood Risk Management 14, e12684.
  * Hao et al. [2022] Hao, Z., Liu, S., Zhang, Y., Ying, C., Feng, Y., Su, H., Zhu, J., 2022.  Physics-informed machine learning: A survey on problems, methods and applications.  arXiv preprint arXiv:2211.08064 .
  * Hrachowitz et al. [2013] Hrachowitz, M., Savenije, H., Blöschl, G., McDonnell, J.J., Sivapalan, M., Pomeroy, J., Arheimer, B., Blume, T., Clark, M., Ehret, U., et al., 2013.  A decade of predictions in ungauged basins (pub)—a review.  Hydrological Sciences Journal 58, 1198–1255.
  * Jagtap et al. [2020] Jagtap, A.D., Kharazmi, E., Karniadakis, G.E., 2020.  Conservative physics-informed neural networks on discrete domains for conservation laws: Applications to forward and inverse problems.  Computer Methods in Applied Mechanics and Engineering 365, 113028.
  * Jamali et al. [2021] Jamali, B., Haghighat, E., Ignjatovic, A., Leitão, J.P., Deletic, A., 2021.  Machine learning for accelerating 2d flood models: Potential and challenges.  Hydrological Processes 35, e14064.
  * Jin et al. [2021] Jin, X., Cai, S., Li, H., Karniadakis, G.E., 2021.  Nsfnets (navier-stokes flow nets): Physics-informed neural networks for the incompressible navier-stokes equations.  Journal of Computational Physics 426, 109951.
  * Kabir et al. [2020] Kabir, S., Patidar, S., Xia, X., Liang, Q., Neal, J., Pender, G., 2020.  A deep convolutional neural network model for rapid prediction of fluvial flood inundation.  Journal of Hydrology 590, 125481.
  * Karniadakis et al. [2021] Karniadakis, G.E., Kevrekidis, I.G., Lu, L., Perdikaris, P., Wang, S., Yang, L., 2021.  Physics-informed machine learning.  Nature Reviews Physics 3, 422–440.
  * Kingma and Ba [2014] Kingma, D.P., Ba, J., 2014.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980 .
  * Kratzert et al. [2019] Kratzert, F., Klotz, D., Herrnegger, M., Sampson, A.K., Hochreiter, S., Nearing, G.S., 2019.  Toward improved predictions in ungauged basins: Exploiting the power of machine learning.  Water Resources Research 55, 11344–11354.
  * Kumar et al. [2023] Kumar, V., Sharma, K.V., Caloiero, T., Mehta, D.J., Singh, K., 2023.  Comprehensive overview of flood modeling approaches: A review of recent advances.  Hydrology 10, 141.
  * LeCun et al. [2015] LeCun, Y., Bengio, Y., Hinton, G., 2015.  Deep learning.  Nature 521, 436–444.
  * Löwe et al. [2021] Löwe, R., Böhm, J., Jensen, D.G., Leandro, J., Rasmussen, S.H., 2021.  U-flood–topographic deep learning for predicting urban pluvial flood water depth.  Journal of Hydrology 603, 126898.
  * Lu et al. [2021] Lu, L., Meng, X., Mao, Z., Karniadakis, G.E., 2021.  Deepxde: A deep learning library for solving differential equations.  SIAM review 63, 208–228.
  * Mahesh et al. [2022] Mahesh, R.B., Leandro, J., Lin, Q., 2022.  Physics informed neural network for spatial-temporal flood forecasting, in: Climate Change and Water Security: Select Proceedings of VCDRR 2021, Springer. pp. 77–91.
  * Nair and Hinton [2010] Nair, V., Hinton, G.E., 2010.  Rectified linear units improve restricted boltzmann machines, in: Proceedings of the 27th international conference on machine learning (ICML-10), pp. 807–814.
  * Ng [2004] Ng, A.Y., 2004.  Feature selection, l 1 vs. l 2 regularization, and rotational invariance, in: Proceedings of the twenty-first international conference on Machine learning, p. 78.
  * Nguyen et al. [2023] Nguyen, H.D., Dang, D.K., Nguyen, Y.N., Pham Van, C., Van Nguyen, T.T., Nguyen, Q.H., Nguyen, X.L., Pham, L.T., Pham, V.T., Bui, Q.T., 2023.  Integration of machine learning and hydrodynamic modeling to solve the extrapolation problem in flood depth estimation.  Journal of Water and Climate Change , jwc2023573.
  * Prestininzi et al. [2011] Prestininzi, P., Di Baldassarre, G., Schumann, G., Bates, P., 2011.  Selecting the appropriate hydraulic model structure using low-resolution satellite imagery.  Advances in Water Resources 34, 38–46.
  * Qian et al. [2019] Qian, K., Mohamed, A., Claudel, C., 2019.  Physics informed data driven model for flood prediction: application of deep learning in prediction of urban flood development.  arXiv preprint arXiv:1908.10312 .
  * Raissi et al. [2019] Raissi, M., Perdikaris, P., Karniadakis, G.E., 2019.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics 378, 686–707.
  * Rumelhart et al. [1986] Rumelhart, D.E., Hinton, G.E., Williams, R.J., 1986.  Learning representations by back-propagating errors.  Nature 323, 533–536.
  * Shorten and Khoshgoftaar [2019] Shorten, C., Khoshgoftaar, T.M., 2019.  A survey on image data augmentation for deep learning.  Journal of Big Data 6, 1–48.
  * Zhu et al. [2019] Zhu, Y., Zabaras, N., Koutsourelakis, P.S., Perdikaris, P., 2019.  Physics-constrained deep learning for high-dimensional surrogate modeling and uncertainty quantification without labeled data.  Journal of Computational Physics 394, 56–81.


# Spectrum Trading for Device-to-Device Communication In Cellular Networks
using Incomplete Information Bandwidth-Auction Game

Mohammad Karimzadeh Farshbafan  Dept. of Electrical and Computer Eng.
College of Eng., University of Tehran
Tehran, Iran, 14176-14418
<EMAIL_ADDRESS>Mohammad Hossein Bahonar  Dept. of Electrical and
Computer Eng.
Isfahan University of Technology
Isfahan, Iran, 84156-83111
<EMAIL_ADDRESS>Farshid Khajehrayeni  Dept. of Electrical and
Computer Eng.
Tarbiat Modares University
Tehran, Iran, 14115-111
<EMAIL_ADDRESS>
###### Abstract

Device-to-device (D2D) communication that allows proximity users to
communicate directly has been recently proposed to improve spectral efficiency
of cellular networks. In this paper, we assume a cellular network consisting
of multiple cellular user equipments (CUEs), which are the primary users, and
a cognitive D2D pair, which is the secondary user. The D2D pair needs a
bandwidth for data transmission that can be obtained via spectrum trading. We
introduce a bandwidth-auction game for the spectrum trading problem. The base
station (BS) and CUEs are able to sell their spectrum or share it with the D2D
pair, which allows the D2D pair to operate in orthogonal sharing, cellular, or
non-orthogonal sharing (NOS) modes. Operation of the D2D pair in the NOS mode
causes interference to the CUEs, which is possible under low interference
condition. In the auction, the D2D pair can buy its required spectrum from
three different service providers (SPs) corresponding to each mode that
operate on different frequency spectrums. The D2D pair bids a price-bandwidth
demand curve and the SPs offer a price-demand supply curve. Since each player
is not aware of the strategy of other players in practical scenarios, the game
is assumed to be an incomplete information repeated one. A best response based
learning method is proposed for the decision making procedure of all players,
the D2D pair and SPs. It is shown that the proposed method converges to the
Nash equilibrium (NE) point of the game more rapidly than the state-of-the-art
methods when the game is played repeatedly. The sensitivity of the proposed
method to the learning rate variable is also less than the state-of-the-art
methods and hence can be considered as a robust one.

###### Index Terms:

Game Theory, Auction Theory, Device-to-Device Communication, Incomplete
Information, Repeated Game.

##  I Introduction

Cognitive radio is a promising techniques that has been proposed in order to
solve the spectrum scarcity problem [1]. In a cognitive radio network,
secondary users (SUs) are able to reuse the radio resources of primary users
(PUs) while the resources are not being used by PUs or the interference caused
to PUs from SUs does not degrade the quality of service of PUs below a
predetermined threshold. In addition, device-to-device (D2D) communication
which corresponds to direct data transmission between two proximity user
equipments (UEs) without going through the base station (BS) has also been
proposed in order to increase the overall spectral efficiency (SE) of cellular
networks or to offload part of the cellular traffic under heavy load of
cellular networks [2]. This communication technique that enables the network
to provide service to a larger number of users has been considered in 3GPP
LTE-Advanced standard.

Due to the potentials of D2D communication, it has received much attention
recently [3, 4, 5, 6]. It is also possible to use D2D communication for
cognitive radio systems. Some researches and technologies such as [6, 7] are
based on the usage of cognitive radio concepts for D2D communication in
cellular networks. D2D communication in the licensed spectrum can be used in
an overlay manner where the resources of D2D pairs are different from that of
cellular user equipments (CUEs) or an underlay manner where D2D pairs and CUEs
use the same set of resources. The resource allocation problem for D2D
communication of D2D pairs has been studied in centralized schemes using
optimization techniques [3] or distributed schemes using game theoretic models
[8]. Since game theoretic based methods are designed in a distributed manner,
they are usually scalable and less complex [9]. Different types of games such
as coalitional game [10], dynamic stackelberg game [11], and auction-based
games [8] has been used to model and solve the problems that are related to
D2D communication in cellular networks.

Auction theory has been employed to model and investigate the spectrum sharing
problem in D2D communication [8] and cognitive radio systems [12]. In [12],
the spectrum access problem of SUs is modeled using a repeated auction game
where a Bayesian learning algorithm is proposed as the solution that lead to
the Nash equilibrium (NE) point of the game. Auction games can also be used
for spectrum trading or spectrum sharing problem in D2D communication. In
[13], a reverse iterative combinatorial auction game has been proposed as the
resource allocation method of D2D communication, which leads to a good
performance of system sum-rate. In [14] a Cournot game is proposed to model
the problem of competitive pricing where a few number of firms compete with
each other in terms of price in order to gain the highest profit. However, it
is assumed that the price is completely in hands of PUs in this model, and it
does not consider the demand curve of the SUs. In [15] a double sided auction
mechanism is proposed to solve the dynamic spectrum sharing problem. The game
considers the supply and bid strategy for suppliers and bidders, respectively.
Based on the proposed model, the bidders and suppliers do not have any
contribution to the amount of bandwidth and price, respectively. In [7] the
authors propose a bandwidth-auction game for resource allocation to one D2D
pair in a cellular network and derive the NE point analytically. Since the
assumption of complete information is not a practical one, a gradient based
learning algorithm is also proposed, which converges to the NE point of the
game.

In this paper, we assume a cellular network consisting of a base station (BS)
and multiple CUEs that are considered as the PUs. A D2D pair, which is
considered as the SU, exists in the cell and needs spectrum resources for data
transmission. The spectrum sharing problem is formulated as a bandwidth-
auction game. The D2D pair can buy its required bandwidth from three types of
service providers (SPs) that are willing to share or resale their spectrum.
The BS and CUEs are able to provide bandwidth for the D2D pair. The D2D pair
can use the spectrum of the SPs that are willing to resale their resource in
an overlay manner where all UEs use orthogonal resources. In this way, the D2D
communication does not cause interference to the other UEs. The D2D pair is
also able to reuse the spectrum of the CUEs in an underlay manner under low
interference conditions. Therefore, CUEs share their resource with the D2D
pair, and D2D communication causes interference to the other CUEs. Hence, the
BS and CUEs are the suppliers of the game while the D2D pair is the consumer
and the BS acts as the brooker in the game. Each palyer aims to maximize its
own payoff. Our system model is mostly based on the system model of [7]. The
D2D pair bids a demand curve and each supplier offers a supply curve. In order
to have a practical and more realistic system model, it is assumed that each
player is not aware of the strategy of the other player and can just be aware
of the bandwidth price. Thus, the game is modeled as a incomplete information
game. In [7], it is proposed to use a gradient base learning method for the
incomplete information. We propose to utilize a best response based learning
method in a repeated manner as the solution of the incomplete information
game. It is shown that the proposed algorithm converges to the NE point of the
game iteratively. The proposed method converges rapidly and quicker than the
state-of-the-art proposed method of [7].

The rest of the paper is organized as follows. In Section II, the system model
is described and the bandwidth-auction game is presented. In Section III the
NE point of the game is calculated and the learning method for the incomplete
game is proposed. The simulation results are presented in Section IV, and
Section V concludes the paper.

##  II System Model

###  II-A Cellular Network Model

In this study, a multi cell network is assumed in which each UE can function
in either cellular or D2D modes. In fact, they communicate with each other
through the BS in the cellular mode. Moreover, they are allowed to work in a
D2D mode on certain conditions. First, the UE must be placed in close
proximity in order to alleviate the incoming interference to the other CUEs,
regardless of being in a cell or not. Second, there should be low amount of
interference on the receiver of the D2D pair, so that its corresponding SINR
remains above a predefined minimum threshold. Third, the network must suffer
from heavy load and prefer to offload a portion of its traffic load. The UE
that satisfies the aforementioned conditions are able to have a D2D
communication. We call the UEs of this D2D pair as $U_{1}^{d}$ and $U_{2}^{d}$
which are the transmitter and receiver of the D2D pair, respectively. Then, a
request is sent to the BS and it allocates the required bandwidth.

Three various sharing options can be considered including non-orthogonal
sharing (NOS), orthogonal sharing (OS), and cellular mode (CM). In the case of
NOS mode, the same bandwidth with the size of $s_{1}$, is reused by the D2D
pair, resulting in experiencing interferences at the other CUEs caused by the
D2D pair and at the receiver side of the D2D pair caused by CUEs. In the case
of OS mode, one or more CUEs can resell their spectrum by providing a portion
of their subcarriers to the D2D pair. Hence, there will be no interference
caused to the D2D pair by the CUE and also to the CUE caused by the D2D pair.
In this case, the CUEs are denoted by $U_{n}^{c},n=2,...,N+1$ where $N$ is the
number of CUEs who are willing to share their spectrum and the corresponding
size of the allocated spectrum is denoted by $s_{j},j=2,...,N+1$. In the case
of CM, the operator itself acts as a SP and sells its spare spectrum
($s_{N+2}$). Thus, BS performs as a relay node and allocates uplink and
downlink bandwidth (denoted by $B_{UL}$ and $B_{DL}$). It should be noted that
the communication between users is based on orthogonal frequency division
multiple access (OFDMA), which is consistent with LTE advanced network. Thus,
in each mode, one or more SPs attempt to sell their spectrum to the D2D
transmitter at price $\alpha_{j}^{S}$.

The game is modeled as a bandwidth-auction game where there exists three kinds
of SPs. The operator, CUEs, and the receiver of the D2D pair in NOS mode act
as SPs which are the suppliers of the game. The D2D pair is the bidder of the
game. Each supplier offers a supply curve and the bidder bids a demand curve.
There are two factors that affect the spectrum demand of the D2D transmitter,
namely SE and the price. The SE of $U_{1}^{d}$ depends on the chosen source
(SP). In the NOS mode, Let ${\xi_{NOS}}$ be the Signal-to-Interference-plus-
Noise-Ratio (SINR) of the D2D pair, $P_{NOS}$ denote the transmit power of
$U_{1}^{d}$, $BER_{tar}$ be the target bit-error-rate, and $g_{NOS}$ indicate
the channel gain of the D2D pair. Then, the SE (throughput) of $U_{1}^{d}$ can
be written as [1]:

$\displaystyle T_{1}$
$\displaystyle=\log_{2}(1+\Theta\xi_{NOS}),\leavevmode\nobreak\ \rm{where}$
(1) $\displaystyle\Theta$
$\displaystyle=\frac{1.5}{\ln(0.2\times(BER_{tar})^{-1})},\leavevmode\nobreak\
\xi_{NOS}=\frac{P_{NOS}g_{NOS}}{\sigma^{2}+I},$

in which $\sigma^{2}$ and $I$ represent the power of additive white Gaussian
noise (AWGN) and interference on the receiver of the D2D pair, respectively.
It must be mentioned that static channels with path-loss model, rayleigh fast
fading, and log-normal slow fading are assumed in this mode. Channel gain
between the receiver and transmitter of D2D pair is denoted by
$g_{NOS}={(K/d_{1,2})}^{p}\psi_{1,2}\gamma_{1,2}$ where $K$ is a constant
depending on the working frequency, $d_{1,2}$ denotes the distance between the
D2D transmitter and receiver, and $p$ (normally between 2 and 6) represents
the path-loss parameter, $\psi_{1,2}$ is the rayleigh fast fading gain of the
link, and $\gamma_{1,2}$ denotes its corresponding log-normal slow fading
gain.

In the OS mode, it is supposed that $\xi_{OS}$ indicates the SNR at the CUE
and $P_{OS}$ is the transmit power of $U_{1}^{d}$. Then, the throughput of the
$U_{1}^{d}$ can be expressed as:

$\displaystyle T_{j}$
$\displaystyle=\log_{2}{(1+\Theta\xi_{OS})},\;\rm{where}$ (2)
$\displaystyle\xi_{OS}$
$\displaystyle=\frac{P_{OS}g_{OS}}{\sigma^{2}},\;\;\text{for}\;\;j=2,3,...,N+1.$

In the CM mode, let $\xi_{CM,1}$ and $\xi_{CM,2}$ be the SNR at the BS and
that of at the receiver of the D2D pair, respectively. Furthermore, we denote
the transmission power emitted from $U_{1}^{d}$ in this mode as $P_{CM}$, the
channel gain between the D2D transmitter and the BS as $g_{CM,1}$, and that of
between the BS and the D2D receiver as $g_{CM,2}$. Then, the throughput of
$U_{1}^{d}$ can be formulated as [1]:

$T_{N+2}=\min\left\\{\begin{array}[]{ll}\frac{B_{UL}}{B_{UL}+B_{DL}}\log_{2}(1+\Theta\xi_{CM,1})\\\
\frac{B_{UL}}{B_{UL}+B_{DL}}{\log_{2}(1+\Theta\xi_{CM,2})}\end{array}\right\\},$
(3)
$\displaystyle\xi_{CM,1}=\frac{P_{CM}g_{CM,1}}{\sigma^{2}},\xi_{CM,2}=\frac{P_{CM}g_{CM,2}}{\sigma^{2}}.$

If we assume the distance between the D2D pair is long enough, the following
formula for the throughput can be derived [1]:

$\displaystyle T_{N+2}$
$\displaystyle\approx\frac{B_{UL}}{B_{UL}+B_{DL}}\log_{2}(1+\Theta\xi_{N+2}),$
(4) $\displaystyle\xi_{N+2}$ $\displaystyle=\frac{P_{CM}g_{CM}}{\sigma^{2}},$

where $g_{CM}$ is the channel gain between the D2D pair and the BS.

As it is previously mentioned, we assume a situation in which there are $N+2$
SPs ($U_{1}^{d}$, $N$ cellular users, and BS) and only one bidder. Each of
these players offers a price based on their strategy functions (supply
functions for SPs and demand functions for the D2D transmitter), and the
players’ desire is to maximize their own payoff. This is a uniform price
auction game whose aim is to set each player’s strategy in a way that they
reach the NE point of the game. In other words, both demand and supply
functions are known to the broker, which is the BS in our model, and its duty
is to calculate the ultimate price by adjusting the bandwidth demand and
supply. Having calculated the price, the broker have a duty to allocate the
determined bandwidth to the D2D user.

###  II-B Bandwidth Acution Game Model

As it was mentioned before, we consider a D2D pair demanding an amount of
bandwidth and $(N+2)$ SPs that compete for providing the requested bandwidth
to the D2D pair. First of all, a demand function for the D2D pair and supply
functions of the SPs are defined. The considered demand and supply functions
for the players of the bandwidth auction game are as follows [1]:

$\displaystyle p$ $\displaystyle=-\Lambda^{D}b^{D}+\lambda^{D},$
$\displaystyle p$
$\displaystyle=\Lambda^{S}_{i}b^{S}_{i}+\lambda^{S}_{i},\;\;\;i=1,2,\ldots,N+2.$
(5)

where $b^{D}$ is the demand bandwidth of the D2D pair, $b^{S}_{i}$ is the
supply bandwidth of the $i^{\text{th}}$ SP and $p$ is the clearing price for
the unit of the bandwidth. Also, $\Lambda^{D}$, $\Lambda^{S}_{i}$,
$\lambda^{D}$ and $\lambda^{S}_{i}$ are the strategies of different players
which are positive numbers. The considered functions in (5) means that D2D
user buys more bandwidth at the lower price and vice-versa. Also, SPs sell
more bandwidth at the higher price and vice-versa.

The payoff function for the players of the bandwidth auction game must be
defined. For the D2D player, it can be defined as [1]:

$\displaystyle\pi^{D}=U(b^{D})-pb^{D},\quad
U(b^{D})=\sum_{i=1}^{N+2}{R_{i}T_{i}b^{S}_{i}},$ (6)

where $U(b^{D})$ is the utility function of the D2D player, $R_{i}$ is its
revenue when using the bandwidth of the $i^{\text{th}}$ SP and $T_{i}$ is the
throughput of the D2D pair when using the bandwidth of the $i^{\text{th}}$ SP
which is introduced in (1)-(3).

The payoff function for the SPs can also be written as

$\displaystyle\pi^{S}_{i}=pb^{S}_{i}-C_{i}(b^{S}_{i}),\quad
C_{i}(b^{S}_{i})=c_{i}b^{S}_{i},$ (7)

where $C_{i}(b^{S}_{i})$ is the cost of providing bandwidth for the
$i^{\text{th}}$ SP. As seen in (7), the $C_{i}(b^{S}_{i})$ for all of the SPs
is considered linearly proportional to the sold bandwidth, $b^{S}_{i}$.
However, the unit cost of providing bandwidth, $c_{i}$, is determined
differently for the SPs. More precisely, for the first SP, the unit cost of
providing bandwidth can be written as $c_{1}=c_{1}^{BS}+c_{1}^{I}$ in which
$c_{1}^{BS}$ is the unit price of the bandwidth that first SP has paid to the
BS and $c_{1}^{I}$ is the unit cost of the imposed interference to the first
SP. As mentioned earlier, the second to the $(N+1)^{\text{th}}$ SPs are the
cellular users. As a result, the $c_{i},\;i=2,3,\ldots,N+1$ is the price of
purchasing bandwidth from the BS for the $(i-1)^{\text{th}}$ cellular user.
Finally, for the last SP which is the BS, $c_{N+2}$ is the price of providing
bandwidth.

##  III Proposed Method for Bandwidth Auction Game

In this section the proposed method is described. We propose a learning method
for the incomplete information model of the game which is more reasonable and
practical in section III-B. In order to evaluate the performance of the
proposed method, the NE point of the complete information game is also
calculated in section III-A.

###  III-A Complete Information Game

In a complete information game, all the players are aware of the strategies,
payoffs, and demand/supply functions of the other players. Using the demand
and supply functions in (5) and also the fact that
$b^{D}=\sum_{i=1}^{N+2}{b_{i}^{S}}$, the clearing price can be written as

$\displaystyle
p=\dfrac{\dfrac{\lambda^{D}}{\Lambda^{D}}+\sum_{i=1}^{N+2}{\dfrac{\lambda^{S}_{i}}{\Lambda^{S}_{i}}}}{\dfrac{1}{\Lambda^{D}}+\sum_{i=1}^{N+2}{\dfrac{1}{\Lambda^{S}_{i}}}}.$
(8)

When the price is cleared using (8), each player can determine the amount of
purchased or sold bandwidth using (5).

It is known that all of the players are willing to maximize their payoff
function which is a function of the clearing price. The players can maximize
their payoff by adjusting their strategy
$\big{(}\Lambda^{D},\lambda^{D},\Lambda^{S}_{i},\lambda^{S}_{i}\big{)}$. On
the other hand, the unit price of bandwidth, $p$ is a function of all players’
strategies. As a result, the strategy of each player will affect the payoff of
the other players. The problem with this characteristic can be modeled as a
game.

As mentioned earlier, the strategy of each player is to adjust the
coefficients of the demand/supply functions which are introduced in (5).
However, when both of the $\Lambda^{D},\lambda^{D}$ for the D2D user and
$\Lambda^{S}_{i},\lambda^{S}_{i}$ for the SPs are considered as the strategies
of the players, the mathematical approach of NE calculation becomes
intractable. As a result, only one of the coefficients on the demand/supply
functions should be selected as the strategy of the players. The authors in
[8] proved that if the $\Lambda^{D}$ for the D2D user and $\Lambda^{S}_{i}$
for the SPs are selected as the strategy of the players and $\lambda^{D}$ and
$\lambda^{S}_{i}$ are considered constant and equal to their optimal values,
the NE point does not exist. Hence, $\lambda^{D}$ and $\lambda^{S}_{i}$ are
considered as the strategy of the players. Moreover, it is assumed that the
$\Lambda^{D}$ and $\Lambda^{S}_{i}$ are positive and have constant values.

Therefore, the NE point of the game while $\lambda^{D}$ and $\lambda^{S}_{i}$
are considered as the strategy of the players and $\Lambda^{D}$ and
$\Lambda^{S}_{i}$ as positive constants can be calculated. The NE point is the
intersection point of the best response of all players. Considering this
method and based on the solution presented at [8], the NE point can be
calculated as follows

$\displaystyle\bm{\mathrm{\lambda}}^{*}=\bm{\mathrm{A}}^{-1}\bm{\mathrm{B}},$
(9)

in which
$\bm{\mathrm{\lambda}}=\left[{\lambda_{1}^{s},\ldots,\lambda_{N+2}^{s},\lambda^{D}}\right]^{T}\in\Re^{N+3}$
is the ultimate price related to the NE point. ordered lexicographically as a
column vector,
$\bm{\mathrm{A}}=\left[{a_{i,j}}\right]\in\Re^{(N+3)\times(N+3)}$ is an
invertible square matrix ($det(\bm{\mathrm{A}})\neq 0$), and
$\bm{\mathrm{B}}=\left[{B_{1},\ldots,B_{N+2},B_{N+3}}\right]^{T}\in\Re^{N+3}$
is a column vector. The following equation illustrates how elements of the
mentioned matrices are calculated:

$\displaystyle a_{i,i}$ $\displaystyle=2(\Phi\Lambda_{i}^{s})^{-1}-2,$
$\displaystyle a_{i,j}$
$\displaystyle=2(\Phi\Lambda_{j}^{s})^{-1}-\Lambda_{i}^{s}(\Lambda_{j}^{s})^{-1},\quad
i\neq j,$ $\displaystyle a_{i,{N+3}}$
$\displaystyle=2(\Phi\Lambda^{D})^{-1}-\Lambda_{i}^{s}(\Lambda^{D})^{-1},$
$\displaystyle a_{{N+3},j}$
$\displaystyle=2(\Phi\Lambda_{j}^{s})^{-1}-\Lambda^{D},$ $\displaystyle
a_{{N+3},{N+3}}$ $\displaystyle=2(\Phi\Lambda^{D})^{-1}-2,$ $\displaystyle
B_{i}$
$\displaystyle=c_{i}(1-\Phi\Lambda_{i}^{s}),\quad\Phi=\sum\limits_{n=1}^{N+2}({\Lambda_{n}^{s}})^{-1}+(\Lambda^{D})^{-1},$
$\displaystyle B_{N+3}$
$\displaystyle=-\Lambda^{D}\sum\limits_{n=1}^{N+2}R_{n}T_{n}(\Lambda_{n}^{s})^{-1},$
$\displaystyle i,j$ $\displaystyle=1,2,\ldots,N+2.$

It should be noted that if matrix $\bm{\mathrm{A}}$ is not an invertible
matrix ($det(\bm{\mathrm{A}}=0$), there will be no NE point defined for the
game.

The most important assumption of using (9) for calculating the NE point is the
complete information assumption of the game. More precisely, for calculation
of the NE point using (9), each player should be informed about the payoff
function and demand/supply functions of the other players which is not
possible in a practical scenario.

###  III-B Incomplete Information Game

The previous sections derive the NE point of the game considering the complete
information assumption which is not a reasonable assumption in a practical
scenario. Practically, each player cannot be aware of the payoff,
demand/supply, or strategy of other players. In fact, the bandwidth
information is the only information each player has. According to these
explanations, the players only can reach the NE point of the game in multiple
steps. More precisely, they should learn to improve their strategy until
reaching the NE point. In [8], a learning method for adjusting the strategy of
the players based on the marginal profit function is introduced. In this
method, the players make their strategy in each step based on a local
estimation of their marginal payoffs. This adjustment mechanism is called
myopic. We propose the best response based learning method for the incomplete
information game model.

In the best response based learning, each player utilizes the best response
function for updating the playing strategy in each step. As a result, the
strategy of each player of the considered game can be written as follows

$\displaystyle\lambda^{D}(t+1)$
$\displaystyle=f^{D}\big{(}\lambda^{S}_{1}(t),\ldots,\lambda^{S}_{N+2}(t)\big{)},$
(10) $\displaystyle\lambda^{S}_{i}(t+1)$ $\displaystyle=$ (11) $\displaystyle
f_{i}^{S}\big{(}\lambda^{D}(t),$
$\displaystyle\lambda^{S}_{1}(t),\ldots,\lambda^{S}_{i-1}(t),\lambda^{S}_{i+1}(t),\ldots,\lambda^{S}_{N+2}(t)\big{)},$

where $f^{D}$ and $f_{i}^{S}$ are the best response function of the D2D user
and the $i^{\text{th}}$ SP, respectively. This updating approach might lead to
large fluctuations in the strategy of each player due to not using the
strategy of the player in the previous step. As a result, converging to the
Nash equilibrium point cannot be guaranteed. Owing to these challenges, the
updating formula in (10)-(11) are revised as follows

$\displaystyle\lambda^{D}(t+1)$
$\displaystyle=(1-\tau^{D})\lambda^{D}(t)+\tau^{D}f^{D}\big{(}\lambda^{S}_{1}(t),\ldots,\lambda^{S}_{N+2}(t)\big{)}$
(12) $\displaystyle\lambda^{S}_{i}(t+1)$
$\displaystyle=(1-\tau^{S}_{i})\lambda^{S}_{i}(t)+$ (13)
$\displaystyle\tau^{S}_{i}f_{i}^{S}\big{(}\lambda^{D}$
$\displaystyle(t),\lambda^{S}_{1}(t),\ldots,\lambda^{S}_{i-1}(t),\lambda^{S}_{i+1}(t),\ldots,\lambda^{S}_{N+2}(t)\big{)},$

where $0<\tau^{D}\leq 1$ and $0<\tau^{S}_{i}\leq 1$ are the learning rates of
the D2D user and the $i^{\text{th}}$ SP, respectively. For determining the
best response function for the D2D user and SPs, the clearing price of the
last step is used. As a result, the updating formula in (12)-(13) can be
written as

$\displaystyle\lambda^{D}(t+1)=$ $\displaystyle(1-\tau^{D})\lambda^{D}(t)+$
(14)
$\displaystyle\tau^{D}\Big{(}p(t)(2-K(t)\Lambda^{D})+\Lambda^{D}C(t)\Big{)},$
$\displaystyle\lambda^{S}_{i}(t+1)=$
$\displaystyle(1-\tau^{S}_{i})\lambda^{S}_{i}(t)+$ (15)
$\displaystyle\tau^{S}_{i}\Big{(}p(t)+(p(t)-\Lambda_{i}^{S})(1-K(t)\Lambda_{i}^{S})\Big{)},$

where $p(t)$ is the last clearing price. Also, $C(t)$ and $K(t)$, which are
unknown parameters to the players, can formulated as

$\displaystyle K(t)\simeq$
$\displaystyle\dfrac{2p(t)+\Lambda^{D}C(t)}{p(t)\Lambda^{D}+(\Lambda^{D})^{2}\phi^{D}(t,t-1)}$
(16) $\displaystyle\phi^{D}(t,t-1)=$
$\displaystyle\dfrac{\pi_{D}(t)-\pi_{D}(t-1)}{\Lambda_{D}(t)-\Lambda_{D}(t-1)}$
(17) $\displaystyle C(t)=$
$\displaystyle\sum_{i=1}^{N+2}{\dfrac{R_{i}T_{i}}{\Lambda_{i}^{S}}}$ (18)

As it can be seen, both of these unknown parameters depend on the strategies
of all players. As a result, each player should estimate them through the
learning procedure. However, as described in (16), $K(t)$ is a function of
$C(t)$. Hence, it is assumed that the value of $K(t)$ is updated in odd time
step using (16) and $C(t)=C(t-1)$. Also, the value of $C(t)$ is updated in
even time step using (19) and $K(t)=K(t-1)$.

$\displaystyle C(t)=$
$\displaystyle\dfrac{K(t)\Big{(}p(t)\Lambda^{D}+(\Lambda^{D})^{2}\phi^{D}(t,t-1)-2p(t)\Big{)}}{\Lambda^{D}}$
(19)

##  IV Simulation Results

In this section, we evaluate the performance of the proposed learning method
using numerical results. We consider a single cell cognitive radio-based
cellular system with one BS, one D2D pair, and one CUE where the BS is placed
at the middle point of the cell. CUEs are uniformly distributed in the cell
using 2D uniform distribution. Due to the short distance between transmitter
and receiver of a D2D pair, the transmitter of D2D pair is uniformly
distributed in the cell while the receiver of the D2D pair is placed uniformly
in a cluster around the transmitter.

To draw a comparison among state-of-the-art methods in the incomplete
information case, we employed the presented simulation setup in [8] that used
a gradient-based learning method. The parameters of SPs and the D2D pair are
presented in (20).

$\displaystyle\Lambda^{S}_{1}=0.6,\Lambda^{S}_{2}=0.45,\leavevmode\nobreak\
\Lambda^{S}_{3}=0.6,\rm{and}\Lambda^{D}=0.5$ (20)

Also the target BER for the users is assumed to be $BER_{tar}=10^{-4}$.
Furthermore, the channel qualities of these SPs in NOS mode, OS mode, and CM
which are represented in terms of SINR are introduced in (21).

$\xi_{NOS}=20\rm{dB},\leavevmode\nobreak\
\xi_{OS}=30\rm{dB},\leavevmode\nobreak\ \rm{and}\leavevmode\nobreak\
\xi_{CM}=10\rm{dB}$ (21)

The channel quality of CM is lower than that of other modes due to the large
distance between a user and the BS. The channel quality of NOS mode is also
lower than that of OS mode due to low amount of interference. Considering the
described simulation parameters, there can be a fair comparison between the
our simulation results and those of [8].

Simulation results verify the convergence of the players’ strategies to the NE
point of the game after sufficient number of time steps. After a sufficient
number of time steps, the strategies of $U_{2}^{d}$, $U^{c}$, BS, and the D2D
user converge to $(4.8656,4.3152,7.6923,24.9058)$. By calculating the NE point
of the game in the complete information case using equation (9), it can be
shown that the NE point of the game $\bm{\mathrm{\lambda}}^{*}$ is equal to
$\bm{\mathrm{\lambda}}^{*}=(4.8656,4.3152,7.6923,24.9058)$. Hence the best
response based proposed method converges to the NE point of the game.

In the following, we evaluate the performance of two learning method including
gradient-based learning and the proposed best response based learning
regarding their convergence to the NE point in the incomplete information
case. In Fig. 1, the strategies of the SPs and the D2D user during the
learning game, using best response method and convergence to the NE point are
indicated. The learning rates of all players are set to $0.6$. As seen in this
figure, the strategies of all players converge to the NEm point of complete
information case. It is worth noting that the convergence occurs only in $5$
time steps. However, as indicated in [8], the convergence of strategies using
gradient-based learning method occurs in $9$ time steps. As a result, the
convergence speed of the proposed algorithm is higher than gradient-based
learning method.

Figure 1: Strategies of the SPs and D2D user during the learning game and
convergence to the NE point when learning rates are set to $0.6$. Figure 2:
Strategies of the SPs and D2D user during the learning game and convergence to
the NE point of the game when learning rates of the players are set to $0.4$.
Figure 3: Strategies of the SPs and D2D user during the learning game and
convergence to the NE point of the game when learning rates of the players are
set to $0.1$.

The effect of different learning rates on the convergence of the players’
strategy is considered in Fig. 2 and Fig. 3. More precisely, Fig. 2 and Fig. 3
demonstrate the strategy of the players for different time steps of the
learning game when the learning rate is set to $0.4$ and $0.1$, respectively.
As it can be seen in these figures, decrease in the learning rate value causes
the strategy of the players to get away from the NE point of the game.
However, the convergence speed of the proposed method has not changed. It is
worth noting that the strategy of the players in gradient-based learning is
very sensitive to the learning rate, as indicated in [8]. So it can be
concluded that the proposed method is more robust compared to the gradient-
based learning method.

##  V Conclusion

In this study, we have simulated a spectrum trading problem in a cellular
network with primary users and a secondary user (a cognitive D2D pair). In the
considered problem, the primary users can share their spectrum with the D2D
pair. They are also able to resell their bandwidth to the D2D pair. Since each
user tries to optimize its own benefit, the problem can be modeled by a
bandwidth-auction game. Owing to being practically operative, the game has
been assumed to be an incomplete information repeated one as no user is aware
of the others’ strategy in a practical scenario. We have addressed the problem
by our proposed best response based learning approach. The NE point of the
game in an complete information scenario is derived and it has been
demonstrated that the proposed method converges to the NE point of the game in
an iterative manner. The simulation results have confirmed that the proposed
approach outperforms the state-of-the-art methods in terms of convergence
speed. It has been shown that the the proposed learning method can be
considered as a robust one regarding the learning rate when it is compared
with the state-of-the-art methods. Thus, the proposed best response learning
based method is an effective approach in terms of performance and convergence
speed.

## References

  * [1] S. Haykin _et al._ , “Cognitive radio: brain-empowered wireless communications,” _IEEE journal on selected areas in communications_ , vol. 23, no. 2, pp. 201–220, 2005.
  * [2] P. Mach, Z. Becvar, and T. Vanek, “In-band device-to-device communication in ofdma cellular networks: A survey and challenges,” _IEEE Communications Surveys & Tutorials_, vol. 17, no. 4, pp. 1885–1922, 2015.
  * [3] D. Feng, L. Lu, Y. Yuan-Wu, G. Y. Li, G. Feng, and S. Li, “Device-to-device communications underlaying cellular networks,” _IEEE Transactions on Communications_ , vol. 61, no. 8, pp. 3541–3551, 2013.
  * [4] M. H. Bahonar and M. J. Omidi, “Centralized qos-aware resource allocation for d2d communications with multiple d2d pairs in one resource block,” in _Electrical Engineering (ICEE), Iranian Conference on_. IEEE, 2018, pp. 643–648.
  * [5] M. Banagar, B. Maham, P. Popovski, and F. Pantisano, “Power distribution of device-to-device communications in underlaid cellular networks,” _IEEE Wireless Communications Letters_ , vol. 5, no. 2, pp. 204–207, 2016.
  * [6] A. H. Sakr and E. Hossain, “Cognitive and energy harvesting-based d2d communication in cellular networks: Stochastic geometry modeling and analysis,” _IEEE Transactions on Communications_ , vol. 63, no. 5, pp. 1867–1880, 2015.
  * [7] F. Pantisano, M. Bennis, W. Saad, and M. Debbah, “Spectrum leasing as an incentive towards uplink macrocell and femtocell cooperation,” _IEEE Journal on Selected Areas in Communications_ , vol. 30, no. 3, pp. 617–630, 2012\.
  * [8] H. Kebriaei, B. Maham, and D. Niyato, “Double-sided bandwidth-auction game for cognitive device-to-device communication in cellular networks,” _IEEE Transactions on Vehicular Technology_ , vol. 65, no. 9, pp. 7476–7487, 2016.
  * [9] L. Song, D. Niyato, Z. Han, and E. Hossain, “Game-theoretic resource allocation methods for device-to-device communication,” _IEEE Wireless Communications_ , vol. 21, no. 3, pp. 136–144, 2014.
  * [10] K. Akkarajitsakul, P. Phunchongharn, E. Hossain, and V. K. Bhargava, “Mode selection for energy-efficient d2d communications in lte-advanced networks: A coalitional game approach,” in _Communication Systems (ICCS), 2012 IEEE International Conference on_. IEEE, 2012, pp. 488–492.
  * [11] K. Zhu and E. Hossain, “Joint mode selection and spectrum partitioning for device-to-device communication: A dynamic stackelberg game,” _IEEE Transactions on Wireless Communications_ , vol. 14, no. 3, pp. 1406–1420, 2015\.
  * [12] Z. Han, R. Zheng, and H. V. Poor, “Repeated auctions with bayesian nonparametric learning for spectrum access in cognitive radio networks,” _IEEE Transactions on Wireless Communications_ , vol. 10, no. 3, pp. 890–900, 2011.
  * [13] C. Xu, L. Song, Z. Han, Q. Zhao, X. Wang, X. Cheng, and B. Jiao, “Efficiency resource allocation for device-to-device underlay communication systems: A reverse iterative combinatorial auction based approach,” _IEEE Journal on Selected Areas in Communications_ , vol. 31, no. 9, pp. 348–358, 2013.
  * [14] D. Niyato and E. Hossain, “Competitive pricing for spectrum sharing in cognitive radio networks: Dynamic game, inefficiency of nash equilibrium, and collusion,” _IEEE journal on selected areas in communications_ , vol. 26, no. 1, 2008.
  * [15] H. Sartono, Y. H. Chew, W. H. Chin, and C. Yuen, “Joint demand and supply auction pricing strategy in dynamic spectrum sharing,” in _Personal, Indoor and Mobile Radio Communications, 2009 IEEE 20th International Symposium on_. IEEE, 2009, pp. 833–837.